全文鏈接：http://tecdat.cn/?p=11724

最近我們被客戶要求撰寫關(guān)于混合效應(yīng)模型的研究報告，包括一些圖形和統(tǒng)計輸出。

文中本教程對多層_回歸_模型進(jìn)行了基本介紹

介紹

本教程期望：

• 多層_回歸_模型的基礎(chǔ)知識 。

• R中編碼的基礎(chǔ)知識。

• 安裝R軟件包??lme4，和??lmerTest。

步驟1：設(shè)定?

如果尚未安裝所有下面提到的軟件包，則可以通過命令安裝它們??install.packages("NAMEOFPACKAGE")。

library(lme4)?#?用于分析library(haven)?#?加載SPSS?.sav文件library(tidyverse)?#?數(shù)據(jù)處理所需。

受歡迎程度數(shù)據(jù)集包含不同班級學(xué)生的特征。本教程的主要目的是找到模型和檢驗關(guān)于這些特征與學(xué)生受歡迎程度（根據(jù)其同學(xué)）之間的關(guān)系的假設(shè)。我們將使用.sav文件，該文件可以在SPSS文件夾中找到。將數(shù)據(jù)下載到工作目錄后，可以使用read_sav()?命令將其打開? 。

GitHub是一個平臺，允許研究人員和開發(fā)人員共享代碼，軟件和研究成果，并在項目上進(jìn)行協(xié)作。

步驟2：數(shù)據(jù)清理

數(shù)據(jù)集中有一些我們不使用的變量，因此我們可以選擇將要使用的變量，并查看前幾個觀察值。

#?我們只選擇將要使用的變量head(populardata)?#?我們來看一下前6個觀察樣本##?#?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：繪制數(shù)據(jù)

在開始分析之前，我們可以繪制外向性和受歡迎程度之間的關(guān)系，而無需考慮數(shù)據(jù)的多層結(jié)構(gòu)。

ggplot(data??=?populardata, ???????aes(x?=?extrav, ???????????y?=?popular))+ ??geom_point(size?=?1.2, ?????????????alpha?=?.8, ?????????????position?=?"jitter")+#?為繪圖目的添加一些隨機(jī)噪聲??theme_minimal()

現(xiàn)在我們可以向該圖添加回歸線。

到目前為止，我們已經(jīng)忽略了數(shù)據(jù)的嵌套多層結(jié)構(gòu)。我們可以通過對不同類進(jìn)行顏色編碼來顯示這種多層結(jié)構(gòu)。

現(xiàn)在我們可以為數(shù)據(jù)中的100個不同類別繪制不同的回歸線

我們清楚地看到，外向性和受歡迎程度之間的關(guān)系在所有層級中并不相同，但平均而言，存在明顯的正向關(guān)系。在本教程中，我們將顯示這些不同斜率的估計值（以及如何解釋這些差異）。?

點擊標(biāo)題查閱往期內(nèi)容

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

左右滑動查看更多

01

02

03

04

我們還可以對最極端的回歸線進(jìn)行顏色編碼。

人氣數(shù)據(jù)：

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)+

步驟4：分析數(shù)據(jù)

截距模型

我們第一個模型是截距模型。

如果我們查看LMER函數(shù)的不同輸入，則：

1. “受歡迎程度”，表示我們要預(yù)測的因變量。

2. 一個“?”，用于表示我們現(xiàn)在給出了其他感興趣的變量。（與回歸方程式的'='相比）。

3. 公式中表示截距的“ 1”。

4. 由于這是僅截距模式，因此我們在這里沒有任何其他自變量。

5. 在方括號之間，我們具有隨機(jī)效果/斜率。同樣，值1表示垂直“ |”的截距和變量右側(cè) 條用于指示分組變量。在這種情況下，類ID。因此，因變量“受歡迎程度”是由截距和該截距的隨機(jī)誤差項預(yù)測的。

6. 最后，我們在data =?命令后指定要使用的數(shù)據(jù)集

summary(interceptonlymodel)?#得到參數(shù)估計.##?通過REML進(jìn)行線性混合模型擬合。t檢驗使用Satterthwaite的方法##?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

如果查看匯總輸出，則在“隨機(jī)效應(yīng)”下可以看到，類別層0.7021上的殘差和第一層（學(xué)生層）上的殘差為1.2218。這意味著類內(nèi)相關(guān)性（ICC）為0.7021 /（1.2218 + 0.7021）=.36。
在“固定效果”下，報告截距的估計值為5.078。
我們還可以輸出計算ICC。

##?#?Intraclass?Correlation?Coefficient##?##??????Adjusted?ICC:?0.365##???Conditional?ICC:?0.365

一層預(yù)測變量

現(xiàn)在我們可以首先添加第一層（學(xué)生）水平的預(yù)測變量。一層預(yù)測因子是性別和外向性?，F(xiàn)在，我們僅將它們添加為固定效果，而不添加為隨機(jī)斜率。在此之前，我們可以繪制兩種性別在效果上的差異。我們發(fā)現(xiàn)性別之間可能存在平均差異，但斜率（回歸系數(shù)）沒有差異。

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

現(xiàn)在的截距為2.14，性別的回歸系數(shù)為1.25，外向回歸系數(shù)為0.44。在輸出的固定效果表的最后一列中，我們看到了P值，這些值表示所有回歸系數(shù)均與0顯著不同。

一層和二層預(yù)測變量

現(xiàn)在，我們（除了重要的1層變量）還在第2層（教師經(jīng)驗）添加了預(yù)測變量。

##?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

結(jié)果表明，層1和層2變量均顯著。但是，我們尚未為任何變量添加隨機(jī)斜率 。
現(xiàn)在，我們還可以與基礎(chǔ)模型相比，計算出第1層和第2層的解釋方差。

• 對于1層，這是（1.2218 – 0.592）/1.2218 = 0.52

• 對于2層，這是（0.7021 – 0.2954）/0.7021 = 0.58

具有隨機(jī)斜率的一層和二層預(yù)測模型（1）

現(xiàn)在我們還想包括隨機(jī)斜率。第1層的兩個預(yù)測變量（性別和外向性）均具有隨機(jī)斜率。要在LMER中完成此操作，只需將隨機(jī)斜率的變量添加到輸入的隨機(jī)部分。??(1|class)變成??(1+sex+extrav |class)。

##?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)

我們可以看到所有固定的回歸斜率仍然很顯著。然而，沒有給出對隨機(jī)效應(yīng)的顯著性檢驗，但是，可變性別的斜率的誤差項（方差）估計很?。?.0024）。這可能意味著類別之間的SEX變量沒有斜率變化，因此可以從下一次分析中刪除隨機(jī)斜率估計。由于沒有針對此方差的直接顯著性檢驗，我們可以使用?軟件包的??ranova()?函數(shù)??lmerTest，提供類似于ANOVA的隨機(jī)效果表。它檢查如果刪除了某種隨機(jī)效應(yīng)（稱為似然比檢驗），則模型是否變得明顯更差，如果不是這種情況，則隨機(jī)效應(yīng)不顯著。

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

我們看到性別的隨機(jī)影響確實不顯著（P = 0.6792），外向的隨機(jī)影響也很顯著（P <.0001）。

具有隨機(jī)斜率的一層和二層預(yù)測模型

我們在忽略性別的隨機(jī)斜率之后繼續(xù)。

##?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

• 老師經(jīng)驗的影響是0.091

• 外向的平均影響為0.453

• 外向斜率的隨機(jī)效應(yīng)為0.035

• 一層殘差為0.552

• 二層的殘差為1.303

具有隨機(jī)斜率和跨水平交互作用的一層和二層預(yù)測?

作為最后一步，我們可以在教師的經(jīng)驗和外向性之間添加跨層的交互作用。換句話說，我們要調(diào)查的是，班上外向與受歡迎程度之間關(guān)系的差異是否可以由該班老師的老師經(jīng)驗來解釋。我們添加了Extraversion和Teacher體驗之間的層級交互項。這意味著我們必須添加TEXP作為EXTRAV系數(shù)的預(yù)測因子。外向性和教師經(jīng)驗之間的跨層級交互作用可以通過“：”符號或乘以符號來創(chuàng)建。
如果將所有這些都以公式形式表示，則得到：

受歡迎程度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īng)驗j +γ21? extraversionij ?經(jīng)驗j + u2j ? extraversionij + u0j + eij

受歡迎程度ij =γ00+γ10? sexij +γ20? extraversionij +γ01?經(jīng)驗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。
從這些結(jié)果中，我們現(xiàn)在還可以通過使用教師經(jīng)驗作為第二層變量來計算解釋的外向斜率方差：（0.03455-0.005409）/0.03455 = .843。因此，外向斜率回歸系數(shù)的方差的84.3％可以由老師的經(jīng)驗來解釋。
外向系數(shù)在受歡迎程度上的截距和斜率均受教師經(jīng)驗的影響。一名具有0年經(jīng)驗的老師的班級中，外向得分為0的男學(xué)生（SEX = 0）的預(yù)期受歡迎度為-1.2096。一名類似的（男）學(xué)生，每增加1分外向度，就將獲得0.8036分，以提高其知名度。當(dāng)教師經(jīng)驗增加時，每年經(jīng)驗的截距也增加0.226。因此，同一個沒有外向性的男學(xué)生與一個有15年經(jīng)驗的老師一起上課，其預(yù)期受歡迎度得分為-1.2096 +（15 x .226）= 2.1804。教師的經(jīng)驗也減輕了外向性對普及的影響。對于具有15年經(jīng)驗的教師，外向的回歸系數(shù)僅為0.8036 –（15 x .0247）= 0.4331（相比之下，具有0年經(jīng)驗的教師班級為0.8036）。
我們還可以清楚地看到，多年的教師經(jīng)驗既影響截距，又影響外向度的回歸系數(shù)。

最后

在本教程結(jié)束，我們將檢查模型的殘差是否正態(tài)分布（在兩個層級上）。除了殘差是正態(tài)分布的之外，多層模型還假設(shè)，對于不同的隨機(jī)效應(yīng)，殘差的方差在組（類）之間是相等的。確實存在跨組的正態(tài)性和方差相等性的統(tǒng)計檢驗。

首先，我們可以通過比較殘差和擬合項來檢查均方差。

我們還可以使用QQ圖檢查殘差的正態(tài)性。該圖確實表明殘差是正態(tài)分布的。

現(xiàn)在，我們還可以檢查100個班級的兩個隨機(jī)效果。同樣，可以看到符合正態(tài)分布。

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《R語言LME4混合效應(yīng)模型研究教師的受歡迎程度》。

點擊標(biāo)題查閱往期內(nèi)容

R語言貝葉斯廣義線性混合（多層次/水平/嵌套）模型GLMM、邏輯回歸分析教育留級影響因素數(shù)據(jù)
R語言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
多水平模型、分層線性模型HLM、混合效應(yīng)模型研究教師的受歡迎程度
R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實例
R語言混合線性模型、多層次模型、回歸模型分析學(xué)生平均成績GPA和可視化
R語言線性混合效應(yīng)模型（固定效應(yīng)&隨機(jī)效應(yīng)）和交互可視化3案例
R語言用lme4多層次（混合效應(yīng)）廣義線性模型（GLM），邏輯回歸分析教育留級調(diào)查數(shù)據(jù)R語言 線性混合效應(yīng)模型實戰(zhàn)案例
R語言混合效應(yīng)邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語言如何用潛類別混合效應(yīng)模型（LCMM）分析抑郁癥狀
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言建立和可視化混合效應(yīng)模型mixed effect model
R語言LME4混合效應(yīng)模型研究教師的受歡迎程度
R語言 線性混合效應(yīng)模型實戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言如何解決線性混合模型中畸形擬合(Singular fit)的問題
基于R語言的lmer混合線性回歸模型
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型
R語言分層線性模型案例
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語言用WinBUGS 軟件對學(xué)術(shù)能力測驗建立層次（分層）貝葉斯模型
SPSS中的多層（等級）線性模型Multilevel linear models研究整容手術(shù)數(shù)據(jù)
用SPSS估計HLM多層（層次）線性模型模型