R語言LME4混合效應模型研究教師的受歡迎程度
原文連結：http://tecdat.cn/?p=11724
介紹
本教程對多級回歸進行了基本介紹 。
本教程期望：
- 多級分析的基礎知識 。
- R中編碼的基礎知識。
- 安裝R軟體包
lme4
lmerTest
步驟1：設定
如果尚未安裝所有下面提到的軟體包，則可以通過指令安裝它們
install.packages("NAMEOFPACKAGE")
。
library(lme4) # for the analysis
library(haven) # to load the SPSS .sav file
library(tidyverse) # needed for data manipulation.
library(RColorBrewer) # needed for some extra colours in one of the graphs
library(lmerTest)# to get p-value estimations that are not part of the standard lme4 packages 受歡迎程度資料集包含不同班級學生的特征。本教程的主要目的是找到模型和檢驗關于這些特征與學生受歡迎程度（根據其同學）之間的關系的假設。 我們将使用.sav檔案，該檔案可以在SPSS檔案夾中找到。将資料下載下傳到工作目錄後，可以使用
read_sav()
指令将其打開 。
GitHub是一個平台，允許研究人員和開發人員共享代碼，軟體和研究成果，并在項目上進行協作。
資料清理
資料集中有一些我們不使用的變量，是以我們可以選擇将要使用的變量，并檢視前幾個觀察值。
popular2data <- select(popular2data, pupil, class, extrav, sex, texp, popular) # we select just the variables we will use
head(popular2data) # we have a look at the first 6 observations ## # A tibble: 6 x 6
## pupil class extrav sex texp popular
## <dbl> <dbl> <dbl> <dbl+lbl> <dbl> <dbl>
## 1 1 1 5 1 [girl] 24 6.3
## 2 2 1 7 0 [boy] 24 4.9
## 3 3 1 4 1 [girl] 24 5.3
## 4 4 1 3 1 [girl] 24 4.7
## 5 5 1 5 1 [girl] 24 6
## 6 6 1 4 0 [boy] 24 4.7 步驟3：繪制資料
在開始分析之前，我們可以繪制外向性和流行度之間的關系，而無需考慮資料的多級結構。
ggplot(data = popular2data,
aes(x = extrav,
y = popular))+
geom_point(size = 1.2,
alpha = .8,
position = "jitter")+# to add some random noise for plotting purposes
theme_minimal()+
labs(title = "Popularity vs. Extraversion") 
現在我們可以向該圖添加回歸線。
到目前為止，我們已經忽略了資料的嵌套多層結構。我們可以通過對不同類進行顔色編碼來顯示這種多級結構。
現在我們可以為資料中的100個不同類别繪制不同的回歸線
我們清楚地看到，外向性和受歡迎程度之間的關系在所有階層中并不相同，但平均而言，存在明顯的正向關系。在本教程中，我們将顯示這些不同斜率的估計值（以及如何解釋這些差異）。
我們還可以對最極端的回歸線進行顔色編碼。
現在我們可以在人氣資料上使用此功能。
f1(data = as.data.frame(popular2data),
x = "extrav",
y = "popular",
grouping = "class",
n.highest = 3,
n.lowest = 3) %>%
ggplot()+
geom_point(aes(x = extrav,
y = popular,
fill = class,
group = class),
size = 1,
alpha = .5,
position = "jitter",
shape = 21,
col = "white")+
geom_smooth(aes(x = extrav,
y = popular,
col = high_and_low,
group = class,
size = as.factor(high_and_low),
alpha = as.factor(high_and_low)),
method = lm,
se = FALSE)+
theme_minimal()+
theme(legend.position = "none")+
scale_fill_gradientn(colours = rainbow(100))+
scale_color_manual(values=c("top" = "blue",
"bottom" = "red",
"none" = "grey40"))+
scale_size_manual(values=c("top" = 1.2,
"bottom" = 1.2,
"none" = .5))+
scale_alpha_manual(values=c("top" = 1,
"bottom" = 1,
"none" =.3))+
labs(title="Linear Relationship Between Popularity and Extraversion for 100 Classes",
subtitle="The 6 with the most extreme relationship have been highlighted red and blue") 步驟4：分析資料
僅截距模型
我們複制的第一個模型是僅截距模型。
如果我們檢視LMER函數的不同輸入，則：
- “流行”，表示我們要預測的因變量。
- 一個“〜”，用于表示我們現在給出了其他感興趣的變量。（與回歸方程式的'='相比）。
- 公式中表示截距的“ 1”。
- 由于這是僅截距模式，是以我們在這裡沒有任何其他自變量。
- 在方括号之間，我們具有随機效果/斜率。同樣，值1表示垂直“ |”的截距和變量右側 條用于訓示分組變量。在這種情況下，類ID。是以，因變量“流行”是由截距和該截距的随機誤差項預測的。
- 最後，我們在
data =
summary(interceptonlymodel) #to get paramater estimates. ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Data: popular2data
##
## REML criterion at convergence: 6330.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5655 -0.6975 0.0020 0.6758 3.3175
##
## Random effects:
## Groups Name Variance Std.Dev.
## class (Intercept) 0.7021 0.8379
## Residual 1.2218 1.1053
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.07786 0.08739 98.90973 58.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 如果檢視彙總輸出，則在“随機效應”下可以看到，類别級别0.7021上的殘差和第一級别（學生級别）上的殘差為1.2218。這意味着類内相關性（ICC）為0.7021 /（1.2218 + 0.7021）=。36。
在“固定效果”下，報告截距的估計值為5.078。
我們還可以輸出計算ICC。
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.365
## Conditional ICC: 0.365 一級預測變量
現在我們可以首先添加第一（學生）水準的預測變量。一級預測因子是性别和外向性。現在，我們僅将它們添加為固定效果，而不添加為随機斜率。在此之前，我們可以繪制兩種性别在效果上的差異。我們發現性别之間可能存在平均差異，但斜率（回歸系數）沒有差異。
summary(model1) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Data: popular2data
##
## REML criterion at convergence: 4948.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2091 -0.6575 -0.0044 0.6732 2.9755
##
## Random effects:
## Groups Name Variance Std.Dev.
## class (Intercept) 0.6272 0.7919
## Residual 0.5921 0.7695
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.141e+00 1.173e-01 3.908e+02 18.25 <2e-16 ***
## sex 1.253e+00 3.743e-02 1.927e+03 33.48 <2e-16 ***
## extrav 4.416e-01 1.616e-02 1.957e+03 27.33 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sex
## sex -0.100
## extrav -0.705 -0.085 預設情況下，lmer函數僅提供測試統計資訊和估計值，而不提供p值。但是，因為我們使用，是以
lmerTest package
确實獲得了P值。現在的截距為2.14，性别的回歸系數為1.25，外向回歸系數為0.44。在輸出的固定效果表的最後一列中，我們看到了P值，這些值表示所有回歸系數均與0顯着不同。
一級和二級預測變量
現在，我們（除了均重要的1級變量）還在第二級（教師經驗）添加了預測變量。
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Data: popular2data
##
## REML criterion at convergence: 4885
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1745 -0.6491 -0.0075 0.6705 3.0078
##
## Random effects:
## Groups Name Variance Std.Dev.
## class (Intercept) 0.2954 0.5435
## Residual 0.5920 0.7694
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.098e-01 1.700e-01 2.264e+02 4.764 3.4e-06 ***
## sex 1.254e+00 3.729e-02 1.948e+03 33.623 < 2e-16 ***
## extrav 4.544e-01 1.616e-02 1.955e+03 28.112 < 2e-16 ***
## texp 8.841e-02 8.764e-03 1.016e+02 10.087 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sex extrav
## sex -0.040
## extrav -0.589 -0.090
## texp -0.802 -0.036 0.139 結果表明，級别1和級别2變量均顯着。但是，我們尚未為任何變量添加随機斜率 。
現在，我們還可以與基礎模型相比，計算出第1級和第2級的解釋方差。
- 對于級别1，這是（1.2218 – 0.592）/1.2218 = 0.52
- 對于2級，這是（0.7021 – 0.2954）/0.7021 = 0.58
具有随機斜率的一級和二級預測器（1）
現在我們還想包括随機斜率。在表2.1的第三欄中，第1級的兩個預測變量（性别和外向性）均具有随機斜率。要在LMER中完成此操作，隻需将我們要為其添加随機斜率的變量添加到輸入的随機部分。這意味着
(1|class)
成為
(1+sex+extrav |class)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.026373
## (tol = 0.002, component 1) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Data: popular2data
##
## REML criterion at convergence: 4833.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1643 -0.6555 -0.0247 0.6711 2.9571
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## class (Intercept) 1.341820 1.15837
## sex 0.002395 0.04894 -0.39
## extrav 0.034738 0.18638 -0.88 -0.10
## Residual 0.551448 0.74260
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.585e-01 1.973e-01 1.811e+02 3.845 0.000167 ***
## sex 1.251e+00 3.694e-02 9.862e+02 33.860 < 2e-16 ***
## extrav 4.529e-01 2.464e-02 9.620e+01 18.375 < 2e-16 ***
## texp 8.952e-02 8.617e-03 1.014e+02 10.389 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sex extrav
## sex -0.062
## extrav -0.718 -0.066
## texp -0.684 -0.039 0.089
## convergence code: 0
## Model failed to converge with max|grad| = 0.026373 (tol = 0.002, component 1) 我們可以看到所有固定的回歸斜率仍然很顯着。然而，沒有給出對随機效應的顯着性檢驗，但是我們确實看到，可變性别的斜率的誤差項（方差）估計很小（0.0024）。這可能意味着類别之間的SEX變量沒有斜率變化，是以可以從下一次分析中删除随機斜率估計。由于沒有針對此方差的直接顯着性檢驗，我們可以使用 軟體包的
ranova()
功能
lmerTest
，這将為我們提供類似于ANOVA的随機效果表。它檢查如果删除了某種随機效應（正式稱為似然比檢驗），則模型是否變得明顯更差，如果不是這種情況，則随機效應不顯着。
ranova(model3) ## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## npar logLik AIC LRT Df
## <none> 11 -2416.6 4855.3
## sex in (1 + sex + extrav | class) 8 -2417.4 4850.8 1.513 3
## extrav in (1 + sex + extrav | class) 8 -2441.9 4899.8 50.506 3
## Pr(>Chisq)
## <none>
## sex in (1 + sex + extrav | class) 0.6792
## extrav in (1 + sex + extrav | class) 6.232e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 我們看到性别的随機影響确實不顯着（P = 0.6792），性外向的随機影響也很顯着（P <.0001）。
具有随機斜率的一級和二級預測器
我們在忽略性别的随機斜率之後繼續。
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Data: popular2data
##
## REML criterion at convergence: 4834.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1768 -0.6475 -0.0235 0.6648 2.9684
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## class (Intercept) 1.30296 1.1415
## extrav 0.03455 0.1859 -0.89
## Residual 0.55209 0.7430
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.362e-01 1.966e-01 1.821e+02 3.745 0.000242 ***
## sex 1.252e+00 3.657e-02 1.913e+03 34.240 < 2e-16 ***
## extrav 4.526e-01 2.461e-02 9.754e+01 18.389 < 2e-16 ***
## texp 9.098e-02 8.685e-03 1.017e+02 10.475 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sex extrav
## sex -0.031
## extrav -0.717 -0.057
## texp -0.688 -0.039 0.086 我們看到：
- 截距是0.736
- 性别的固定影響是1.252
- 老師經驗的影響是0.091
- 外向性的平均影響為0.453
- 外傾斜率的随機效應為0.035
- 一級殘差為0.552
- 第二級的殘差為1.303
具有随機斜率和跨水準互動作用的一級和二級預測
作為最後一步，我們可以在教師的經驗和外向性之間添加跨層次的互動作用（因為這具有很大的随機效應，我們也許可以解釋）。換句話說，我們要調查的是，班上外向性與受歡迎程度之間關系的差異是否可以由該班老師的老師經驗來解釋。 我們添加了Extraversion和Teacher體驗之間的 層次互動。這意味着我們必須添加TEXP作為EXTRAV系數的預測因子。外向性和教師經驗之間的跨層次互動作用術語可以通過“：”符号或乘以術語來建立。
如果将所有這些都以公式形式表示，則得到：
流行度ij =β0j+β1* genderij +β2j* extraversionij + eij流行度ij =β0j+β1* genderij +β2j* extraversionij + eij。
其中β0j=γ00+γ01∗ experiencej +u0jβ0j=γ00+γ01∗ experiencej + u0j和β2j=γ20+γ21∗ experiencej +u2jβ2j=γ20+γ21∗ experiencej + u2j
合并得到：
流行度ij =γ00+γ10∗ sexij +γ20∗ extraversionij +γ01∗經驗j +γ21∗ extraversionij ∗經驗j + u2j ∗ extraversionij + u0j + eij流行度ij =γ00+γ10∗ sexij +γ20∗ extraversionij +γ01∗經驗j +γ21∗ extraij u2j * extraversionij + u0j + eij
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Data: popular2data
##
## REML criterion at convergence: 4780.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.12872 -0.63857 -0.01129 0.67916 3.05006
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## class (Intercept) 0.478639 0.69184
## extrav 0.005409 0.07355 -0.64
## Residual 0.552769 0.74348
## Number of obs: 2000, groups: class, 100
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.210e+00 2.719e-01 1.093e+02 -4.449 2.09e-05 ***
## sex 1.241e+00 3.623e-02 1.941e+03 34.243 < 2e-16 ***
## extrav 8.036e-01 4.012e-02 7.207e+01 20.031 < 2e-16 ***
## texp 2.262e-01 1.681e-02 9.851e+01 13.458 < 2e-16 ***
## extrav:texp -2.473e-02 2.555e-03 7.199e+01 -9.679 1.15e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sex extrav texp
## sex 0.002
## extrav -0.867 -0.065
## texp -0.916 -0.047 0.801
## extrav:texp 0.773 0.033 -0.901 -0.859 互動項用
extrav:texp
下标 表示，
Fixed effects
并估計為-0.025。
從這些結果中，我們現在還可以通過使用教師經驗作為第二級變量來計算解釋的外傾斜率方差：（0.03455-0.005409）/0.03455 = .843（這些結果與本書和HLM略有不同，即因為使用了不同的估算和舍入方法）。是以，外傾斜率回歸系數的方差的84.3％可以由老師的經驗來解釋。
外推系數在受歡迎程度上的截距和斜率均受教師經驗的影響。一名具有0年經驗的老師的班級中，外向得分為0的男學生（SEX = 0）的預期受歡迎度為-1.2096（這些值當然是不可能的，居中是防止這些無法實作的結果的好政策）。一名類似的（男）學生，每增加1分外向度，就将獲得0.8036分，以提高其知名度。當教師經驗增加時，每年經驗的截距也增加0.226。是以，同一個沒有外向性的男學生與一個有15年經驗的老師一起上課，其預期受歡迎度得分為-1.2096 +（15 x .226）= 2.1804。教師的經驗也減輕了外向性對普及的影響。對于具有15年經驗的教師，外向度的回歸系數僅為0.8036 –（15 x .0247）= 0.4331（相比之下，具有0年經驗的教師班級為0.8036）。
我們還可以清楚地看到，多年的教師經驗既影響截距，又影響外向度的回歸系數。
最後
在本教程結束時，我們将檢查模型的殘差是否正态分布（在兩個級别上）。除了殘差是正态分布的之外，多級模型還假設，對于不同的随機效應，殘差的方差在組（類）之間是相等的。确實存在跨組的正态性和方差相等性的統計檢驗，但是本教程僅限于視覺檢查。
首先，我們可以通過比較殘差和拟合項來檢查均方差。
我們還可以使用QQ圖檢查殘差的正态性。該圖确實表明殘差是正态分布的。
現在，我們還可以檢查它是否具有100個類别的兩個随機效果（攔截）。
▍關注我們
【大資料部落】第三方資料服務提供商,提供全面的統計分析與資料挖掘咨詢服務,為客戶定制個性化的資料解決方案與行業報告等。
▍咨詢連結：http://y0.cn/teradat
▍聯系郵箱：[email protected]