Öncelikle aşağıdaki kodla kullanacağımız veriyi ve daha sonra çok seviyeli modelleme için gerekli kütüphaneleri yükleyelim.
install.packages("mlmRev")
library(mlmRev)
data(Exam)
install.packages("lme4")
library(lme4)
Bu verideki değişkenler:
normexam: sınav sonucu (normalize edilmiş)
school = okul kodu
standLRT = LRT sınav sonucu
schavg = okula ortalama kabul notu
Örneğimizde normexam bağımlı değişkenini, standLRT bağımsız değişkeni cinsinden açıklamaya çalışacağız. Bu veri setinde öğrenciler okullar bazında gruplanmıştır. Aynı okula giden öğrencilerin sınav sonuçlarının birbirinden bağımsız olduğunu varsayamayız, dolayısıyla lineer modelleme bu veri için uygun değildir. Onun yerine her bir okuldaki öğrencilerin bir grup olduğu çok seviyeli bir model kuracağız. Yani öğrenciler birinci seviye, okullar ikinci seviye olacak. İlk modelimiz:
m0 <- lmer(normexam ~ (1 | school), data=Exam)
Bu modeldee bağımsız değişken yok, sadece sınav sonuçlarının ortalamasının her okulda farklı olduğunu ifade eden modeli kurduk. Buna boş model (null model) diyeceğiz. Bu modelin matematiksel ifadesi $$normexam_{ij} = \beta_{0j} + \epsilon_{ij}$$ olur. Burada
- $normexam_{ij}$, $j$ numaralı okuldaki $i$ numaralı öğrencinin sınav notu,
- $\beta_{0j}$, $j$ numaralı okuldaki öğrencilerin ortalama sınav notu,
- $\epsilon_{ij}$, modelin $j$ numaralı okuldaki $i$ numaralı öğrencinin sınav notunu tahminde yaptığı hata
olur. Modelin katsayılarını
coefficients(m0)
kodu ile görüntüleyebilirsiniz. Gördüğünüz değerler (intercept adıyla) $\beta_{0j}$ katsayıları olacaktır.
Bir sonraki modelimiz:
m1 <- lmer(normexam ~ standLRT + (1 | school), data=Exam)
Bu model bağımlı değişkenin, standLRT değişkeni ile belirlendiğini ve her okulda öğrencilerin ortalama sınav sonuçlarının farklı olduğunu ifade ediyor. Bu modelin matematiksel ifadesi $$normexam_{ij} = \beta_{0j} + \beta_1standLRT_{ij}+\epsilon_{ij}$$ olur. Burada yukarıdakilere ek olarak
- $standLRT_{ij}$, $j$ numaralı okuldaki $i$ numaralı öğrencinin LRT sınav sonucu,
- $\beta_1$, her bir öğrenci için standLRT değişkeninin katsayısı
olur. Bu modelde her bir okul için farklı bir ortalama sınav değeri hesaplanırken, bütün okullarda standLRT için tek bir katsayı hesaplandığına dikkatinizi çekerim. Buna rassal orta değer- sabit eğim (random intercept-fixed slope) modeli denir, yani her bir okul için farklı bir orta değer ve sabit bir eğim vardır.
Son modelimiz:
m2 <- lmer(normexam ~ standLRT + (1+standLRT | school), data=Exam)
Bu modelin matematiksel ifadesi $$normexam_{ij} = \beta_{0j} + \beta_{1j}standLRT_{ij}+\epsilon_{ij}$$ olur. Burada $\beta_{1j}$, $j$ numaralı okulda standLRT değişkeninin eğimidir. Dolayısıyla bu modelde her bir okulun hem orta noktası (intercept) hem de standLRT değişkeninin katsayısı birbirinden farklıdır.
Bu kurduğumuz modelleri
anova(m0,m1,m2)
koduyla birbiriyle karşılaştırabilirsiniz. Çıktı olarak aldığınız değerlerin en sağında verilen $p$ değerleri bir sonraki modelin, bir öncekinden farkının istatistiki olarak ne kadar anlamlı olduğunu gösterir. Bu örnekte $m1$, veriyi $m0$'dan daha iyi tarif etmekte ve $m2$, $m1$'den daha iyi tarif etmekte.