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ğimTü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
# 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
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.261
## Unadjusted ICC: 0.261
## 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))# 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
## 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()# 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")# 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()# 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()## 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
## # 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
## (Intercept) ses cinsiyet teacher_exp ses:teacher_exp
## 61.146974 5.801368 2.475116 1.000544 -0.139132
## 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
## # R2 for Mixed Models
##
## Conditional R2: 0.476
## Marginal R2: 0.301
# 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")# 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")# 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)## $school_id
## Yorumlar
## 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
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.261
## Unadjusted ICC: 0.261
## 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.
## 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
## 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.
## 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
## 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.
## 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
## 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.
## 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
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.