使用R語言跑多層次模型:模型比較

實務中，我們很難決定要使用哪一種多層次模型更好，這時候可以透過模型比較方式給我們有用的建議。本文章將介紹嵌套模型，還有如何透過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)