Hiyerarşik Lineer Modelleme

HLM Uygulama — Popülerlik (Popularity) Verisi

Giriş

Sosyal ve eğitsel araştırmalarda veri yapıları genellikle hiyerarşik (yuvalanmış) bir doğaya sahiptir: öğrenciler sınıflara, sınıflar okullara, okullar bölgelere yuvalanır. Bu rapor, 100 sınıf ve 2.000 öğrenciden oluşan Popülerlik (Popularity) veri setini temel alarak, hiyerarşik lineer modellerin (HLM) karmaşık varyans yapılarını çözümlemedeki kapasitesini ele almaktadır.

Bu veride öğrenciler (Düzey-1) sınıflara (Düzey-2) yuvalanmıştır. Aynı sınıftaki öğrenciler aynı öğretmeni, aynı sınıf iklimini ve aynı akran grubunu paylaştıkları için gözlemleri birbirinden bağımsız değildir. Geleneksel regresyon bu bağımlılığı göz ardı ettiğinden standart hataları yanlı kestirir; HLM ise varyansı düzeylere ayırarak bu sorunu çözer.

Araştırma Soruları

Analizlerde ele alınacak sorular aşağıdaki gibidir:

• Sınıflar ortalama popülerlik puanlarında ne kadar farklılık gösterir?
• Öğrencinin dışadönüklük (extraversion) düzeyi ile popülerliği arasındaki ilişkinin gücü bütün sınıflarda benzer midir, yoksa sınıftan sınıfa değişir mi?
• Ortalama dışadönüklük düzeyi yüksek olan sınıfların popülerlik düzeyi de mi yüksektir?
• Öğretmen deneyimi kontrol altına alındıktan sonra sınıflar, hem ortalama popülerlik açısından hem de dışadönüklük–popülerlik ilişkisi açısından nasıl karşılaştırılır?

Değişkenlerimiz

• popülerlik : Öğrencinin akran değerlendirmesine (sosyometrik) dayalı popülerlik puanı (Düzey-1, Bağımlı değişken).
• dışadönüklük : Öğrencinin dışadönüklük (extraversion) düzeyi (Düzey-1 yordayıcısı).
• merkezli_dışadönüklük : dışadönüklük değişkeninin sınıf ortalamasına göre merkezlenmiş hali (Group-mean centered).
• ortalama_dışadönüklük : Sınıfın ortalama dışadönüklük değeri (Düzey-2 yordayıcısı, bağlamsal değişken).
• öğretmen_deneyimi : Öğretmenin yıl cinsinden deneyimi (Düzey-2 yordayıcısı).
• cinsiyet : Öğrencinin cinsiyeti (0 = Erkek, 1 = Kız) — betimsel olarak incelenecek Düzey-1 kategorik değişkeni.

---

Paketler, Veri ve Merkezleme

Analizden önce iki hazırlık adımı şarttır. (1) dışadönüklük değişkenini sınıf ortalamasına göre merkezlemek (group-mean centering). Bu sayede Düzey-1 eğimi “öğrencinin kendi sınıfının ortalamasına göre ne kadar dışadönük olduğunun etkisini ölçer ve bağlamsal etkiden (ortalama_dışadönüklük) ayrışır. (2) Sınıf ortalamasını (ortalama_dışadönüklük) ayrı bir Düzey-2 değişkeni olarak veriye eklemek; çünkü”ortalama dışadönüklük yüksek sınıflarda popülerlik de yüksek mi?” sorusunu ancak bu bağlamsal değişkenle yanıtlayabiliriz.

library(tidyverse)
library(haven)        # SPSS .sav dosyasını okumak için
library(sjmisc)       # Frekans ve betimsel istatistikler
library(car)          # Alternatif model özetleri için (S fonksiyonu)
library(lme4)         # HLM / Mixed modelleri kurmak için
library(interactions) # Etkileşim (interaction) grafikleri için
library(broom.mixed)  # Karma model çıktılarını düzenli (tidy) almak için

veri <- read_sav("popular2.sav") %>%
  rename(
    sınıf             = class,
    dışadönüklük      = extrav,
    öğretmen_deneyimi = texp,
    popülerlik        = popular,
    cinsiyet          = sex
  )

veri <- veri %>%
  mutate(cinsiyet = factor(cinsiyet, levels = c(0, 1), labels = c("Erkek", "Kız")))

# Veri Manipülasyonu: Sınıf ortalaması ve Merkezleme (Group-Mean Centering)
veri <- veri %>%
  group_by(sınıf) %>%
  mutate(ortalama_dışadönüklük = mean(dışadönüklük, na.rm = TRUE)) %>%
  ungroup() %>%                                                       # gruplamayı kaldır
  mutate(merkezli_dışadönüklük = dışadönüklük - ortalama_dışadönüklük) %>%  # sınıf ort. göre merkezle
  relocate(merkezli_dışadönüklük, .after = dışadönüklük)

---

Betimsel İstatistikler

Modelleme öncesi her değişkenin dağılımını, ölçeğini ve kayıp veri durumunu görmek; uç değer/kodlama hatası olup olmadığını saptamak ve sonraki yorumlara zemin hazırlamak için betimsel istatistikleri inceleriz. Hiyerarşik veride bunu iki düzeyde ayrı ayrı yaparız: önce öğrenci (Düzey-1), sonra sınıf (Düzey-2) düzeyi.

Öğrenci Düzeyi Kategorik Değişkenler

veri %>%
  select(cinsiyet) %>%
  sjmisc::frq()

## cinsiyet <categorical> 
## # total N=2000 valid N=2000 mean=1.51 sd=0.50
## 
## Value |    N | Raw % | Valid % | Cum. %
## ---------------------------------------
## Erkek |  989 | 49.45 |   49.45 |  49.45
## Kız   | 1011 | 50.55 |   50.55 | 100.00
## <NA>  |    0 |  0.00 |    <NA> |   <NA>

Öğrenci Düzeyi Sürekli Değişkenler

veri %>%
  select(dışadönüklük, merkezli_dışadönüklük, popülerlik) %>%
  descr(show = c("n", "NA.prc", "mean", "sd", "range"))

Sınıf (Düzey-2) Değişkenleri

Düzey-2 değişkenleri sınıf içinde sabittir (örn. bir sınıfın tek bir öğretmeni vardır). Bunları öğrenci satırları üzerinden betimlersek her sınıfı n kez sayar ve istatistikleri yanlı kestiririz. Bu nedenle önce sınıf başına tek satır içeren bir özet veri çerçevesi (veri.S) oluşturup betimlemeleri bunun üzerinden yaparız.

veri.S <-
  veri %>%
  group_by(sınıf) %>%
  summarise(n = n(),
            öğretmen_deneyimi     = unique(öğretmen_deneyimi),
            ortalama_dışadönüklük = unique(ortalama_dışadönüklük))

veri.S %>%
  select(n, öğretmen_deneyimi, ortalama_dışadönüklük) %>%
  descr(show = c("n", "NA.prc", "mean", "sd", "range"))

---

Model 1 — Tek Yönlü ANOVA: Rastgele Etkiler Modeli

Analizin ilk aşaması, hiçbir açıklayıcı değişken içermeyen “Tam Koşulsuz Model”in (Unconditional / Boş Model) kurulmasıdır. Bu model, popülerlikteki toplam varyansın ne kadarının sınıflar arası farklardan kaynaklandığını belirlemek için kullanılan referans noktasıdır. ICC’yi buradan hesaplar, daha sonra eklenecek yordayıcıların açıkladığı varyansı (R²) bu modele kıyasla yorumlarız. Eğer sınıflar arası varyans önemsizse HLM’e gerek kalmaz; bu model o kararı verdirir.

Bu model “Rastgele Etkiler Tek Yönlü ANOVA” olarak da adlandırılır ve hiçbir yordayıcı içermez. Amacımız, sınıfların popülerlik ortalamalarının birbirinden ne kadar farklılaştığını görmektir.

Matematiksel Formülasyon:

• Level-1 (Öğrenci Düzeyi): \(Y_{ij} = \beta_{0j} + r_{ij}\)
• Level-2 (Sınıf Düzeyi): \(\beta_{0j} = \gamma_{00} + u_{0j}\)

Burada \(r_{ij} \sim N(0, \sigma^2)\) ve \(u_{0j} \sim N(0, \tau_{00})\) varsayılmaktadır. lmer() içindeki (1 | sınıf) ifadesi “Rastgele Kesişim” (Random Intercept) anlamına gelir.

mod01 <- lmer(popülerlik ~ 1 + (1 | sınıf), data = veri)
print(summary(mod01), correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: popülerlik ~ 1 + (1 | sınıf)
##    Data: veri
## 
## REML criterion at convergence: 6330.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5655 -0.6975  0.0020  0.6758  3.3175 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sınıf    (Intercept) 0.7021   0.8379  
##  Residual             1.2218   1.1053  
## Number of obs: 2000, groups:  sınıf, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.07786    0.08739    58.1

mod01 %>% S()

## Linear mixed model fit by REML 
## Call: lmer(formula = popülerlik ~ 1 + (1 | sınıf), data = veri)
## 
## Estimates of Fixed Effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.07786    0.08739    58.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimates of Random Effects (Covariance Components):
##  Groups   Name        Std.Dev.
##  sınıf    (Intercept) 0.8379  
##  Residual             1.1053  
## 
## Number of obs: 2000, groups:  sınıf, 100
## 
##   logLik       df      AIC      BIC 
## -3165.25        3  6336.51  6353.31

tidy(mod01)

Güvenirlik (Reliability) Kestirimi

HLM programlarının çıktısındaki “Reliability Estimate” değeri, her sınıfın gözlenen ortalamasının o sınıfın gerçek ortalamasını ne kadar isabetli temsil ettiğini gösterir. Güvenirlik, basitçe “Gerçek Varyansın, Toplam Varyansa Oranı” olarak tanımlanır: \(\lambda_{qj} = \tau_{00} / (\tau_{00} + v_{qj})\). Hata varyansı (\(v_{qj}\)) küçüldükçe (büyük sınıf, kesin kestirim) güvenirlik 1’e yaklaşır. lmer bu değeri doğrudan vermediği için elle hesaplıyoruz.

# 1. Modelden varyans bileşenlerini çekelim
var_comp <- as.data.frame(VarCorr(mod01))
tau_00   <- var_comp$vcov[var_comp$grp == "sınıf"]    # Sınıflar arası varyans
sigma_sq <- var_comp$vcov[var_comp$grp == "Residual"] # Sınıf içi hata varyansı

# 2. Her bir sınıfın örneklem büyüklüğünü (nj) bulalım
n_j <- veri %>%
  count(sınıf) %>%
  pull(n)

# 3. Her sınıf için kesişim güvenirliğini (lambda) hesaplayalım
reliability_j <- tau_00 / (tau_00 + (sigma_sq / n_j))

# 4. "Reliability Estimate" değeri: bu güvenirliklerin ortalamasıdır
mean_reliability <- mean(reliability_j)

print(paste("Kesişim Güvenirliği (Reliability Estimate):", round(mean_reliability, 3)))

## [1] "Kesişim Güvenirliği (Reliability Estimate): 0.919"

Boş Modelle Karşılaştırma (Olabilirlik Oranı Testi)

Sınıflar arası varyansın (rastgele kesişimin) istatistiksel olarak gerekli olup olmadığını sınamak için, rastgele etki içeren modeli (mod01) yalnızca sabit terim içeren sıradan regresyonla (mod1_sl) karşılaştırırız. Anlamlı bir ki-kare, “sınıf yapısı önemlidir, HLM kurmak gereklidir” demektir.

# sadece sabit terim içeren (yuvalanmayı yok sayan) model
mod1_sl <- lm(popülerlik ~ 1, data = veri)

# karşılaştır
anova(mod01, mod1_sl)

Model 1 Yorumu

Bu modelde Düzey-1 hata terimi \(r_{ij} \sim N(0,\sigma^2)\) ve Düzey-2 hata terimi \(u_{0j} \sim N(0,\tau_{00})\) varsayıldığında eşitliklerimiz şöyledir:

\[ \text{popülerlik}_{ij} = \beta_{0j} + r_{ij} \] \[ \beta_{0j} = \gamma_{00} + u_{0j} \]

Varyans Bileşenleri (Variance Components): Çıktıdaki Random effects tablosunda sınıf (Intercept) satırı sınıflar arası varyansı (\(\tau_{00}\)), Residual satırı ise sınıf içi (öğrenciler arası) varyansı (\(\sigma^2\)) verir.

Sınıflararası Korelasyon (ICC): Varyans bileşenlerinden ICC’yi hesaplarız:

\[ \rho = \frac{\tau_{00}}{\tau_{00} + \sigma^2} \]

Bu değer, öğrenci popülerliğindeki toplam varyansın yüzde kaçının sınıflar arası farklardan kaynaklandığını gösterir. (Bu veride ICC oldukça yüksektir — popülerlik güçlü biçimde sınıfa bağlı bir olgudur.)

Olası Değerler Aralığı (Plausible Values Range): Sınıf ortalamalarının %95’inin düştüğü aralık \(\gamma_{00} \pm 1.96\sqrt{\tau_{00}}\) ile bulunur.

Güvenirlik: Yukarıda hesaplanan değer, sınıf ortalamalarının gerçek farkları ayırt etme gücünü özetler; 1’e yakın olması ortalamaların güvenilir kestirildiğini gösterir.

Model Karşılaştırması: anova çıktısındaki ki-kare anlamlıysa (p < .001), sınıflar arasında popülerlik bakımından anlamlı farklılıklar vardır ve rastgele etkiler modeli bu farklılıkları başarıyla yakalamaktadır; dolayısıyla çok düzeyli modelleme gereklidir.

---

Model 2 — Bağımlı Değişken Olarak Ortalamalarla Regresyon (Means-as-Outcomes)

Boş model sınıflar arasında varyans olduğunu kanıtladıktan sonra, bu varyansı açıklamaya başlarız. İlk adımda sınıf düzeyindeki bir bağlamsal değişkeni — ortalama dışadönüklük (ortalama_dışadönüklük) — modele katarak “sınıfların popülerlik ortalamalarındaki farkı, sınıfların ortalama dışadönüklük düzeyi açıklıyor mu?” sorusunu yanıtlarız. Bu, kesişimi (sınıf ortalamasını) bir Düzey-2 değişkeniyle yordayan modeldir.

Matematiksel Formülasyon:

• Düzey-1 Modeli: \[ \text{popülerlik}_{ij} = \beta_{0j} + r_{ij} \]
• Düzey-2 Modeli: \[ \beta_{0j} = \gamma_{00} + \gamma_{01}(\text{ortalama dışadönüklük}_j) + u_{0j} \]

Tek karma eşitlik halinde: \[ \text{popülerlik}_{ij} = \gamma_{00} + \gamma_{01}(\text{ortalama dışadönüklük}_j) + u_{0j} + r_{ij} \]

mod02 <- lmer(popülerlik ~ ortalama_dışadönüklük + (1 | sınıf), data = veri)
print(summary(mod02), correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: popülerlik ~ ortalama_dışadönüklük + (1 | sınıf)
##    Data: veri
## 
## REML criterion at convergence: 6332.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5652 -0.6979  0.0018  0.6761  3.3182 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sınıf    (Intercept) 0.7097   0.8424  
##  Residual             1.2218   1.1053  
## Number of obs: 2000, groups:  sınıf, 100
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            5.19713    0.67209   7.733
## ortalama_dışadönüklük -0.02287    0.12779  -0.179

mod02 %>% S()

## Linear mixed model fit by REML 
## Call: lmer(formula = popülerlik ~ ortalama_dışadönüklük + (1 | sınıf), data =
##            veri)
## 
## Estimates of Fixed Effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            5.19713    0.67209   7.733 1.05e-14 ***
## ortalama_dışadönüklük -0.02287    0.12779  -0.179    0.858    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimates of Random Effects (Covariance Components):
##  Groups   Name        Std.Dev.
##  sınıf    (Intercept) 0.8424  
##  Residual             1.1053  
## 
## Number of obs: 2000, groups:  sınıf, 100
## 
##   logLik       df      AIC      BIC 
## -3166.38        4  6340.76  6363.16

tidy(mod02)

# Koşullu güvenirlik (conditional reliability)
var_comp <- as.data.frame(VarCorr(mod02))
tau_00   <- var_comp$vcov[var_comp$grp == "sınıf"]
sigma_sq <- var_comp$vcov[var_comp$grp == "Residual"]

n_j <- veri %>% count(sınıf) %>% pull(n)
reliability_j   <- tau_00 / (tau_00 + (sigma_sq / n_j))
mean_reliability <- mean(reliability_j)

print(paste("Koşullu Kesişim Güvenirliği:", round(mean_reliability, 3)))

## [1] "Koşullu Kesişim Güvenirliği: 0.92"

Model 2 Yorumu

Sabit etkiler: \(\gamma_{00}\) sabit terimi, ortalama_dışadönüklük = 0 olduğunda beklenen popülerliği; \(\gamma_{01}\) ise sınıfın ortalama dışadönüklüğündeki 1 birimlik artışın ortalama popülerliğe yansımasını verir. Etkinin işareti, büyüklüğü ve p-değeri yorumlanır.

Açıklanan Varyans Oranı (\(R^2_2\)): Düzey-2 değişkeninin sınıflar arası varyansı ne kadar azalttığını gösterir:

\[ R^2_2 = \frac{\tau_{00}^{(boş)} - \tau_{00}^{(mod02)}}{\tau_{00}^{(boş)}} \]

Not : Bu veride sınıf-ortalaması dışadönüklük tek başına sınıflar arası popülerlik farkını neredeyse hiç açıklamayabilir (etki anlamsız, \(R^2_2 \approx 0\) çıkabilir). Bu, HLM örneklerinde her bağlamsal ortalamanın güçlü bir yordayıcı olmadığını gösteren bir bulgudur. İlginç biçimde, ortalama_dışadönüklük’ün etkisi ancak öğretmen deneyimi kontrol edildikten sonra (Model 4) belirginleşebilir.

Koşullu ICC: ortalama_dışadönüklük eklendikten sonra kalan sınıflar arası bağımlılığı ölçer: \(\rho_{koşullu} = \tau_{00}^{(mod02)} / (\tau_{00}^{(mod02)} + \sigma^2)\).

Koşullu Güvenirlik: Açıklayıcı değişken eklendikçe parametre (gerçek) varyansının bir kısmı açıklanır; bu nedenle koşullu güvenirlik genellikle koşulsuz modeldekinden bir miktar düşer.

---

Model 3 — Rastgele Katsayı Modeli (Random-Coefficient Model)

Buraya kadar yalnızca sınıf ortalamalarının (kesişimlerin) değiştiğini varsaydık. Şimdi asıl ilgilendiğimiz soruya geliyoruz: “Dışadönüklük ile popülerlik arasındaki ilişkinin gücü (eğim) bütün sınıflarda aynı mı, yoksa sınıftan sınıfa değişiyor mu?” Bu modelde Düzey-1 yordayıcısı merkezli_dışadönüklük’ün eğimini rastgele bırakırız ((merkezli_dışadönüklük | sınıf)). Bunu, her sınıf için ayrı bir regresyon doğrusu (100 sınıf → 100 doğru) tahmin edip bu doğruların ortalamasını, varyansını ve kesişim–eğim korelasyonunu kestirmek olarak düşünebiliriz.

Bu model şu temel soruları motive eder:

• Bu 100 regresyon denkleminin ortalaması (ortalama kesişim ve ortalama eğim) nedir?
• Regresyon denklemlerinin (özellikle eğimlerin) sınıftan sınıfa ne kadar farklılaştığı (varyansı) nedir?
• Kesişimler ile eğimler arasındaki korelasyon nedir?

Matematiksel Formülasyon:

• Düzey-1 Modeli: \[ \text{popülerlik}_{ij} = \beta_{0j} + \beta_{1j}(\text{dışadönüklük}_{ij} - \text{ortalama dışadönüklük}_j) + r_{ij} \]
• Düzey-2 Modeli: \[ \beta_{0j} = \gamma_{00} + u_{0j} \qquad \beta_{1j} = \gamma_{10} + u_{1j} \]

mod03 <- lmer(popülerlik ~ merkezli_dışadönüklük +          # Düzey-1 yordayıcı
                           (merkezli_dışadönüklük | sınıf),  # rastgele kesişim + rastgele eğim
              data = veri)
print(summary(mod03), correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## popülerlik ~ merkezli_dışadönüklük + (merkezli_dışadönüklük |  
##     sınıf)
##    Data: veri
## 
## REML criterion at convergence: 5778
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0866 -0.7351  0.0089  0.6835  3.1729 
## 
## Random effects:
##  Groups   Name                  Variance Std.Dev. Corr 
##  sınıf    (Intercept)           0.71910  0.8480        
##           merkezli_dışadönüklük 0.03118  0.1766   -0.64
##  Residual                       0.89291  0.9449        
## Number of obs: 2000, groups:  sınıf, 100
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            5.07809    0.08742   58.09
## merkezli_dışadönüklük  0.49955    0.02693   18.55

mod03 %>% S()

## Linear mixed model fit by REML 
## Call: lmer(formula = popülerlik ~ merkezli_dışadönüklük +
##            (merkezli_dışadönüklük | sınıf), data = veri)
## 
## Estimates of Fixed Effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            5.07809    0.08742   58.09   <2e-16 ***
## merkezli_dışadönüklük  0.49955    0.02693   18.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimates of Random Effects (Covariance Components):
##  Groups   Name                  Std.Dev. Corr 
##  sınıf    (Intercept)           0.8480        
##           merkezli_dışadönüklük 0.1766   -0.64
##  Residual                       0.9449        
## 
## Number of obs: 2000, groups:  sınıf, 100
## 
##   logLik       df      AIC      BIC 
## -2888.98        6  5789.96  5823.57

tidy(mod03)

Kesişim ve Eğim Güvenirlikleri

Artık iki rastgele etkimiz var: kesişim (sınıf ortalaması) ve eğim (dışadönüklük–popülerlik ilişkisi). Her ikisinin güvenirliğini ayrı hesaplarız. Kritik fark: eğim güvenirliği, sınıf içindeki yordayıcının yayılımına (kareler toplamı) bağlıdır. Dışadönüklük bakımından homojen sınıflarda eğim daha az kesin kestirilir; bu yüzden eğim güvenirliği genellikle kesişim güvenirliğinden düşüktür.

# 1. Varyans bileşenlerini çekelim
var_comp <- as.data.frame(VarCorr(mod03))

# tau_00: sınıflar arası kesişim varyansı
tau_00 <- var_comp$vcov[var_comp$var1 == "(Intercept)" & is.na(var_comp$var2)]
# tau_11: sınıflar arası eğim (dışadönüklük slope) varyansı
tau_11 <- var_comp$vcov[var_comp$var1 == "merkezli_dışadönüklük" & is.na(var_comp$var2)]
# sigma_sq: Düzey-1 (sınıf içi) hata varyansı
sigma_sq <- var_comp$vcov[var_comp$grp == "Residual"]

# 2. Her sınıf için gerekli istatistikler (nj ve dışadönüklük yayılımı)
sinif_stats <- veri %>%
  group_by(sınıf) %>%
  summarise(
    n_j  = n(),
    dd_var = var(dışadönüklük, na.rm = TRUE),
    dd_ss  = (n_j - 1) * dd_var   # X'in kareler toplamı: sum((X - mean(X))^2)
  ) %>%
  filter(n_j > 1, dd_ss > 0)      # eğim hesaplanabilir sınıflar

# 3. Tahmin hataları (v_qj) ve güvenirlikler (lambda)
sinif_stats <- sinif_stats %>%
  mutate(
    v_0j      = sigma_sq / n_j,              # kesişim tahmin hatası
    lambda_0j = tau_00 / (tau_00 + v_0j),    # kesişim güvenirliği
    v_1j      = sigma_sq / dd_ss,            # eğim tahmin hatası (dd_ss'e bölünür)
    lambda_1j = tau_11 / (tau_11 + v_1j)     # eğim güvenirliği
  )

# 4. Ortalama güvenirlikler
mean_rel_intercept <- mean(sinif_stats$lambda_0j, na.rm = TRUE)
mean_rel_slope     <- mean(sinif_stats$lambda_1j, na.rm = TRUE)

cat("Kesişim (Sınıf Ortalaması) Güvenirliği:", round(mean_rel_intercept, 3), "\n")

## Kesişim (Sınıf Ortalaması) Güvenirliği: 0.941

cat("Eğim (Dışadönüklük-Popülerlik) Güvenirliği:", round(mean_rel_slope, 3), "\n")

## Eğim (Dışadönüklük-Popülerlik) Güvenirliği: 0.419

Model 3 Yorumu

Ortalama eğim (\(\gamma_{10}\)): Tüm sınıflar genelinde, öğrencinin kendi sınıf ortalamasının üzerindeki her 1 birim dışadönüklük artışının popülerliğe ortalama etkisidir. İşareti ve anlamlılığı yorumlanır.

Eğim heterojenliği: Random effects tablosundaki merkezli_dışadönüklük varyansı (\(\tau_{11}\)) sıfırdan büyük ve anlamlıysa, dışadönüklük–popülerlik ilişkisinin gücü sınıftan sınıfa anlamlı biçimde değişmektedir. Bu, 2. araştırma sorusuna doğrudan yanıttır.

95% Olası Değer Aralıkları:

• Sınıf Ortalamaları (Kesişimler) için: \(\gamma_{00} \pm 1.96\sqrt{\tau_{00}}\)
• Dışadönüklük–Popülerlik Eğimleri için: \(\gamma_{10} \pm 1.96\sqrt{\tau_{11}}\)

Eğim aralığı çoğu sınıfta ilişkinin pozitif olduğunu, ancak bazı sınıflarda güçlü bazılarında zayıf olduğunu gösterir.

Düzey-1’de Açıklanan Varyans (\(R^2_1\)): Öğrenci düzeyi dışadönüklüğün sınıf içi popülerlik varyansını ne kadar azalttığını gösterir:

\[ R^2_1 = \frac{\sigma^2_{(boş)} - \sigma^2_{(mod03)}}{\sigma^2_{(boş)}} \]

Kesişim–Eğim Korelasyonu: Random effects tablosundaki Corr değeri, sınıf ortalaması ile eğim arasındaki ilişkiyi verir. (Bu veride negatif çıkabilir: ortalama popülerliği yüksek sınıflarda dışadönüklüğün ek katkısı zayıflama eğiliminde olabilir — tavan etkisi olarak yorumlanabilir.)

---

Model 4 — Kesişim ve Eğimlerin Bağımlı Değişken Olduğu Model

(Intercepts- and Slopes-as-Outcomes Model)

Neden bu basamak? Model 3’te eğimlerin sınıftan sınıfa değiştiğini tespit ettik. Şimdi bu değişkenliği açıklayacağız: “Neden bazı sınıflarda ortalama popülerlik daha yüksek? Neden bazı sınıflarda dışadönüklüğün etkisi daha güçlü?” Bunun için hem kesişimi hem de eğimi, Düzey-2 değişkenleri ortalama_dışadönüklük ve öğretmen_deneyimi ile modelleriz. merkezli_dışadönüklük:ortalama_dışadönüklük ve merkezli_dışadönüklük:öğretmen_deneyimi terimleri çapraz düzey etkileşimleridir ve eğimin Düzey-2 değişkenlerine göre nasıl değiştiğini gösterir.

Matematiksel Formülasyon (Düzey-2):

\[ \beta_{0j} = \gamma_{00} + \gamma_{01}(\text{ortalama dışadönüklük}_j) + \gamma_{02}(\text{öğretmen deneyimi}_j) + u_{0j} \] \[ \beta_{1j} = \gamma_{10} + \gamma_{11}(\text{ortalama dışadönüklük}_j) + \gamma_{12}(\text{öğretmen deneyimi}_j) + u_{1j} \]

mod04 <- lmer(popülerlik ~ merkezli_dışadönüklük +                          # Düzey-1 yordayıcısı
                           ortalama_dışadönüklük + öğretmen_deneyimi +      # Düzey-2 yordayıcıları
                           merkezli_dışadönüklük:ortalama_dışadönüklük +    # Çapraz düzey etkileşimi 1
                           merkezli_dışadönüklük:öğretmen_deneyimi +        # Çapraz düzey etkileşimi 2
                           (1 + merkezli_dışadönüklük | sınıf),             # Rastgele kesişim + eğim
              data = veri)
S(mod04)

## Linear mixed model fit by REML 
## Call: lmer(formula = popülerlik ~ merkezli_dışadönüklük + ortalama_dışadönüklük
##            + öğretmen_deneyimi + merkezli_dışadönüklük:ortalama_dışadönüklük +
##            merkezli_dışadönüklük:öğretmen_deneyimi + (1 + merkezli_dışadönüklük | sınıf),
##            data = veri)
## 
## Estimates of Fixed Effects:
##                                              Estimate Std. Error z value
## (Intercept)                                 -0.842301   0.891496  -0.945
## merkezli_dışadönüklük                        0.918281   0.262321   3.501
## ortalama_dışadönüklük                        0.802684   0.139663   5.747
## öğretmen_deneyimi                            0.121345   0.014603   8.309
## merkezli_dışadönüklük:ortalama_dışadönüklük -0.004989   0.041082  -0.121
## merkezli_dışadönüklük:öğretmen_deneyimi     -0.027562   0.004278  -6.442
##                                             Pr(>|z|)    
## (Intercept)                                 0.344752    
## merkezli_dışadönüklük                       0.000464 ***
## ortalama_dışadönüklük                       9.07e-09 ***
## öğretmen_deneyimi                            < 2e-16 ***
## merkezli_dışadönüklük:ortalama_dışadönüklük 0.903335    
## merkezli_dışadönüklük:öğretmen_deneyimi     1.18e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimates of Random Effects (Covariance Components):
##  Groups   Name                  Std.Dev. Corr 
##  sınıf    (Intercept)           0.64056       
##           merkezli_dışadönüklük 0.02599  -1.00
##  Residual                       0.94357       
## 
## Number of obs: 2000, groups:  sınıf, 100
## 
##   logLik       df      AIC      BIC 
## -2851.44       10  5722.89  5778.90

print(summary(mod04), correlation = FALSE)

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## popülerlik ~ merkezli_dışadönüklük + ortalama_dışadönüklük +  
##     öğretmen_deneyimi + merkezli_dışadönüklük:ortalama_dışadönüklük +  
##     merkezli_dışadönüklük:öğretmen_deneyimi + (1 + merkezli_dışadönüklük |  
##     sınıf)
##    Data: veri
## 
## REML criterion at convergence: 5702.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2851 -0.7304  0.0130  0.6902  3.1959 
## 
## Random effects:
##  Groups   Name                  Variance  Std.Dev. Corr 
##  sınıf    (Intercept)           0.4103132 0.64056       
##           merkezli_dışadönüklük 0.0006754 0.02599  -1.00
##  Residual                       0.8903219 0.94357       
## Number of obs: 2000, groups:  sınıf, 100
## 
## Fixed effects:
##                                              Estimate Std. Error t value
## (Intercept)                                 -0.842301   0.891496  -0.945
## merkezli_dışadönüklük                        0.918281   0.262321   3.501
## ortalama_dışadönüklük                        0.802684   0.139663   5.747
## öğretmen_deneyimi                            0.121345   0.014603   8.309
## merkezli_dışadönüklük:ortalama_dışadönüklük -0.004989   0.041082  -0.121
## merkezli_dışadönüklük:öğretmen_deneyimi     -0.027562   0.004278  -6.442
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

Olası uyarı: Düzey-2 değişkenleri eğim varyansının büyük kısmını açıkladığında, kalan eğim varyansı sıfıra çok yaklaşabilir.

Birim Bazlı Katsayı Tahmini ve Büzülme (Shrinkage)

Neden bu basamak? Tek bir sınıfı sadece kendi verisiyle (OLS) tahmin etmek, özellikle küçük sınıflarda gürültülü ve uç değerler üretir. HLM’in Empirical Bayes (EB) yaklaşımı, zayıf verili birimleri genel ortalamaya doğru “büzerek” (shrinkage) ve diğer sınıflardan “güç ödünç alarak” daha kararlı tahminler verir. Aşağıda küçük bir sınıf (n=16) için OLS eğimini hesaplayıp havuzlanmış (pooled) eğimle karşılaştırarak bu mekanizmanın neden gerekli olduğunu gösteriyoruz.

# Küçük bir sınıfı tek başına (OLS ile) tahmin edelim
sinif_39 <- veri %>% filter(sınıf == 39)
lm(popülerlik ~ merkezli_dışadönüklük, data = sinif_39) %>% tidy()

# Tüm veriyi tek bir regresyona sıkıştırırsak (yuvalanmayı yok sayan pooled OLS)
lm(popülerlik ~ merkezli_dışadönüklük, data = veri) %>% tidy()

Küçük sınıfın OLS eğimi, havuzlanmış eğimden hayli uzakta ve güvensizdir; EB bu değeri genel eğime doğru çeker. Bu süreçte model, ilgili birimin zayıf verisini diğer 99 sınıftan güç ödünç alarak stabilize eder. Bunun temelinde Değiştirilebilirlik (Exchangeability) varsayımı yatar: sınıf düzeyi belirleyiciler kontrol edildikten sonra sınıfların artık etkilerinin aynı dağılımdan geldiği varsayılır.

Stratejik uyarılar: (1) Model Yanlış Kurulumu: EB tahminleri Düzey-2 modelinin doğruluğuna aşırı duyarlıdır; önemli bir sınıf düzeyi değişken dışarıda kalırsa büzülme yanlı bir noktaya yönelebilir. (2) Varyans Daralması: EB tahminleri doğası gereği gerçek parametre varyansından daha az yayılım gösterir; bu “aşırı büzülme” birimler arası farkların olduğundan küçük görünmesine yol açabilir.

Model Karşılaştırması

Neden bu basamak? Daha karmaşık modelin (mod04) daha basit modele (mod03) göre veriyi anlamlı ölçüde daha iyi açıklayıp açıklamadığını sınamak için olabilirlik oranı (ki-kare) testi kullanırız. Anlamlı sonuç, eklenen Düzey-2 yordayıcılarının ve çapraz düzey etkileşimlerinin modele gerçek katkı yaptığını gösterir.

anova(mod04, mod03, test = "Chisq")

Etkileşim Grafikleri

Neden bu basamak? Çapraz düzey etkileşim katsayıları soyuttur; “öğretmen deneyimi arttıkça dışadönüklüğün etkisi nasıl değişiyor?” sorusunu görselleştirerek somutlaştırırız. Eğim doğrularının öğretmen_deneyimi düzeyine göre ayrışması, etkileşimin yönünü ve büyüklüğünü tek bakışta anlaşılır kılar.

mod04 %>% interact_plot(pred = merkezli_dışadönüklük, modx = öğretmen_deneyimi, interval = TRUE)

mod04 %>% interact_plot(pred = merkezli_dışadönüklük, modx = öğretmen_deneyimi,
                        mod2 = ortalama_dışadönüklük, interval = TRUE)

Model 4 Yorumu

Yorum, kendi çıktınızdaki değerlerle tabloların altına yazılacaktır.

Parametre kestirimleri yerine konduğunda Düzey-2 denklemleri şu biçimi alır:

\[ \beta_{0j} = \gamma_{00} + \gamma_{01}(\text{ortalama dışadönüklük}_j) + \gamma_{02}(\text{öğretmen deneyimi}_j) + u_{0j} \] \[ \beta_{1j} = \gamma_{10} + \gamma_{11}(\text{ortalama dışadönüklük}_j) + \gamma_{12}(\text{öğretmen deneyimi}_j) + u_{1j} \]

Kesişim denkleminin yorumu (sınıf ortalamasını ne belirler?):

• Sabit terim (\(\gamma_{00}\)): Referans noktası — öğretmen_deneyimi = 0, ortalama_dışadönüklük = 0 ve kendi sınıf ortalamasında bir öğrencinin (merkezli_dışadönüklük = 0) beklenen popülerliği.
• öğretmen_deneyimi etkisi (\(\gamma_{02}\)): Öğretmen deneyimindeki her 1 yıllık artışın sınıf ortalama popülerliğine etkisi. (Bu veride güçlü ve pozitif çıkması beklenir — deneyimli öğretmenlerin sınıflarında popülerlik daha yüksektir.)
• ortalama_dışadönüklük etkisi (\(\gamma_{01}\)): Bağlamsal etki. (Model 2’de anlamsızken burada öğretmen_deneyimi kontrol edildikten sonra anlamlı hale gelebilir — bu, baskılama/karıştırıcı etkisinin öğretici bir örneğidir.)

Eğim denkleminin yorumu (ilişkinin gücünü ne değiştirir?):

• Temel eğim (\(\gamma_{10}\)): Referans bir sınıfta dışadönüklüğün popülerliğe etkisi.
• merkezli_dışadönüklük:öğretmen_deneyimi etkileşimi (\(\gamma_{12}\)): En kritik bulgu. Öğretmen deneyimi arttıkça dışadönüklük–popülerlik ilişkisinin nasıl değiştiğini gösterir. (Bu veride negatif çıkması beklenir: deneyimli öğretmenlerin sınıflarında dışadönüklüğün popülerlik üzerindeki etkisi zayıflar — yani deneyimli öğretmenler, kişilik kaynaklı popülerlik farklarını eşitleyici bir rol oynar.)
• merkezli_dışadönüklük:ortalama_dışadönüklük etkileşimi (\(\gamma_{11}\)): Sınıf ortalama dışadönüklüğü, bireysel dışadönüklüğün etkisini değiştiriyor mu?

Eğim varyansının anlamlılığı: Düzey-2 değişkenleri eklendikten sonra merkezli_dışadönüklük rastgele eğim varyansı (\(\tau_{11}\)) ve p-değeri incelenir. Varyans sıfıra yaklaşıp anlamsızsa (\(p > .05\)), “eğim farklılıkları öğretmen deneyimi ve ortalama dışadönüklük tarafından yeterince açıklanmıştır; rastgele eğim artık gerekli olmayabilir” sonucuna varılır (singular fit uyarısının anlamı budur).

Açıklanan Varyans Oranları (Model 3’e kıyasla):

• Sınıf ortalamaları (kesişim) için: \[ R^2_{kesişim} = \frac{\tau_{00}^{(mod03)} - \tau_{00}^{(mod04)}}{\tau_{00}^{(mod03)}} \]
• Dışadönüklük–popülerlik eğimleri için: \[ R^2_{eğim} = \frac{\tau_{11}^{(mod03)} - \tau_{11}^{(mod04)}}{\tau_{11}^{(mod03)}} \]

Bunlar, öğretmen deneyimi ve ortalama dışadönüklüğün sınıflar arası kesişim ve eğim farklılıklarının yüzde kaçını açıkladığını verir.

---

Sayıltılar (Varsayımların İncelenmesi)

Neden bu basamak? Hiyerarşik modellerin geçerliliği, kestirimlerin kurulmasıyla bitmez; varsayımların karşılanıp karşılanmadığının titizlikle sınanmasına bağlıdır. İhlaller (normallikten sapma, varyans heterojenliği, rastgele etkilerin çarpık dağılımı) standart hataları ve çıkarımları yanlı hale getirebilir. Aşağıda Düzey-1 ve Düzey-2 varsayımlarını sistematik olarak inceliyoruz.

Temel varsayımlar:

• Düzey-1 hata terimi \(r_{ij} \sim N(0,\sigma^2)\): bağımsız, normal dağılımlı ve sabit varyanslı (homoskedastik).
• Düzey-2 rastgele etkileri \(u_{qj}\): ortalaması sıfır olan çok değişkenli normal dağılımlı ve birbirinden bağımsız.

Aşağıdaki tanılamalar için referans modelimizi kuralım:

library(performance)  # Sayıltı kontrolü (en kritik paket)

model_hlm <- lmer(popülerlik ~ merkezli_dışadönüklük + öğretmen_deneyimi +
                               (merkezli_dışadönüklük | sınıf), data = veri)

Düzey-1 Hatalarının Normalliği

Neden? Düzey-1 kalıntılarının normalliği, sabit etkilere ilişkin çıkarımların geçerliliği için gereklidir. QQ-grafiğindeki noktaların referans doğruya yakınlığı normalliği destekler.

check_normality(model_hlm) %>% plot()

Varyansların Homojenliği (Homoscedasticity)

Neden? Kalıntı varyansının yordanan değer boyunca sabit olması varsayımıdır. Huni/yelpaze biçimli bir dağılım heteroskedastisiteye işaret eder ve robust standart hata kullanımını gündeme getirir.

check_heteroscedasticity(model_hlm) %>% plot()

Rastgele Etkilerin (Düzey-2) Normalliği

Neden? Düzey-2 rastgele etkilerinin (sınıf kesişim ve eğimleri) normal dağıldığı varsayılır. Bu varsayımın ihlali, varyans bileşeni kestirimlerini ve EB tahminlerini etkiler.

check_normality(model_hlm, effects = "random") %>% plot()

## [[1]]

Rastgele Etkilerin Dağılımı (Caterpillar Plot)

Neden? Her sınıfın rastgele etkisini güven aralıklarıyla birlikte göstererek hangi sınıfların genel ortalamadan anlamlı biçimde saptığını (uç birimleri) görselleştirir.

library(lattice)
dotplot(ranef(model_hlm, condVar = TRUE))

## $sınıf

Toplu Tanılama (All-in-One Diagnostic)

Neden? performance::check_model, yukarıdaki tüm varsayımları tek bir panelde toplayarak modelin veriye uyumunun bütüncül bir resmini verir.

check_model(model_hlm)

---

Terimler Özeti

• Sınıf İçi Korelasyon (ICC): Toplam varyansın gruplar arası oranını belirleyen katsayı.
• Güvenilirlik (Reliability): Bir grubun örneklem ortalamasının, o grubun gerçek parametresini temsil etme hassasiyeti.
• Büzülme (Shrinkage): Güvenilirliği düşük olan birim tahminlerinin, model tarafından öngörülen beklenen değere çekilmesi.
• Değiştirilebilirlik (Exchangeability): Birimlerin (sınıfların) artık etkilerinin, belirli özellikler kontrol edildikten sonra benzer bir dağılım sergilediği varsayımı.
• Çapraz Düzeyli Etkileşim: Üst düzey bir değişkenin (örn. öğretmen deneyimi), alt düzeydeki bir ilişkiyi (örn. dışadönüklük–popülerlik eğimi) değiştirmesi.

ÖĞRENME GÜNLÜĞÜ

HLM başlangıçta anlaşılır gibi başlayan ancak sonrasında kafa karıştıran bir konu. Ödevi yaparkende çok zorlandım, biraz yardım da aldım. Ama galiba sonunda biraz anladım :) Problemimiz üzerinden düşünürsek; elimizde 2000 öğrenci var ama bunlar 100 sınıfa dağılmış. Aynı sınıftaki öğrenciler aynı öğretmeni ve aynı ortamı paylaştığı için birbirlerine benziyorlar; başka sınıftaki öğrencilere göre daha çok ortak özellik taşıyorlar. Sıradan regresyon ile çözsek tüm öğrencileri bağımsız kabul etmemiz gerekir. Bu kümelenmeyi göz ardı ettiği için standart hataları olduğundan küçük hesaplamış oluruz; bu da sonuçların gerçekte olduğundan daha güvenilir görünmesine yol açıyor.Bu durumda HLM, katmanlı verilerde her gruba kendi doğrusunu veren ama bu doğruları birbirine bağlayarak az veriye sahip grupları genel eğilime yaklaştıran bir yöntem olarak çözüm buluyor. Böylece hem sınıflar arası gerçek farkları yakalıyor hem de tek tek sınıflardan gelen verileri değerlendiriyor. Hala tekrar etmem gereken bir konu olduğuna inanıyorum:)