R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度|附代碼數(shù)據(jù)
?全文鏈接：
文中本教程對(duì)多層_回歸_模型進(jìn)行了基本介紹（點(diǎn)擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)）?。??
介紹
本教程期望：
多層_回歸_模型的基礎(chǔ)知識(shí) 。
R中編碼的基礎(chǔ)知識(shí)。
安裝R軟件包??
lme4，和??lmerTest。
步驟1：設(shè)定?
如果尚未安裝所有下面提到的軟件包，則可以通過(guò)命令安裝它們??install.packages("NAMEOFPACKAGE")。
library(lme4)?#?用于分析 library(haven)?#?加載SPSS?.sav文件 library(tidyverse)?#?數(shù)據(jù)處理所需。
受歡迎程度數(shù)據(jù)集包含不同班級(jí)學(xué)生的特征。本教程的主要目的是找到模型和檢驗(yàn)關(guān)于這些特征與學(xué)生受歡迎程度（根據(jù)其同學(xué)）之間的關(guān)系的假設(shè)。我們將使用.sav文件，該文件可以在SPSS文件夾中找到。將數(shù)據(jù)下載到工作目錄后，可以使用read_sav()?命令將其打開(kāi)? 。
GitHub是一個(gè)平臺(tái)，允許研究人員和開(kāi)發(fā)人員共享代碼，軟件和研究成果，并在項(xiàng)目上進(jìn)行協(xié)作。
步驟2：數(shù)據(jù)清理
數(shù)據(jù)集中有一些我們不使用的變量，因此我們可以選擇將要使用的變量，并查看前幾個(gè)觀察值。
#?我們只選擇將要使用的變量 head(populardata)?#?我們來(lái)看一下前6個(gè)觀察樣本##?#?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ù)
在開(kāi)始分析之前，我們可以繪制外向性和受歡迎程度之間的關(guān)系，而無(wú)需考慮數(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)。我們可以通過(guò)對(duì)不同類(lèi)進(jìn)行顏色編碼來(lái)顯示這種多層結(jié)構(gòu)。
現(xiàn)在我們可以為數(shù)據(jù)中的100個(gè)不同類(lèi)別繪制不同的回歸線
我們清楚地看到，外向性和受歡迎程度之間的關(guān)系在所有層級(jí)中并不相同，但平均而言，存在明顯的正向關(guān)系。在本教程中，我們將顯示這些不同斜率的估計(jì)值（以及如何解釋這些差異）。?
點(diǎn)擊標(biāo)題查閱往期內(nèi)容
轉(zhuǎn)存失敗重新上傳取消
左右滑動(dòng)查看更多
轉(zhuǎn)存失敗重新上傳取消
01
02
03
04
我們還可以對(duì)最極端的回歸線進(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ù)
截距模型
我們第一個(gè)模型是截距模型。
如果我們查看LMER函數(shù)的不同輸入，則：
“受歡迎程度”，表示我們要預(yù)測(cè)的因變量。
一個(gè)“?”，用于表示我們現(xiàn)在給出了其他感興趣的變量。（與回歸方程式的'='相比）。
公式中表示截距的“ 1”。
由于這是僅截距模式，因此我們?cè)谶@里沒(méi)有任何其他自變量。
在方括號(hào)之間，我們具有隨機(jī)效果/斜率。同樣，值1表示垂直“ |”的截距和變量右側(cè) 條用于指示分組變量。在這種情況下，類(lèi)ID。因此，因變量“受歡迎程度”是由截距和該截距的隨機(jī)誤差項(xiàng)預(yù)測(cè)的。
最后，我們?cè)?code>data =?命令后指定要使用的數(shù)據(jù)集
summary(interceptonlymodel)?#得到參數(shù)估計(jì).##?通過(guò)REML進(jìn)行線性混合模型擬合。t檢驗(yàn)使用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)”下可以看到，類(lèi)別層0.7021上的殘差和第一層（學(xué)生層）上的殘差為1.2218。這意味著類(lèi)內(nèi)相關(guān)性（ICC）為0.7021 /（1.2218 + 0.7021）=.36。
在“固定效果”下，報(bào)告截距的估計(jì)值為5.078。
我們還可以輸出計(jì)算ICC。
##?#?Intraclass?Correlation?Coefficient ##? ##??????Adjusted?ICC:?0.365 ##???Conditional?ICC:?0.365
一層預(yù)測(cè)變量
現(xiàn)在我們可以首先添加第一層（學(xué)生）水平的預(yù)測(cè)變量。一層預(yù)測(cè)因子是性別和外向性?，F(xiàn)在，我們僅將它們添加為固定效果，而不添加為隨機(jī)斜率。在此之前，我們可以繪制兩種性別在效果上的差異。我們發(fā)現(xiàn)性別之間可能存在平均差異，但斜率（回歸系數(shù)）沒(méi)有差異。
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ù)測(cè)變量
現(xiàn)在，我們（除了重要的1層變量）還在第2層（教師經(jīng)驗(yàn)）添加了預(yù)測(cè)變量。
##?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ǔ)模型相比，計(jì)算出第1層和第2層的解釋方差。
對(duì)于1層，這是（1.2218 – 0.592）/1.2218 = 0.52
對(duì)于2層，這是（0.7021 – 0.2954）/0.7021 = 0.58
具有隨機(jī)斜率的一層和二層預(yù)測(cè)模型（1）
現(xiàn)在我們還想包括隨機(jī)斜率。第1層的兩個(gè)預(yù)測(cè)變量（性別和外向性）均具有隨機(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)
我們可以看到所有固定的回歸斜率仍然很顯著。然而，沒(méi)有給出對(duì)隨機(jī)效應(yīng)的顯著性檢驗(yàn)，但是，可變性別的斜率的誤差項(xiàng)（方差）估計(jì)很?。?.0024）。這可能意味著類(lèi)別之間的SEX變量沒(méi)有斜率變化，因此可以從下一次分析中刪除隨機(jī)斜率估計(jì)。由于沒(méi)有針對(duì)此方差的直接顯著性檢驗(yàn)，我們可以使用?軟件包的??ranova()?函數(shù)??lmerTest，提供類(lèi)似于ANOVA的隨機(jī)效果表。它檢查如果刪除了某種隨機(jī)效應(yīng)（稱(chēng)為似然比檢驗(yàn)），則模型是否變得明顯更差，如果不是這種情況，則隨機(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ī)影響確實(shí)不顯著（P = 0.6792），外向的隨機(jī)影響也很顯著（P <.0001）。
具有隨機(jī)斜率的一層和二層預(yù)測(cè)模型
我們?cè)诤雎孕詣e的隨機(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)驗(yàn)的影響是0.091
外向的平均影響為0.453
外向斜率的隨機(jī)效應(yīng)為0.035
一層殘差為0.552
二層的殘差為1.303
具有隨機(jī)斜率和跨水平交互作用的一層和二層預(yù)測(cè)?
作為最后一步，我們可以在教師的經(jīng)驗(yàn)和外向性之間添加跨層的交互作用。換句話說(shuō)，我們要調(diào)查的是，班上外向與受歡迎程度之間關(guān)系的差異是否可以由該班老師的老師經(jīng)驗(yàn)來(lái)解釋。我們添加了Extraversion和Teacher體驗(yàn)之間的層級(jí)交互項(xiàng)。這意味著我們必須添加TEXP作為EXTRAV系數(shù)的預(yù)測(cè)因子。外向性和教師經(jīng)驗(yàn)之間的跨層級(jí)交互作用可以通過(guò)“：”符號(hào)或乘以符號(hào)來(lái)創(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)驗(yàn)j +γ21? extraversionij ?經(jīng)驗(yàn)j + u2j ? extraversionij + u0j + eij
受歡迎程度ij =γ00+γ10? sexij +γ20? extraversionij +γ01?經(jīng)驗(yàn)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
交互項(xiàng)用extrav:texp?表示，??Fixed effects?并估計(jì)為-0.025。
從這些結(jié)果中，我們現(xiàn)在還可以通過(guò)使用教師經(jīng)驗(yàn)作為第二層變量來(lái)計(jì)算解釋的外向斜率方差：（0.03455-0.005409）/0.03455 = .843。因此，外向斜率回歸系數(shù)的方差的84.3％可以由老師的經(jīng)驗(yàn)來(lái)解釋。
外向系數(shù)在受歡迎程度上的截距和斜率均受教師經(jīng)驗(yàn)的影響。一名具有0年經(jīng)驗(yàn)的老師的班級(jí)中，外向得分為0的男學(xué)生（SEX = 0）的預(yù)期受歡迎度為-1.2096。一名類(lèi)似的（男）學(xué)生，每增加1分外向度，就將獲得0.8036分，以提高其知名度。當(dāng)教師經(jīng)驗(yàn)增加時(shí)，每年經(jīng)驗(yàn)的截距也增加0.226。因此，同一個(gè)沒(méi)有外向性的男學(xué)生與一個(gè)有15年經(jīng)驗(yàn)的老師一起上課，其預(yù)期受歡迎度得分為-1.2096 +（15 x .226）= 2.1804。教師的經(jīng)驗(yàn)也減輕了外向性對(duì)普及的影響。對(duì)于具有15年經(jīng)驗(yàn)的教師，外向的回歸系數(shù)僅為0.8036 –（15 x .0247）= 0.4331（相比之下，具有0年經(jīng)驗(yàn)的教師班級(jí)為0.8036）。
我們還可以清楚地看到，多年的教師經(jīng)驗(yàn)既影響截距，又影響外向度的回歸系數(shù)。
最后
在本教程結(jié)束，我們將檢查模型的殘差是否正態(tài)分布（在兩個(gè)層級(jí)上）。除了殘差是正態(tài)分布的之外，多層模型還假設(shè)，對(duì)于不同的隨機(jī)效應(yīng)，殘差的方差在組（類(lèi)）之間是相等的。確實(shí)存在跨組的正態(tài)性和方差相等性的統(tǒng)計(jì)檢驗(yàn)。
首先，我們可以通過(guò)比較殘差和擬合項(xiàng)來(lái)檢查均方差。
我們還可以使用QQ圖檢查殘差的正態(tài)性。該圖確實(shí)表明殘差是正態(tài)分布的。
現(xiàn)在，我們還可以檢查100個(gè)班級(jí)的兩個(gè)隨機(jī)效果。同樣，可以看到符合正態(tài)分布。
點(diǎn)擊文末“閱讀原文”
獲取全文完整資料。
本文選自《R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度》。
本文中的受歡迎數(shù)據(jù)分享到會(huì)員群，掃描下面二維碼即可加群！
點(diǎn)擊標(biāo)題查閱往期內(nèi)容
