拓端tecdat|R語言用線性混合效應(yīng)（多水平/層次/嵌套）模型分析聲調(diào)高低與禮貌態(tài)度的關(guān)
原文鏈接：http://tecdat.cn/?p=23681?
原文出處：拓端數(shù)據(jù)部落公眾號
?
定義
線性混合效應(yīng)模型與我們已經(jīng)知道的線性模型有什么不同？
線性混合模型（有時被稱為 "多層次模型 "或 "層次模型"，取決于上下文）是一種回歸模型，它同時考慮了（1）被感興趣的自變量（如lm()）所解釋的變化--固定效應(yīng)，以及（2）不被感興趣的自變量解釋的變化--隨機效應(yīng)。由于該模型包括固定效應(yīng)和隨機效應(yīng)的混合，所以被稱為混合模型。這些隨機效應(yīng)本質(zhì)上賦予誤差項?結(jié)構(gòu)。
固定效應(yīng)和隨機效應(yīng)的定義可能會有所不同，所以要注意你在文獻中的解釋；但是，對于大多數(shù)目的來說，如果從所有感興趣的層面收集了數(shù)據(jù)，你可以把一個變量視為固定效應(yīng)因素（例如。?性別：男/女，條件：易/中/難，劑量：低/高），如果變量有一堆可能的水平，但你只對一個隨機的集合（如受試者、刺激物、教室）進行采樣，盡管這些樣本會有一些特異性，但你一般不會關(guān)心它們，目的是對更廣泛的人群進行概括（如所有的人、所有的場景、所有的教室）。
例子
比方說，你對語言感興趣，更確切地說，是對聲音的高低與禮貌態(tài)度的關(guān)系感興趣。你要求你的受試者對假設(shè)的場景（IV，受試者內(nèi)部）做出反應(yīng)，這些場景要么是需要禮貌態(tài)度的正式場合（例如，給教授一個遲到的借口），要么是比較非正式的場合（例如，向朋友解釋你為什么遲到），并測量他們的音調(diào)（DV）。每個受試者都會得到一份所有場景的清單，因此每個受試者都會給出多個禮貌態(tài)度的或非正式的回答。你還注意到每個受試者的性別（IV，受試者之間），因為這是對聲調(diào)的另一個重要影響。
在迄今為止我們所看到的線性模型中，我們將建立這樣的模型。
聲調(diào)=禮貌態(tài)度+性別+?
其中最后一項是我們的誤差項。這個誤差項代表了由于我們無法在實驗中控制的 "隨機 "因素而導(dǎo)致的與我們預(yù)測的偏差。
對于這種數(shù)據(jù)，由于每個受試者都給出了多個反應(yīng)（"重復(fù)測量 "設(shè)計），我們可以看到，這將違反線性建模中重要的獨立性假設(shè)：同一受試者的多個反應(yīng)不能被視為彼此獨立。在我們的方案中，每個人的聲調(diào)都略有不同，這將成為影響同一受試者所有反應(yīng)的特異性因素，從而使這些不同的反應(yīng)相互依賴（相關(guān)）而非獨立。
隨機效應(yīng)
我們要處理這種情況的方法是為主體添加一個隨機效應(yīng)。這使我們能夠通過為每個受試者假設(shè)不同的 "基準(zhǔn) "音高值來解決這種非獨立性。因此，受試者1在不同的話語中可能有233赫茲的平均聲調(diào)，而受試者2可能有210赫茲的平均聲調(diào)。在我們的模型中，我們通過對受試者的隨機效應(yīng)來解釋這些聲調(diào)的個體差異。
我們將一些數(shù)據(jù)為例進行分析。
table(subject)
把數(shù)據(jù)可視化。
qplot(condition, pitch, facets = . ~ subject)
受試者 "F#"為女性受試者。對象 "M#"是男性對象。你馬上就會發(fā)現(xiàn)，男性的聲音比女性低（這是可以預(yù)期的）。但除此之外，在男性和女性群體中，你會看到很多個體差異，一些人的性別值相對較高，而另一些人的性別值相對較低。
來自同一主體的樣本的相關(guān)性
另一種說法是，在受試者內(nèi)部，不同條件下的音高存在著相關(guān)性。讓我們把它形象化。
?
用隨機截距對個體平均值進行建模
我們可以通過為每個參與者假設(shè)不同的隨機截距來建立這些個體差異的模型；每個參與者都被分配了不同的截距值（即不同的平均聲調(diào)），而混合模型基本上是為你估計這些截距。
回過頭來看我們的模型，我們以前的公式是。
聲調(diào)=截距+禮貌+性別+?
我們更新后的公式是這樣的。
聲調(diào)=截距+禮貌+性別+(1|個體)+?
"(1|subject) "是隨機截距的R語法。這句話的意思是 "假設(shè)每個主體的截距都不同"......而 "1 "代表這里的截距。你可以認(rèn)為這個公式是告訴你的模型，它應(yīng)該期望每個受試者會有多個反應(yīng)，而這些反應(yīng)將取決于每個受試者的基準(zhǔn)水平。這就有效地解決了因同一受試者有多個反應(yīng)而產(chǎn)生的非獨立性問題。
請注意，該公式仍然包含一個一般誤差項?。這是必要的，因為即使我們考慮到了每個主體的變化，同一主體的不同音高之間仍然會存在 "隨機 "差異。
對不同條件下的不同參與者的平均值有一個概念。
aggregate(pitch ~ subject, FUN = "mean")
現(xiàn)在用lmer() ，我們可以估計每個參與者的平均值。為了做到這一點，我們將為每個受試者包含一個隨機截距，然后看一下估計的截距。?
coef(lmer(pitch ~ (1 | subject))
#固定效應(yīng)+隨機效應(yīng)的主體['（截距）'] + subject
請注意，估計值與實際平均音高相當(dāng)接近，我們可以看到，各受試者的實際平均音高是估計值（Intercept），而各受試者平均音高的標(biāo)準(zhǔn)差是隨機效應(yīng)的標(biāo)準(zhǔn)差（Std.Dev）。
# 使用原始數(shù)據(jù)mean
## [1] 193sd## [1] 63.47# 使用每個子項目的估計截距
mean(subject[1][,'(Intercept)'])
## [1] 193sd## [1] 62.4# 這也是模型輸出中的總結(jié)
summary(res1)
包括固定效應(yīng)
由于我們預(yù)測假設(shè)狀態(tài)的條件（"非正式 "與 "禮貌態(tài)度"）會影響音調(diào)（也許在非正式狀態(tài)下音調(diào)會更高），此外還有受試者的性別（女性的音調(diào)可能會更高），讓我們把這些條件納入模型，同時也考慮到每個受試者的隨機截距（讓截距因受試者而異）。
lmer(音調(diào)~禮貌+性別+(1|個體))
音調(diào)j=截距+截距j+狀態(tài)+性別
音調(diào) A=截距+A截距的變動〔Intercept?Shift〕+狀態(tài)+性別
ggplot(x=condition, y=mean)
在這里，我們還將對比狀態(tài)和性別，這樣我們就可以看到條件在女性和男性之間的 "平均值 "的影響，以及性別在 "非正式 "和 "不禮貌 "之間的平均值的影響。
可以看到，我們的平均音調(diào)是192.88，非正式狀態(tài)的音調(diào)比禮貌狀態(tài)的音調(diào)高，b=9.71，t=3.03，女性的音調(diào)比男性高，b=54.10，t=5.14。我們可以使用一個經(jīng)驗法則，如果t大于2，就可能是顯著的。
更多模型信息
我們可以用許多不同類型的信息來評估模型。在比較模型的時候，這些信息可能很有用 一個有用的衡量標(biāo)準(zhǔn)是AIC，即偏差+2?(p+1)，其中p是模型中的參數(shù)數(shù)量（這里，我們將參數(shù)分解，所以1是估計的殘差，p是所有其他參數(shù)，例如，固定效應(yīng)系數(shù)+估計的隨機效應(yīng)的方差等）。較低的AIC比較好，因為較高的偏差意味著模型不能很好地擬合數(shù)據(jù)。由于AIC隨著p的增加而增加，所以AIC會因為更多的參數(shù)而受到懲罰。
?
deviance = -2*logLikelihood; deviance## [1] 789.5
# 用手計算AICp = 4 # 參數(shù) = 3（固定效應(yīng)）+1（隨機效應(yīng)deviance + 2*(p+1) # 總參數(shù) = 4 + 1 用于估計殘差
## [1] 799.5AIC(res2)## [1] 799.5
提取所有系數(shù)
ggplot(aes(x=factor(subject), y=mean_pitch))
coef(res2)
在這里，我們可以看到，這個模型為每個受試者產(chǎn)生了一個單獨的截距，此外還有一個參數(shù)估計值/斜率的條件和性別，在各受試者中是恒定的。從這里，我們可以嘗試根據(jù)這些系數(shù)來估計一個給定樣本的平均音高 。例如，讓我們嘗試用他們估計的截距和作為女性的影響來估計受試者F1的平均數(shù)（xˉ=232.0357）。
179.3003 + 0*(9.7) + 1*(54.10244)## [1] 233.4
相當(dāng)接近
現(xiàn)在是M3（xˉ=168.9786）:
220.3196 + 0*(9.7) + -1*(54.10244)## [1] 166.2
還不錯
隨機斜率
在上面的模型中，我們假設(shè)禮貌態(tài)度的影響對所有被試都是一樣的，因此，禮貌態(tài)度的系數(shù)是一個。然而，禮貌態(tài)度的影響對于不同的受試者主體可能是不同的；也就是說，可能存在著受試者主體禮貌態(tài)度的相互作用。例如，我們可以預(yù)期，有些人在有需要禮儀的場景下更有禮貌，有些人則不那么有禮貌。因此，我們需要的是一個隨機斜率模型，在這個模型中，不僅允許主體有不同的截距，而且還允許它們對禮貌的影響有不同的斜率（即狀態(tài)對音調(diào)的不同影響）。
?
讓我們開始將數(shù)據(jù)可視化。
根據(jù)這幅圖，看起來各受試者的斜率是否很不穩(wěn)定？
現(xiàn)在加入隨機斜率。
coef(res3)
anova(res2, res3, refit=FALSE)
增加每個受試者的隨機斜率又占用了2個自由度，并沒有明顯改善模型擬合，χ2(2)=0.024,P=0.99。注意df=2，因為我們同時加入了斜率方差和截距與斜率之間的相關(guān)關(guān)系。看一下AIC值，更復(fù)雜的模型的AIC值更高，所以我們想用不太復(fù)雜（更簡明）的模型。如果我們看一下估計的斜率，我們可以看到它們在不同的樣本中是非常相似的。所以，我們不需要在模型中加入條件的隨機斜率。
然而，其他人會認(rèn)為，我們應(yīng)該保持我們的模型最全面，并且應(yīng)該包括隨機斜率，即使它們沒有改善模型!?
測試顯著性
雖然對是否應(yīng)該獲得lmer()模型的p值有一些爭論（例如，這個；大多數(shù)爭論圍繞著如何計算dfs），但你可以使用{lmerTest}包獲得df的近似值（以及因此獲得p值）。
獲取P值
summary(res3b)
將模型輸出與SS/Kenward-Roger appox進行比較
anova
anova(res2b)
模型比較
另一方面，有些人認(rèn)為，用似然比檢驗進行模型比較是檢驗一個參數(shù)是否顯著的更好方法。也就是說，如果在你的模型中加入該參數(shù)能顯著提高模型的擬合度，那么該參數(shù)就應(yīng)該被納入模型中。
似然比檢驗本質(zhì)上告訴我們，數(shù)據(jù)在更復(fù)雜模型下的可能性比在簡單模型下的可能性大多少（這些模型需要嵌套！）：
?
D的分布大約是χ2，自由度為df2-df1。我們要么 "手動 "做這個計算，要么就直接使用anova()函數(shù)!
anova(res4, res4b)
# 手動計算dev0 <- -2*logLik(res4) # 較簡單的偏差模型devdiff <- (dev0-dev1) # 偏差差值dfdiff ?# 參數(shù)的差異（使用dfs）。
在這里，我們比較了兩個嵌套模型，一個沒有條件，另一個有條件。通過模型比較，我們得出結(jié)論，在我們的模型中加入條件是有必要的，因為它明顯改善了模型的擬合，χ2(1)=8.79,P<0.01。
REML與ML
讓我們從一個統(tǒng)計模型開始，指定（i）固定效應(yīng)和（ii）各種隨機效應(yīng)的正態(tài)分布的變異和協(xié)方差。
在ML（最大似然）估計中，我們計算上述（i）和（ii）組中任意選擇的參數(shù)值的數(shù)據(jù)的對數(shù)（似然）（LL）。然后，我們尋找能使L最大化（或最小化-L）的參數(shù)值。這些最佳參數(shù)值被稱為ML參數(shù)估計值。L的最大值（-2倍）被稱為模型的偏差。對于某些目的，如描述數(shù)據(jù)，我們關(guān)注ML參數(shù)估計值；對于其他目的，如模型比較，我們關(guān)注偏差。在比較固定效應(yīng)不同的模型時，你應(yīng)該使用ML，而且你必須包括lmer(, REML=FALSE)。此外，如果你要比較一個lm()和lmer()模型（即測試是否有必要使用任何隨機效應(yīng)），你也應(yīng)該使用ML估計。
在REML（REstricted ML）估計中，我們的主要興趣是估計隨機效應(yīng)，而不是固定效應(yīng)。想象一下，我們通過將上面第(i)組中的固定效應(yīng)參數(shù)設(shè)置為某些合理的值來限制我們的參數(shù)空間。在這個限制的空間里，我們尋找集(ii)中的隨機效應(yīng)參數(shù)值，使數(shù)據(jù)的LL值最大化；同時注意LL值的最大值。然后多次重復(fù)這個過程。然后對固定效應(yīng)參數(shù)值、隨機效應(yīng)參數(shù)的估計值和LL的最大值進行平均。這種平均法可以得到REML參數(shù)估計值和REML偏差值。因為這個過程對固定效應(yīng)參數(shù)的關(guān)注度很低，所以它不應(yīng)該被用來比較固定效應(yīng)結(jié)構(gòu)不同的模型。你應(yīng)該在比較隨機效應(yīng)不同的模型時使用這個方法。要做到這一點，你應(yīng)該養(yǎng)成在運行模型比較時包括lmer(, REML=TRUE)，并使用anova(, refit=FALSE)的習(xí)慣。
例子
測試隨機效應(yīng)是否有必要
dev1 <- -2*logLik(res5)dev0 <- -2*logLik(res5b)devdiff <- as.numeric(dev0-dev1); devdiff
# 比較AICsAIC(res5b, res5)
AICtab(res5b, res5) # ?AIC
??用anova直接比較lm和nlme模型
看起來，是的，納入受試者的隨機截距是有道理的，χ2(1) = 19.51, P < 0.001.
在這里，χ2分布并不總是一個很好的無效分布的近似值（在測試一些隨機效應(yīng)時過于保守，而在測試一些固定效應(yīng)時不夠保守）。
場景效應(yīng)
不同的場景可能會引起不同的 "音調(diào) "值；因此，在不同的受試者之間，甚至在一個受試者內(nèi)部，在禮貌和非正式的狀態(tài)下，特定場景的音調(diào)可能是相關(guān)的。我們可以把這作為一個隨機效應(yīng)的模型。
ggplot(d, aes(x=scenario, y=pitch)
anova
coef(res4)
與受試者的隨機截距相似，現(xiàn)在我們有了每種場景下的平均音高水平。
那么，每個場景的斜率變化呢？
ggplot(byscenario, aes(x=condition, y=pitch))
anova
沒有改善模型。這說明我們的方案在非正式和禮貌狀態(tài)下的差異方面可能是相似的。
混合模型說明
這里有幾件重要的事情要講。你可能會問 "我應(yīng)該指定哪些隨機斜率？" 甚至是 "隨機斜率到底有沒有必要？"
從概念上講，將隨機斜率和隨機截距一起包括進來是非常有意義的。畢竟，你可以認(rèn)為人們對實驗操縱的反應(yīng)不同！同樣，你可以認(rèn)為人們對實驗操縱的反應(yīng)不同。同樣，你總是可以認(rèn)為，實驗操縱的效果對實驗中的所有項目都不一樣。
在上述模型中，我們的整個研究關(guān)鍵在于說明有關(guān)禮貌的東西。我們對性別差異不感興趣，但它們是很值得控制的。這就是為什么我們對禮貌態(tài)度的影響有隨機斜率（按被試和項目），而不是性別。換句話說，在禮貌態(tài)度對音調(diào)的影響方面，我們只模擬了按主體和按項目的變化。
在線性模型背景下討論的一切都直接適用于混合模型。因此，你還必須擔(dān)心共線性和異常值。你還得擔(dān)心同方差（方差相等）和潛在的正態(tài)性缺失問題。
獨立性，作為最重要的假設(shè)，需要一個特殊的詞。我們轉(zhuǎn)向混合模型而不是僅僅使用線性模型的主要原因之一是為了解決我們數(shù)據(jù)中的非依賴性。然而，混合模型仍然可以違反獨立性。如果你缺少重要的固定或隨機效應(yīng)。因此，例如，如果我們用一個不包括隨機效應(yīng) "主體 "的模型來分析我們的數(shù)據(jù)，那么我們的模型就不會 "知道 "每個主體有多個反應(yīng)。這就相當(dāng)于違反了獨立假設(shè)。因此，要仔細(xì)選擇你的固定效應(yīng)和隨機效應(yīng)，解決非獨立性問題。
其他一些說明。
如果你的因變量是...
連續(xù)：使用混合效應(yīng)的線性回歸模型
二元：使用混合效應(yīng)的Logistic回歸模型
函數(shù)lmer用于擬合線性混合模型，函數(shù)glmer用于擬合廣義（非高斯）線性混合模型。
?
最受歡迎的見解
1.基于R語言的lmer混合線性回歸模型
2.R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
3.R語言線性混合效應(yīng)模型實戰(zhàn)案例
4.R語言線性混合效應(yīng)模型實戰(zhàn)案例2
5.R語言線性混合效應(yīng)模型實戰(zhàn)案例
6.線性混合效應(yīng)模型Linear Mixed-Effects Models的部分折疊Gibbs采樣
7.R語言LME4混合效應(yīng)模型研究教師的受歡迎程度
8.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測GDP增長
9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
