1 Hiyerarşik Lineer Modelleme

Vakalar için toplanan verilerin birden fazla düzeyde düzenlendiği araştırma desenlerinde çok düzeyli (hiyerarşik) doğrusal modelleme kullanılır (Tabachnick & Fidell, 2020).

library(lme4)
library(lmerTest)
library(performance)
library(ggplot2)
library(dplyr)
library(lattice)

set.seed(123)

Veri seti 750 gözlemden oluşmaktadır. Veri seti 30 okul, her okulda 25 öğrenci olacak biçimde simule edilmiştir.

# Veri üretimi
n_schools  <- 30   # grup sayısı (Seviye-2)
n_per      <- 25   # grup başına gözlem (Seviye-1)

school_id    <- rep(1:n_schools, each = n_per)
teacher_exp  <- rep(round(rnorm(n_schools, 10, 3)), each = n_per)  # Svy-2
ses          <- rnorm(n_schools * n_per, 0, 1)                      # Svy-1
cinsiyet     <- rbinom(n_schools * n_per, 1, 0.5)                   # Svy-1

u0j <- rep(rnorm(n_schools, 0, 5), each = n_per)   # rastgele kesişim
u1j <- rep(rnorm(n_schools, 0, 1.5), each = n_per) # rastgele eğim

Tüm değişkenler sıfır olsa bile her öğrencinin başladığı temel puanı 65 olarak alınmıştır.

4 × ses: Bir öğrencinin sosyoekonomik statüsü 1 birim arttığında matematik puanı 4 puan artıyor. Analizde bunu 4.44 olarak bulduk. 2 × cinsiyet: Kız öğrenciler erkeklere kıyasla başlangıçta 2 puan daha yüksek. Analizde bunu 2.46 olarak bulduk. 0.7 × teacher_exp: Bir okulun öğretmen deneyimi 1 yıl arttığında o okuldaki tüm öğrencilerin puanı 0.7 puan artıyor. Analizde bunu 1.09 bulduk. Gerçek değerden biraz yüksek çıktı, çünkü rastgele varyasyon tahmine karışıyor.

u0j Rastgele Kesişim: Her okulun kendi genel ortalaması 65’ten yukarı veya aşağı sapıyor. Bu sapma N(0, 5) dağılımından çekildi. Analizde okul varyansını 25.69 (standart sapma 5.07) bulduk, gerçek değerle birebir örtüşüyor. u1j×ses Rastgele Eğim: SES’in etkisi her okulda aynı değil. Bazı okullarda SES etkisi 4’ün çok üzerinde, bazılarında çok altında. Bu sapma N(0, 1.5) dağılımından çekildi. Analizde SES eğimi varyansını 1.44 (standart sapma 1.20) bulduk. Gerçek değerle örtüşüyor.

N(0, 7) Bireysel Hata: Aynı okulda, aynı SES ve cinsiyetteki iki öğrenci bile birbirinden farklı puan alır. Bu bireysel ölçülemez farklılıkları temsil eder. Analizde bireysel varyansı 51.31 (standart sapma 7.16) bulduk. Gerçek standart sapma olan 7 ile neredeyse aynı.

math <- 65 + 4*ses + 2*cinsiyet + 0.7*teacher_exp + u0j + u1j*ses +
        rnorm(n_schools * n_per, 0, 7)

df <- data.frame(school_id = factor(school_id),
                 teacher_exp, ses, cinsiyet, math)
head(df)
##   school_id teacher_exp        ses cinsiyet     math
## 1         1           8  0.4264642        1 81.42117
## 2         1           8 -0.2950715        0 54.27096
## 3         1           8  0.8951257        1 67.92123
## 4         1           8  0.8781335        0 68.10540
## 5         1           8  0.8215811        1 66.79001
## 6         1           8  0.6886403        0 78.08594

1.1 Tek Yönlü ANOVA-Rastgele Etkiler Modeli

# Null model: sadece rastgele kesişim, hiç bağımsız değişken yok
model0 <- lmer(math ~ 1 + (1 | school_id), data = df, REML = TRUE)
summary(model0)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ 1 + (1 | school_id)
##    Data: df
## 
## REML criterion at convergence: 5407.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.84256 -0.58849 -0.01924  0.66449  2.87074 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  school_id (Intercept) 25.69    5.068   
##  Residual              72.57    8.519   
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  72.2763     0.9762 29.0000   74.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC
# ICC > 0.05 okullar arası varyasyon anlamlı ise HLM gerekli, sonuç anlamlı
icc(model0)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.261
##   Unadjusted ICC: 0.261
# Varyans bileşenleri
VarCorr(model0)
##  Groups    Name        Std.Dev.
##  school_id (Intercept) 5.0681  
##  Residual              8.5190
# Görsel: okul bazında puan dağılımı
ggplot(df, aes(x = school_id, y = math)) +
  geom_boxplot(fill = "steelblue", alpha = 0.6) +
  labs(title = "Okullara Göre Matematik Puanı (Null Model)",
       x = "Okul", y = "Puan") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 7))

1.2 Bağımlı Değişken Olarak Ortalamalar

# Seviye-2 değişken (teacher_exp) eklenir, Seviye-1 yok
# Amaç: grup ortalamasını Svy-2 değişkenlerle açıklamak
model1 <- lmer(math ~ teacher_exp + (1 | school_id), data = df, REML = FALSE)
summary(model1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: math ~ teacher_exp + (1 | school_id)
##    Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    5405.6    5424.1   -2698.8    5397.6       746 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8973 -0.5891 -0.0268  0.6897  2.8876 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  school_id (Intercept) 15.70    3.962   
##  Residual              72.57    8.519   
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  62.0899     2.7828 30.0000  22.312  < 2e-16 ***
## teacher_exp   1.0394     0.2724 30.0000   3.816 0.000631 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## teacher_exp -0.959
# Model karşılaştırması
anova(model0, model1)
## Data: df
## Models:
## model0: math ~ 1 + (1 | school_id)
## model1: math ~ teacher_exp + (1 | school_id)
##        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## model0    3 5415.4 5429.3 -2704.7    5409.4                         
## model1    4 5405.6 5424.1 -2698.8    5397.6 11.873  1  0.0005696 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Görsel: öğretmen deneyimi-okul ortalaması ilişkisi
df_school <- df %>%
  group_by(school_id, teacher_exp) %>%
  summarise(mean_math = mean(math), .groups = "drop")

ggplot(df_school, aes(x = teacher_exp, y = mean_math)) +
  geom_point(color = "coral", size = 3) +
  geom_smooth(method = "lm", se = TRUE, color = "darkred") +
  labs(title = "Means-as-Outcomes: Öğretmen Deneyimi ~ Okul Ortalaması",
       x = "Öğretmen Deneyimi (yıl)", y = "Okul Matematik Ortalaması") +
  theme_minimal()

1.3 Rastgele Katsayı Modeli

# Svy-1 değişkenler eklenir+SES'in eğimi okuldan okula farklılaşıyor
model2 <- lmer(math ~ ses + cinsiyet + teacher_exp + (1 + ses | school_id),
               data = df, REML = FALSE)
summary(model2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: math ~ ses + cinsiyet + teacher_exp + (1 + ses | school_id)
##    Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    5175.0    5212.0   -2579.5    5159.0       742 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.03817 -0.67077  0.01928  0.72719  3.02564 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  school_id (Intercept) 15.76    3.970         
##            ses          1.44    1.200    0.51 
##  Residual              51.31    7.163         
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  60.2982     2.6366  31.8116  22.870  < 2e-16 ***
## ses           4.4410     0.3502  29.7022  12.682 1.62e-13 ***
## cinsiyet      2.4647     0.5341 717.5194   4.614 4.67e-06 ***
## teacher_exp   1.0871     0.2562  31.0198   4.243 0.000185 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ses    cinsyt
## ses          0.090              
## cinsiyet    -0.094  0.017       
## teacher_exp -0.951 -0.006 -0.009
# Rastgele eğimleri görselleştirme
df$fit_m2 <- fitted(model2)

ggplot(df, aes(x = ses, y = math, group = school_id, color = school_id)) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.5) +
  labs(title = "Rastgele Katsayı Modeli: SES Eğimi Okullara Göre Farklılaşıyor",
       x = "SES", y = "Matematik Puanı") +
  theme_minimal() +
  theme(legend.position = "none")

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

# Hem kesişim hem SES eğimi Svy-2 değişkenle (teacher_exp) açıklanıyor
# Cross-level interaction: ses*teacher_exp
model3 <- lmer(math ~ ses + cinsiyet + teacher_exp + ses:teacher_exp +
                 (1 + ses | school_id), data = df, REML = FALSE)
summary(model3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: math ~ ses + cinsiyet + teacher_exp + ses:teacher_exp + (1 +  
##     ses | school_id)
##    Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    5175.6    5217.2   -2578.8    5157.6       741 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04020 -0.64941  0.03313  0.72492  3.03034 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  school_id (Intercept) 15.722   3.965         
##            ses          1.306   1.143    0.52 
##  Residual              51.300   7.162         
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##                 Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)      61.1470     2.7324  30.4632  22.378  < 2e-16 ***
## ses               5.8014     1.2037  29.7041   4.819 3.97e-05 ***
## cinsiyet          2.4751     0.5340 717.7921   4.635 4.23e-06 ***
## teacher_exp       1.0005     0.2665  30.0374   3.755 0.000745 ***
## ses:teacher_exp  -0.1391     0.1177  29.7060  -1.182 0.246535    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ses    cinsyt tchr_x
## ses          0.280                     
## cinsiyet    -0.085  0.024              
## teacher_exp -0.954 -0.269 -0.014       
## ses:tchr_xp -0.266 -0.958 -0.020  0.279
# Etkileşimi görselleştirme
df$teacher_grp <- ifelse(df$teacher_exp < 10, "Az Deneyimli", "Deneyimli")

ggplot(df, aes(x = ses, y = math, color = teacher_grp)) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Kesişim Eğim Modeli: SES × Öğretmen Deneyimi Etkileşimi",
       x = "SES", y = "Matematik Puanı", color = "Öğretmen") +
  theme_minimal()

1.5 Shrinkage ve Birim Bazlı Katsayı Tahmini

# Rastgele etkiler: OLS tahminleri vs. EB (Empirical Bayes) tahminleri
ranef_vals <- ranef(model3)$school_id
head(ranef_vals)
##   (Intercept)         ses
## 1 -0.57686210 -0.04105376
## 2  0.98694943  0.01083050
## 3  0.01028476  0.41592699
## 4 -4.38758468  0.09107951
## 5 -3.39060897 -0.88067773
## 6  0.92250125 -0.40268289
# Okul bazında EB kesişimleri
school_effects <- data.frame(
  school_id = 1:n_schools,
  eb_intercept = ranef_vals[,1],
  eb_slope_ses = ranef_vals[,2]
)

# OLS kesişimlerini hesapla (büzülme olmadan)
ols_ints <- df %>%
  group_by(school_id) %>%
  summarise(ols_int = coef(lm(math ~ ses))[1], .groups = "drop") %>%
  mutate(school_id = as.integer(school_id))

shrink_df <- left_join(school_effects, ols_ints, by = "school_id") %>%
  mutate(grand_mean = fixef(model3)["(Intercept)"],
         eb_total = grand_mean + eb_intercept)

# Shrinkage görselleştirme
ggplot(shrink_df, aes(x = ols_int, y = eb_total)) +
  geom_point(color = "steelblue", size = 3) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
  geom_smooth(method = "lm", se = FALSE, color = "coral") +
  labs(title = "Shrinkage: OLS vs. Empirical Bayes Kesişim Tahminleri",
       x = "OLS Kesişimi", y = "EB (Büzülmüş) Kesişimi") +
  theme_minimal()

1.6 Teknik Değerlendirme ve Sabit Etkiler

# Tüm modelleri karşılaştır
anova(model0, model1, model2, model3)
## Data: df
## Models:
## model0: math ~ 1 + (1 | school_id)
## model1: math ~ teacher_exp + (1 | school_id)
## model2: math ~ ses + cinsiyet + teacher_exp + (1 + ses | school_id)
## model3: math ~ ses + cinsiyet + teacher_exp + ses:teacher_exp + (1 + ses | school_id)
##        npar    AIC    BIC  logLik -2*log(L)    Chisq Df Pr(>Chisq)    
## model0    3 5415.4 5429.3 -2704.7    5409.4                           
## model1    4 5405.6 5424.1 -2698.8    5397.6  11.8729  1  0.0005696 ***
## model2    8 5175.0 5212.0 -2579.5    5159.0 238.5598  4  < 2.2e-16 ***
## model3    9 5175.6 5217.2 -2578.8    5157.6   1.3753  1  0.2409108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model uyum istatistikleri
model_performance(model3)
## # Indices of model performance
## 
## AIC    |   AICc |    BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
## --------------------------------------------------------------------------
## 5177.4 | 5177.7 | 5219.0 |      0.476 |      0.301 | 0.250 | 6.991 | 7.162
# Sabit etkiler tablosu
fixef(model3)
##     (Intercept)             ses        cinsiyet     teacher_exp ses:teacher_exp 
##       61.146974        5.801368        2.475116        1.000544       -0.139132
confint(model3, method = "Wald")  # Güven aralıkları
##                      2.5 %      97.5 %
## .sig01                  NA          NA
## .sig02                  NA          NA
## .sig03                  NA          NA
## .sigma                  NA          NA
## (Intercept)     55.7915269 66.50242015
## ses              3.4421042  8.16063274
## cinsiyet         1.4285500  3.52168281
## teacher_exp      0.4782528  1.52283465
## ses:teacher_exp -0.3698142  0.09155028
# R² (conditional ve marginal)
r2(model3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.476
##      Marginal R2: 0.301

1.7 Sayıltıların İncelenmesi

# Düzey-1 artık normalliği
residuals_L1 <- residuals(model3)
qqnorm(residuals_L1, main = "QQ Plot: Düzey-1 Artıklar")
qqline(residuals_L1, col = "red")

hist(residuals_L1, breaks = 30, col = "lightblue",
     main = "Düzey-1 Artık Dağılımı", xlab = "Artık")

# Düzey-2 rastgele etki normalliği
re <- ranef(model3)$school_id
qqnorm(re[,1], main = "QQ Plot: Rastgele Kesişimler (Düzey-2)")
qqline(re[,1], col = "red")

qqnorm(re[,2], main = "QQ Plot: Rastgele Eğimler (SES)")
qqline(re[,2], col = "red")

# Artık vs. Uyum değerleri (Homoskedastisite)
plot(fitted(model3), residuals(model3),
     xlab = "Uyum Değerleri", ylab = "Artıklar",
     main = "Artık ~ Uyum: Varyans Homojenliği")
abline(h = 0, col = "red", lty = 2)

# Leverage ve etki analizi
plot(model3)          # lme4'ün dahili tanısal grafikleri

dotplot(ranef(model3, condVar = TRUE))  # Katman güven aralıkları
## $school_id

## Yorumlar

summary(model0)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: math ~ 1 + (1 | school_id)
##    Data: df
## 
## REML criterion at convergence: 5407.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.84256 -0.58849 -0.01924  0.66449  2.87074 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  school_id (Intercept) 25.69    5.068   
##  Residual              72.57    8.519   
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  72.2763     0.9762 29.0000   74.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
icc(model0)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.261
##   Unadjusted ICC: 0.261
VarCorr(model0)
##  Groups    Name        Std.Dev.
##  school_id (Intercept) 5.0681  
##  Residual              8.5190

Matematik puanındaki toplam varyasyonun 25.69 puanlık kısmı okuldan, 72.57 puanlık kısmı bireyden kaynaklanıyor.

ICC değeri güçlüdür. Matematik puanındaki toplam varyasyonun %26.1’i okullar arasındaki farklılıklardan kaynaklanmaktadır. Yani hangi okula gittiğin, puanını ciddi ölçüde belirliyor.

Tüm okullar ve öğrenciler bir arada düşünüldüğünde ortalama matematik puanı 72.28’dir. Bu, sonraki modellerde referans noktamızdır.

Medyanın sıfıra çok yakın olması ve Min/Max’ın simetrik görünmesi artıkların yaklaşık normal dağıldığına işaret ediyor. Sayıltılar açısından olumlu bir işaret olarak düşünülebilir.

summary(model1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: math ~ teacher_exp + (1 | school_id)
##    Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    5405.6    5424.1   -2698.8    5397.6       746 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8973 -0.5891 -0.0268  0.6897  2.8876 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  school_id (Intercept) 15.70    3.962   
##  Residual              72.57    8.519   
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  62.0899     2.7828 30.0000  22.312  < 2e-16 ***
## teacher_exp   1.0394     0.2724 30.0000   3.816 0.000631 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## teacher_exp -0.959
anova(model0, model1)
## Data: df
## Models:
## model0: math ~ 1 + (1 | school_id)
## model1: math ~ teacher_exp + (1 | school_id)
##        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## model0    3 5415.4 5429.3 -2704.7    5409.4                         
## model1    4 5405.6 5424.1 -2698.8    5397.6 11.873  1  0.0005696 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Öğretmen deneyimi sıfır yıl olsaydı beklenen ortalama matematik puanı 62.09 olurdu.

Öğretmenin deneyimi her 1 yıl arttığında, o okulun ortalama matematik puanı 1.04 puan artıyor. Bu etki istatistiksel olarak çok anlamlı (p < 0.001).

Öğretmen deneyimi, okullar arasındaki varyasyonun %38.9’unu açıklıyor. Bu güçlü bir Seviye-2 etkisi. Bireysel varyasyon hiç değişmedi çünkü bu modele henüz Seviye-1 değişken girmedik.

Model 1, null modele göre istatistiksel olarak anlamlı biçimde daha iyi uyum sağlıyor. AIC’nin düşmesi de bunu destekliyor.

Bu yüksek korelasyon bir çoklu doğrusallık uyarısı değil, kesişimin yorumlandığı nokta (deneyim = 0) verinin uzağında olduğu için oluşan matematiksel bir yapı. İlerleyen adımlarda teacher_exp’i ortalamadan sapmaya çevirirsek bu korelasyon düşer ve kesişim daha anlamlı yorumlanabilir.

summary(model2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: math ~ ses + cinsiyet + teacher_exp + (1 + ses | school_id)
##    Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    5175.0    5212.0   -2579.5    5159.0       742 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.03817 -0.67077  0.01928  0.72719  3.02564 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  school_id (Intercept) 15.76    3.970         
##            ses          1.44    1.200    0.51 
##  Residual              51.31    7.163         
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  60.2982     2.6366  31.8116  22.870  < 2e-16 ***
## ses           4.4410     0.3502  29.7022  12.682 1.62e-13 ***
## cinsiyet      2.4647     0.5341 717.5194   4.614 4.67e-06 ***
## teacher_exp   1.0871     0.2562  31.0198   4.243 0.000185 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ses    cinsyt
## ses          0.090              
## cinsiyet    -0.094  0.017       
## teacher_exp -0.951 -0.006 -0.009
anova(model1, model2)
## Data: df
## Models:
## model1: math ~ teacher_exp + (1 | school_id)
## model2: math ~ ses + cinsiyet + teacher_exp + (1 + ses | school_id)
##        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## model1    4 5405.6 5424.1 -2698.8    5397.6                         
## model2    8 5175.0 5212.0 -2579.5    5159.0 238.56  4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ses (4.44), SES 1 standart sapma arttığında matematik puanı 4.44 puan artıyor. Seviye-1’in en güçlü yordayıcısı. cinsiyet (2.46), kız öğrenciler (1), erkek öğrencilere (0) kıyasla ortalamada 2.46 puan daha yüksek alıyor. teacher_exp (1.09), Model 1’deki etkiyle tutarlı, Seviye-1 değişkenler eklendikten sonra da anlamlılığını koruyor. Güçlü bir Seviye-2 etkisi.

Rastgele kesişim (15.76) Okullar, ortalama puanları bakımından birbirinden farklılaşmaya devam ediyor. Rastgele eğim-ses (1.44) SES’in matematik üzerindeki etkisi okuldan okula anlamlı biçimde değişiyor. Standart sapma 1.20 yani bazı okullarda SES etkisi 4.44 + 1.20 = 5.64, bazılarında 4.44 − 1.20 = 3.24 civarında. SES’in eğimi sabit değil, okula göre farklılaşıyor. Korelasyon (0.51) Ortalama puanı yüksek olan okullar, SES etkisini de daha güçlü yaşıyor. Yani iyi okullar hem genel olarak daha başarılı hem de SES’e daha duyarlı.

Seviye-1 değişkenler (ses + cinsiyet) bireysel varyasyonun %29.3’ünü açıkladı. Okul varyansı Model 1’e göre neredeyse değişmedi çünkü Seviye-2 değişken (teacher_exp) zaten modelde bulunuyor.

Model 2, Model 1’e göre çok büyük bir iyileşme sağladı. Hem SES ve cinsiyet değişkenleri hem de rastgele eğim modele ciddi katkı yapıyor.

teacher_exp ile kesişim arasındaki yüksek korelasyon hâlâ var. Grand-mean centering uygulamak bu sorunu çözebilir.

summary(model3)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: math ~ ses + cinsiyet + teacher_exp + ses:teacher_exp + (1 +  
##     ses | school_id)
##    Data: df
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##    5175.6    5217.2   -2578.8    5157.6       741 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04020 -0.64941  0.03313  0.72492  3.03034 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  school_id (Intercept) 15.722   3.965         
##            ses          1.306   1.143    0.52 
##  Residual              51.300   7.162         
## Number of obs: 750, groups:  school_id, 30
## 
## Fixed effects:
##                 Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)      61.1470     2.7324  30.4632  22.378  < 2e-16 ***
## ses               5.8014     1.2037  29.7041   4.819 3.97e-05 ***
## cinsiyet          2.4751     0.5340 717.7921   4.635 4.23e-06 ***
## teacher_exp       1.0005     0.2665  30.0374   3.755 0.000745 ***
## ses:teacher_exp  -0.1391     0.1177  29.7060  -1.182 0.246535    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ses    cinsyt tchr_x
## ses          0.280                     
## cinsiyet    -0.085  0.024              
## teacher_exp -0.954 -0.269 -0.014       
## ses:tchr_xp -0.266 -0.958 -0.020  0.279
anova(model2, model3)
## Data: df
## Models:
## model2: math ~ ses + cinsiyet + teacher_exp + (1 + ses | school_id)
## model3: math ~ ses + cinsiyet + teacher_exp + ses:teacher_exp + (1 + ses | school_id)
##        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
## model2    8 5175.0 5212.0 -2579.5    5159.0                     
## model3    9 5175.6 5217.2 -2578.8    5157.6 1.3753  1     0.2409

SES’in etkisi, öğretmen deneyimine göre farklılaşıyor mu? p = 0.247 > 0.05 etkileşim istatistiksel olarak anlamlı değil.Yani öğretmen deneyimi, SES’in matematik puanı üzerindeki eğimini anlamlı biçimde değiştirmiyor. SES etkisi okuldan okula farklılaşıyor ama bu farklılaşmanın nedeni öğretmen deneyimi değil.

Etkileşim modele girince ses katsayısı 4.44’ten 5.80’e yükseldi. Bu, etkileşim teriminin ses ile yüksek korelasyonundan kaynaklanıyor. Diğer değişkenler anlamlılığını koruyor.

Etkileşim terimi varyans bileşenlerini neredeyse hiç değiştirmedi.

Model 3, Model 2’ye göre anlamlı bir iyileşme sağlamıyor. Hem AIC hem BIC aynı kararı veriyor, etkileşim terimi modelde kalmamalı.

ses:teacher_exp ile ses arasındaki -0.958 korelasyonu ciddi çoklu doğrusallığa işaret ediyor. Bu da etkileşimin anlamlı çıkmamasının bir nedeni olabilir. Grand-mean centering uygulanırsa bu korelasyon önemli ölçüde düşer.

shrink_df
##    school_id eb_intercept eb_slope_ses  ols_int grand_mean eb_total
## 1          1  -0.57686210  -0.04105376 69.40768   61.14697 60.57011
## 2          2   0.98694943   0.01083050 72.57197   61.14697 62.13392
## 3          3   0.01028476   0.41592699 77.21573   61.14697 61.15726
## 4          4  -4.38758468   0.09107951 67.78313   61.14697 56.75939
## 5          5  -3.39060897  -0.88067773 68.64071   61.14697 57.75636
## 6          6   0.92250125  -0.40268289 79.12361   61.14697 62.06947
## 7          7  -1.38889219   0.29256448 71.56417   61.14697 59.75808
## 8          8   1.94021336   1.16067898 70.87395   61.14697 63.08719
## 9          9   0.70377014  -0.07823065 71.05700   61.14697 61.85074
## 10        10  -0.73776889  -0.17813212 71.00591   61.14697 60.40920
## 11        11  -9.65540613  -1.40288224 65.68177   61.14697 51.49157
## 12        12   2.98458883   0.50551837 76.74230   61.14697 64.13156
## 13        13   4.67224422   0.87312177 78.47978   61.14697 65.81922
## 14        14  -5.16532101  -1.66744333 67.26604   61.14697 55.98165
## 15        15   7.11451394   0.70870920 78.41171   61.14697 68.26149
## 16        16   1.70420708   1.38455350 78.64268   61.14697 62.85118
## 17        17   7.49035442   0.12790045 81.98674   61.14697 68.63733
## 18        18  -3.33678360  -0.53915021 61.94617   61.14697 57.81019
## 19        19   1.53994142  -0.41284036 76.64140   61.14697 62.68691
## 20        20   1.05203989   0.51260232 72.88432   61.14697 62.19901
## 21        21   3.50744931   1.69936401 72.36213   61.14697 64.65442
## 22        22  -1.15377839  -0.62324938 69.98956   61.14697 59.99320
## 23        23   2.02993980   0.44426351 71.63384   61.14697 63.17691
## 24        24  -6.98670320  -1.05315211 62.19387   61.14697 54.16027
## 25        25  -2.33821846  -0.50120492 67.76365   61.14697 58.80876
## 26        26   1.93700689  -0.06825036 69.98199   61.14697 63.08398
## 27        27  -1.73472792  -0.84091273 73.18976   61.14697 59.41225
## 28        28   1.22190890   0.54636574 73.47578   61.14697 62.36888
## 29        29  -2.81164191  -0.82799333 66.51071   61.14697 58.33533
## 30        30   3.84638382   0.74437678 80.82714   61.14697 64.99336

2 Genel Değerlendirme ve Öğrenme Yansıması

Model kurma sürecinde her adımın bir öncekini test ettiğini deneyimledim. Sadece sabit etkiler eklemek yetmez, rastgele etkilerin de sınanması gerekir. AIC, BIC ve likelihood ratio testi üçü birden aynı kararı verdiğinde modele güven artar. Biri farklı söylediğinde ise dikkatli yorumlama şarttır.Grand-mean centering’in neden önemli olduğunu da somut biçimde gördüm. Kesişim ile teacher_exp arasındaki -0.954 korelasyonu, merkezleme yapılmadan kurulan modellerde yorumun nasıl yanıltıcı olabileceğini gösterdi.