散文網 » 生活 »日常 » R語言線性混合效應模型實戰(zhàn)案例

R語言線性混合效應模型實戰(zhàn)案例

2020-12-03 11:07 作者:拓端tecdat 0人讀過 | 我要投稿

原文鏈接：http://tecdat.cn/?p=3059

?

介紹

處理分組數據和復雜層次結構的分析師，從嵌入在參與者中的測量，嵌套在州內的縣或嵌套在教室內的學生，經常發(fā)現(xiàn)他們需要建模工具來反映他們數據的這種結構。在R中，有兩種主要的方法來擬合多級模型，這些模型考慮了數據中的這種結構。這些教程將向用戶展示如何使用lme4R中的包來擬合線性和非線性混合效果模型，以及如何使用rstan以完全適合貝葉斯多級模型。這里的重點是如何使模型適合R而不是模型背后的理論。有關多級建模的背景知識，請參閱參考資料。?

本教程將介紹如何lme4? 設置和運行一些基本模型，其中包括：

在R中構造變化的截距，變化的斜率以及變化的斜率和截距模型
從混合效應模型中生成預測和解釋參數
廣義和非線性多層次模型
完全貝葉斯多級模型適合rstan或其他MCMC方法

設置環(huán)境

在R中開始多級建模很簡單。lme4是在R中實現(xiàn)多級模型的規(guī)范包，盡管有許多包依賴并增強其功能集，包括貝葉斯擴展。lme4?最近已被重寫以提高速度并整合C ++代碼庫，因此封裝的功能有些不斷變化。?

要安裝lme4，我們只需運行：

# 主要版本
install.packages("lme4")
# 或安裝開發(fā)版本
library(devtools)
install_github("lme4", user = "lme4")

讀入數據

多級模型適用于特定類型的數據結構，其中單元嵌套在組內（通常為5個以上組），并且我們希望對數據的組結構進行建模。對于我們的介紹性示例，我們將從lme4文檔中的一個簡單示例開始，并解釋模型正在執(zhí)行的操作。?

library(lme4) ?# 加載庫
library(arm) ?# R中用于回歸的函數
# summary(lmm.data)
head(lmm.data)

## ? id extro ?open agree social class school
## 1 ?1 63.69 43.43 38.03 ?75.06 ? ? d ? ? IV
## 2 ?2 69.48 46.87 31.49 ?98.13 ? ? a ? ? VI
## 3 ?3 79.74 32.27 40.21 116.34 ? ? d ? ? VI
## 4 ?4 62.97 44.41 30.51 ?90.47 ? ? c ? ? IV
## 5 ?5 64.25 36.86 37.44 ?98.52 ? ? d ? ? IV
## 6 ?6 50.97 46.26 38.83 ?75.22 ? ? d ? ? ?I

?模型

讓我們首先擬合一個簡單的OLS回歸。?

OLSexamp <- lm(extro ~ open + agree + social, data = lmm.data)
display(OLSexamp)

## lm(formula = extro ~ open + agree + social, data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 57.84 ? ? 3.15
## open ? ? ? ? 0.02 ? ? 0.05
## agree ? ? ? ?0.03 ? ? 0.05
## social ? ? ? 0.01 ? ? 0.02
## ---
## n = 1200, k = 4
## residual sd = 9.34, R-Squared = 0.00

?R模型接口非常簡單，首先指定因變量，然后是?~符號，預測變量，每個都被命名。加法符號表明這些被建模為加性效應。最后，我們指定要計算模型的數據。這里我們使用該lm函數執(zhí)行OLS回歸，但R中還有許多其他選項。

如果我們想要提取諸如AIC之類的度量。

MLexamp <- glm(extro ~ open + agree + social, data = lmm.data)
display(MLexamp)

## glm(formula = extro ~ open + agree + social, data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 57.84 ? ? 3.15
## open ? ? ? ? 0.02 ? ? 0.05
## agree ? ? ? ?0.03 ? ? 0.05
## social ? ? ? 0.01 ? ? 0.02
## ---
## ? n = 1200, k = 4
## ? residual deviance = 104378.2, null deviance = 104432.7 (difference = 54.5)
## ? overdispersion parameter = 87.3
## ? residual sd is sqrt(overdispersion) = 9.34

AIC(MLexamp)## [1] 8774

這導致模型擬合較差?，F(xiàn)在讓我們看一個簡單的模型。

擬合不同的模型

我們的下一步可能是使用分組變量（如學校或班級）來擬合不同的模型。?

MLexamp.2 <- glm(extro ~ open + agree + social + class, data = lmm.data)
display(MLexamp.2)

## glm(formula = extro ~ open + agree + social + class, data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 56.05 ? ? 3.09
## open ? ? ? ? 0.03 ? ? 0.05
## agree ? ? ? -0.01 ? ? 0.05
## social ? ? ? 0.01 ? ? 0.02
## classb ? ? ? 2.06 ? ? 0.75
## classc ? ? ? 3.70 ? ? 0.75
## classd ? ? ? 5.67 ? ? 0.75
## ---
## ? n = 1200, k = 7
## ? residual deviance = 99187.7, null deviance = 104432.7 (difference = 5245.0)
## ? overdispersion parameter = 83.1
## ? residual sd is sqrt(overdispersion) = 9.12

AIC(MLexamp.2)## [1] 8719anova(MLexamp, MLexamp.2, test = "F")

## Analysis of Deviance Table
##
## Model 1: extro ~ open + agree + social
## Model 2: extro ~ open + agree + social + class
## ? Resid. Df Resid. Dev Df Deviance ? ?F ?Pr(>F)
## 1 ? ? ?1196 ? ? 104378
## 2 ? ? ?1193 ? ? ?99188 ?3 ? ? 5190 20.8 3.8e-13 ***
## ---
## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

這通常被稱為固定效應。?

MLexamp.3 <- glm(extro ~ open + agree + social + school, data = lmm.data)
display(MLexamp.3)

## glm(formula = extro ~ open + agree + social + school, data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 45.02 ? ? 0.92
## open ? ? ? ? 0.01 ? ? 0.01
## agree ? ? ? ?0.03 ? ? 0.02
## social ? ? ? 0.00 ? ? 0.00
## schoolII ? ? 7.91 ? ? 0.27
## schoolIII ? 12.12 ? ? 0.27
## schoolIV ? ?16.06 ? ? 0.27
## schoolV ? ? 20.43 ? ? 0.27
## schoolVI ? ?28.05 ? ? 0.27
## ---
## ? n = 1200, k = 9
## ? residual deviance = 8496.2, null deviance = 104432.7 (difference = 95936.5)
## ? overdispersion parameter = 7.1
## ? residual sd is sqrt(overdispersion) = 2.67

AIC(MLexamp.3)## [1] 5774anova(MLexamp, MLexamp.3, test = "F")

## Analysis of Deviance Table
##
## Model 1: extro ~ open + agree + social
## Model 2: extro ~ open + agree + social + school
## ? Resid. Df Resid. Dev Df Deviance ? ?F Pr(>F)
## 1 ? ? ?1196 ? ? 104378
## 2 ? ? ?1191 ? ? ? 8496 ?5 ? ?95882 2688 <2e-16 ***
## ---
## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

學校因子大大提高了我們的模型擬合度。但是，我們如何解釋這些影響呢？

table(lmm.data$school, lmm.data$class)

##
## ? ? ? ?a ?b ?c ?d
## ? I ? 50 50 50 50
## ? II ?50 50 50 50
## ? III 50 50 50 50
## ? IV ?50 50 50 50
## ? V ? 50 50 50 50
## ? VI ?50 50 50 50

在這里，我們可以看到我們有一個完美平衡的設計，在課堂和學校的每個組合中有50個觀察結果。

讓我們嘗試對這些獨特的因素進行建模。?

MLexamp.4 <- glm(extro ~ open + agree + social + school:class, data = lmm.data)
display(MLexamp.4)

## glm(formula = extro ~ open + agree + social + school:class, data = lmm.data)
## ? ? ? ? ? ? ? ? ?coef.est coef.se
## (Intercept) ? ? ? 80.36 ? ? 0.37
## open ? ? ? ? ? ? ? 0.01 ? ? 0.00
## agree ? ? ? ? ? ? -0.01 ? ? 0.01
## social ? ? ? ? ? ? 0.00 ? ? 0.00
## schoolI:classa ? -40.39 ? ? 0.20
## schoolII:classa ?-28.15 ? ? 0.20
## schoolIII:classa -23.58 ? ? 0.20
## schoolIV:classa ?-19.76 ? ? 0.20
## schoolV:classa ? -15.50 ? ? 0.20
## schoolVI:classa ?-10.46 ? ? 0.20
## schoolI:classb ? -34.60 ? ? 0.20
## schoolII:classb ?-26.76 ? ? 0.20
## schoolIII:classb -22.59 ? ? 0.20
## schoolIV:classb ?-18.71 ? ? 0.20
## schoolV:classb ? -14.31 ? ? 0.20
## schoolVI:classb ? -8.54 ? ? 0.20
## schoolI:classc ? -31.86 ? ? 0.20
## schoolII:classc ?-25.64 ? ? 0.20
## schoolIII:classc -21.58 ? ? 0.20
## schoolIV:classc ?-17.58 ? ? 0.20
## schoolV:classc ? -13.38 ? ? 0.20
## schoolVI:classc ? -5.58 ? ? 0.20
## schoolI:classd ? -30.00 ? ? 0.20
## schoolII:classd ?-24.57 ? ? 0.20
## schoolIII:classd -20.64 ? ? 0.20
## schoolIV:classd ?-16.60 ? ? 0.20
## schoolV:classd ? -12.04 ? ? 0.20
## ---
## ? n = 1200, k = 27
## ? residual deviance = 1135.9, null deviance = 104432.7 (difference = 103296.8)
## ? overdispersion parameter = 1.0
## ? residual sd is sqrt(overdispersion) = 0.98

AIC(MLexamp.4)## [1] 3396

這非常有用，但如果我們想了解學校的影響和課堂的影響，以及學校和班級的影響，該怎么辦？?

MLexamp.5 <- glm(extro ~ open + agree + social + school * class - 1, data = lmm.data)
display(MLexamp.5)

## glm(formula = extro ~ open + agree + social + school * class -
## ? ? 1, data = lmm.data)
## ? ? ? ? ? ? ? ? ?coef.est coef.se
## open ? ? ? ? ? ? ?0.01 ? ? 0.00
## agree ? ? ? ? ? ?-0.01 ? ? 0.01
## social ? ? ? ? ? ?0.00 ? ? 0.00
## schoolI ? ? ? ? ?39.96 ? ? 0.36
## schoolII ? ? ? ? 52.21 ? ? 0.36
## schoolIII ? ? ? ?56.78 ? ? 0.36
## schoolIV ? ? ? ? 60.60 ? ? 0.36
## schoolV ? ? ? ? ?64.86 ? ? 0.36
## schoolVI ? ? ? ? 69.90 ? ? 0.36
## classb ? ? ? ? ? ?5.79 ? ? 0.20
## classc ? ? ? ? ? ?8.53 ? ? 0.20
## classd ? ? ? ? ? 10.39 ? ? 0.20
## schoolII:classb ?-4.40 ? ? 0.28
## schoolIII:classb -4.80 ? ? 0.28
## schoolIV:classb ?-4.74 ? ? 0.28
## schoolV:classb ? -4.60 ? ? 0.28
## schoolVI:classb ?-3.87 ? ? 0.28
## schoolII:classc ?-6.02 ? ? 0.28
## schoolIII:classc -6.54 ? ? 0.28
## schoolIV:classc ?-6.36 ? ? 0.28
## schoolV:classc ? -6.41 ? ? 0.28
## schoolVI:classc ?-3.65 ? ? 0.28
## schoolII:classd ?-6.81 ? ? 0.28
## schoolIII:classd -7.45 ? ? 0.28
## schoolIV:classd ?-7.24 ? ? 0.28
## schoolV:classd ? -6.93 ? ? 0.28
## schoolVI:classd ? 0.06 ? ? 0.28
## ---
## ? n = 1200, k = 27
## ? residual deviance = 1135.9, null deviance = 4463029.9 (difference = 4461894.0)
## ? overdispersion parameter = 1.0
## ? residual sd is sqrt(overdispersion) = 0.98

AIC(MLexamp.5)## [1] 3396

探索隨機斜率

另一種選擇是為每個學校和班級組合建立一個單獨的模型。如果我們認為我們的變量之間的關系可能高度依賴于學校和班級組合，我們可以簡單地擬合一系列模型并探索它們之間的參數變化：

require(plyr)
display(modellist[[1]])

## glm(formula = extro ~ open + agree + social, data = x)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 35.87 ? ? 5.90
## open ? ? ? ? 0.05 ? ? 0.09
## agree ? ? ? ?0.02 ? ? 0.10
## social ? ? ? 0.01 ? ? 0.03
## ---
## ? n = 50, k = 4
## ? residual deviance = 500.2, null deviance = 506.2 (difference = 5.9)
## ? overdispersion parameter = 10.9
## ? residual sd is sqrt(overdispersion) = 3.30

display(modellist[[2]])

## glm(formula = extro ~ open + agree + social, data = x)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 47.96 ? ? 2.16
## open ? ? ? ?-0.01 ? ? 0.03
## agree ? ? ? -0.03 ? ? 0.03
## social ? ? ?-0.01 ? ? 0.01
## ---
## ? n = 50, k = 4
## ? residual deviance = 47.9, null deviance = 49.1 (difference = 1.2)
## ? overdispersion parameter = 1.0
## ? residual sd is sqrt(overdispersion) = 1.02

我們將在未來的教程中更深入地討論此策略，包括如何在此命令中生成的模型列表中進行性能推斷。

建立不同的斜率模型

雖然上述所有技術都是解決這一問題的有效方法，但當我們明確感興趣的是群體之間的變化時，它們并不一定是最好的方法。這是混合效果建?？蚣苡杏玫牡胤健，F(xiàn)在我們使用lmer具有熟悉的公式接口的函數，使用特殊語法指定組級變量：(1|school) ，使lmer擬合具有變量截距組效果的線性模型school。

display(MLexamp.6)

## lmer(formula = extro ~ open + agree + social + (1 | school),
## ? ? data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 59.12 ? ? 4.10
## open ? ? ? ? 0.01 ? ? 0.01
## agree ? ? ? ?0.03 ? ? 0.02
## social ? ? ? 0.00 ? ? 0.00
##
## Error terms:
## ?Groups ? Name ? ? ? ?Std.Dev.
## ?school ? (Intercept) 9.79
## ?Residual ? ? ? ? ? ? 2.67
## ---
## number of obs: 1200, groups: school, 6
## AIC = 5836.1, DIC = 5789
## deviance = 5806.5

我們可以使用多個組來擬合多個組效果。

display(MLexamp.7)

## lmer(formula = extro ~ open + agree + social + (1 | school) +
## ? ? (1 | class), data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 60.20 ? ? 4.21
## open ? ? ? ? 0.01 ? ? 0.01
## agree ? ? ? -0.01 ? ? 0.01
## social ? ? ? 0.00 ? ? 0.00
##
## Error terms:
## ?Groups ? Name ? ? ? ?Std.Dev.
## ?school ? (Intercept) 9.79
## ?class ? ?(Intercept) 2.41
## ?Residual ? ? ? ? ? ? 1.67
## ---
## number of obs: 1200, groups: school, 6; class, 4
## AIC = 4737.9, DIC = 4683
## deviance = 4703.6

最后，我們可以通過以下語法擬合嵌套：

display(MLexamp.8)

## lmer(formula = extro ~ open + agree + social + (1 | school/class),
## ? ? data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 60.24 ? ? 4.01
## open ? ? ? ? 0.01 ? ? 0.00
## agree ? ? ? -0.01 ? ? 0.01
## social ? ? ? 0.00 ? ? 0.00
##
## Error terms:
## ?Groups ? ? ? Name ? ? ? ?Std.Dev.
## ?class:school (Intercept) 2.86
## ?school ? ? ? (Intercept) 9.69
## ?Residual ? ? ? ? ? ? ? ? 0.98
## ---
## number of obs: 1200, groups: class:school, 24; school, 6
## AIC = 3568.6, DIC = 3508
## deviance = 3531.1

在這里(1|school/class)，我們想要為1|學校和學校內的課程設置不同截距的混合效應。

用lmer擬合變化的斜率模型

但是，如果我們想要探索不同學生水平指標的影響，因為它們因教室而異。我們可以擬合不同的斜率模型，而不是按學校（或學校/班級）擬合模型。在這里，我們修改我們的隨機效應項，在分組術語之前包含變量：(1 + open|school/class)告訴R擬合變化的斜率和不同的學校和學校類別的截距模型，并允許open變量的斜率因學校而異。

display(MLexamp.9)

## lmer(formula = extro ~ open + agree + social + (1 + open | school/class),
## ? ? data = lmm.data)
## ? ? ? ? ? ? coef.est coef.se
## (Intercept) 60.26 ? ? 3.93
## open ? ? ? ? 0.01 ? ? 0.01
## agree ? ? ? -0.01 ? ? 0.01
## social ? ? ? 0.00 ? ? 0.00
##
## Error terms:
## ?Groups ? ? ? Name ? ? ? ?Std.Dev. Corr
## ?class:school (Intercept) 2.62
## ? ? ? ? ? ? ? open ? ? ? ?0.01 ? ? 1.00
## ?school ? ? ? (Intercept) 9.51
## ? ? ? ? ? ? ? open ? ? ? ?0.00 ? ? 1.00
## ?Residual ? ? ? ? ? ? ? ? 0.98
## ---
## number of obs: 1200, groups: class:school, 24; school, 6
## AIC = 3574.7, DIC = 3506
## deviance = 3529.3

結論

在R語言和生態(tài)系統(tǒng)中，擬合混合效應模型和探索組變異非常容易。在以后的教程中，我們將探索模型的比較，使用混合效果模型進行推理，以及創(chuàng)建混合效果模型的圖形表示了解它們的效果。

附錄

## Platform: x86_64-w64-mingw32/x64 (64-bit)
##
## attached base packages:
## [1] stats ? ? graphics ?grDevices utils ? ? datasets ?methods ? base
##
## other attached packages:
## [1] plyr_1.8 ? ? ? ?arm_1.6-10 ? ? ?MASS_7.3-29 ? ? lme4_1.0-5
## [5] Matrix_1.1-0 ? ?lattice_0.20-24 knitr_1.5
##
## loaded via a namespace (and not attached):
## ?[1] abind_1.4-0 ? ?coda_0.16-1 ? ?evaluate_0.5.1 formatR_0.10
## ?[5] grid_3.0.1 ? ? minqa_1.2.1 ? ?nlme_3.1-113 ? splines_3.0.1
## ?[9] stringr_0.6.2 ?tools_3.0.1

?

非常感謝您閱讀本文，有任何問題請在下面留言！

參考文獻

?

1.基于R語言的lmer混合線性回歸模型

?

2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）

?

3.R語言線性混合效應模型實戰(zhàn)案例

?

4.R語言線性混合效應模型實戰(zhàn)案例2

?

5.R語言線性混合效應模型實戰(zhàn)案例

?

6.線性混合效應模型Linear Mixed-Effects Models的部分折疊Gibbs采樣

?

7.R語言LME4混合效應模型研究教師的受歡迎程度

?

8.R語言中基于混合數據抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長

?

9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM請選中你要保存的內容，粘貼到此文本框

標簽：