實務中，我們很難決定要使用哪一種多層次模型更好，這時候可以透過模型比較方式給我們有用的建議。本文章將介紹嵌套模型，還有如何透過R語言，使用ranova()和anova()來比較不同的多層次模型，並提供實際範例。
0.嵌套模型
嵌套模型比較是指不同嵌套模型進行比較，以確定哪個模型更適合數據。嵌套模型是指一個模型包含另一個模型的情況。常見的策略是，一個多層次模型中模型的將A因子設定為隨機效應，是另一個模型則將A因子設定固定效應。兩個模型比較將可以得知A因子是否的隨機效果是否成立。假設我們有一個數據集，其中包含學生的數學成績和學生的班級。我們想研究班級對學生數學成績的影響是否具有的隨機斜率效果。
我們可以建立兩個嵌套模型來研究這個問題：
- 模型 1：班級斜率作為固定效應。
- 模型 2：班級斜率作為隨機效應。
模型 1 是模型 2 的嵌套模型，因為模型 2 的隨機效應是模型 1 的固定效應。
1.加載相關package
library(lme4)
library(lmerTest)
library(ggplot2)
library(interactions)
library(haven)
2.加載資料
hsbfull<- read_sav("hsbfull with meanses and sesdev.sav")
3.跑Random Coefficients 模型
CroLev1<- lmer(mathach ~ 1 + sesdev +
( 1+sesdev|schoolid), data = hsbfull)
summary(CroLev1)
這段程式碼是使用 hsbfull 資料集。模型包含：
sesdev: 社會經濟地位 (SES)
模型也包含2個隨機效果：
(1|schoolid): 學校層級的隨機截距(sesdev|schoolid): 學校層級的 SES 隨機斜率
模型的預測變量是學生的數學成績 (mathach)。
summary() 函式會輸出模型的結果，其中包含估計的參數、標準誤差和 p 值。
4.使用ranova()比較模型
ranova()會自動將目前模型和減少一個隨機效果的模型做比較，他會優先減少斜率再來才是截距，例如:
- 原本模型
(1 | 層次2id)，則比較模型是直接將隨機截距拿掉，變成傳統迴歸分析，1代表截距 - 原本模型
(1+因子1 | 層次2id)，則比較模型為(1+ | 層次2id)，1代表截距 - 原本模型
(1+因子1+因子2 | 層次2id)，則分別比較模型為(1+因子1 | 層次2id)和(1+因子2 | 層次2id) - 因子可以是連續或類別
本範例中，輸入ranova(CroLev1) ，結果如下:
Model:
mathach ~ sesdev + (1 + sesdev | schoolid)
npar logLik AIC LRT Df Pr(>Chisq)
<none> 6 -23357 46726
sesdev in (1 + sesdev | schoolid) 4 -23362 46732 9.7617 2 0.007591 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
npar代表該模型用的參數<none>就是原本的模型，其隨機截距+斜率效果為(1 + sesdev | schoolid)sesdev in (1 + sesdev | schoolid)表示刪除sesdev隨機效應項的模型。也就是只有隨機截距效果(1 | schoolid)- LRT為9.7617，P為0.007591 **，代表兩個模型有顯著差異，至於哪個比較好，可以從
logLik和AIC哪個比較小來判斷。
5.使用anova()比較模型
也可以用anova()語法來比較兩個模型，但兩個模型都要透過lmer()先建立一個Random_Intercepts跟Random_Coefficients 模型 。然後透過anova(random_intercepts, Random_Coefficients)進行比較兩個模型。比較結果如下:
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
random_intercepts 4 46728 46756 -23360 46720
Random_Coefficients 10 46572 46640 -23276 46552 168.76 6 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Chisq為168.76，P為2.2e-16，代表兩個模型有顯著差異，至於哪個比較好，可以從deviance/AIC/BIC哪個比較小來判斷。此做法比較自由，可以比較任何兩種模型。但也比較麻煩，要自己建立模型。
6.使用ranova()比較模型:兩隨機斜率模型
random_twoslopes<- lmer(mathach ~ 1 + sesdev + ses+
(1+ sesdev+hsbfull$ses|schoolid), data = hsbfull)
