全文下載鏈接：http://tecdat.cn/?p=23681

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

線性混合效應模型與我們已經(jīng)知道的線性模型有什么不同？

線性混合模型（有時被稱為 "多層次模型 "或 "層次模型"，取決于上下文）是一種回歸模型，它同時考慮了（1）被感興趣的自變量（如lm()）所解釋的變化--固定效應，以及（2）不被感興趣的自變量解釋的變化--隨機效應。由于該模型包括固定效應和隨機效應的混合，所以被稱為混合模型。這些隨機效應本質(zhì)上賦予誤差項?結(jié)構。

固定效應和隨機效應的定義可能會有所不同，所以要注意你在文獻中的解釋；但是，對于大多數(shù)目的來說，如果從所有感興趣的層面收集了數(shù)據(jù)，你可以把一個變量視為固定效應因素（例如。?性別：男/女，條件：易/中/難，劑量：低/高），如果變量有一堆可能的水平，但你只對一個隨機的集合（如受試者、刺激物、教室）進行采樣，盡管這些樣本會有一些特異性，但你一般不會關心它們，目的是對更廣泛的人群進行概括（如所有的人、所有的場景、所有的教室）。

例子

比方說，你對語言感興趣，更確切地說，是對聲音的高低與禮貌態(tài)度的關系感興趣。你要求你的受試者對假設的場景（IV，受試者內(nèi)部）做出反應，這些場景要么是需要禮貌態(tài)度的正式場合（例如，給教授一個遲到的借口），要么是比較非正式的場合（例如，向朋友解釋你為什么遲到），并測量他們的音調(diào)（DV）。每個受試者都會得到一份所有場景的清單，因此每個受試者都會給出多個禮貌態(tài)度的或非正式的回答。你還注意到每個受試者的性別（IV，受試者之間），因為這是對聲調(diào)的另一個重要影響。

在迄今為止我們所看到的線性模型中，我們將建立這樣的模型。

聲調(diào)=禮貌態(tài)度+性別+?

其中最后一項是我們的誤差項。這個誤差項代表了由于我們無法在實驗中控制的 "隨機 "因素而導致的與我們預測的偏差。

對于這種數(shù)據(jù)，由于每個受試者都給出了多個反應（"重復測量 "設計），我們可以看到，這將違反線性建模中重要的獨立性假設：同一受試者的多個反應不能被視為彼此獨立。在我們的方案中，每個人的聲調(diào)都略有不同，這將成為影響同一受試者所有反應的特異性因素，從而使這些不同的反應相互依賴（相關）而非獨立。

隨機效應

我們要處理這種情況的方法是為主體添加一個隨機效應。這使我們能夠通過為每個受試者假設不同的 "基準 "音高值來解決這種非獨立性。因此，受試者1在不同的話語中可能有233赫茲的平均聲調(diào)，而受試者2可能有210赫茲的平均聲調(diào)。在我們的模型中，我們通過對受試者的隨機效應來解釋這些聲調(diào)的個體差異。

我們將一些數(shù)據(jù)為例進行分析。

table(subject)

把數(shù)據(jù)可視化。

qplot(condition,?pitch,?facets?=?.?~?subject)

受試者 "F#"為女性受試者。對象 "M#"是男性對象。你馬上就會發(fā)現(xiàn)，男性的聲音比女性低（這是可以預期的）。但除此之外，在男性和女性群體中，你會看到很多個體差異，一些人的性別值相對較高，而另一些人的性別值相對較低。

來自同一主體的樣本的相關性

另一種說法是，在受試者內(nèi)部，不同條件下的音高存在著相關性。讓我們把它形象化。

用隨機截距對個體平均值進行建模

我們可以通過為每個參與者假設不同的隨機截距來建立這些個體差異的模型；每個參與者都被分配了不同的截距值（即不同的平均聲調(diào)），而混合模型基本上是為你估計這些截距。

回過頭來看我們的模型，我們以前的公式是。

聲調(diào)=截距+禮貌+性別+?

我們更新后的公式是這樣的。

聲調(diào)=截距+禮貌+性別+(1|個體)+?

"(1|subject) "是隨機截距的R語法。這句話的意思是 "假設每個主體的截距都不同"......而 "1 "代表這里的截距。你可以認為這個公式是告訴你的模型，它應該期望每個受試者會有多個反應，而這些反應將取決于每個受試者的基準水平。這就有效地解決了因同一受試者有多個反應而產(chǎn)生的非獨立性問題。

請注意，該公式仍然包含一個一般誤差項?。這是必要的，因為即使我們考慮到了每個主體的變化，同一主體的不同音高之間仍然會存在 "隨機 "差異。

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

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

左右滑動查看更多

01

02

03

04

對不同條件下的不同參與者的平均值有一個概念。

?aggregate(pitch?~?subject,?FUN?=?"mean")

現(xiàn)在用lmer() ，我們可以估計每個參與者的平均值。為了做到這一點，我們將為每個受試者包含一個隨機截距，然后看一下估計的截距。

coef(lmer(pitch?~?(1?|?subject))

#固定效應+隨機效應的主體['（截距）']?+?subject

請注意，估計值與實際平均音高相當接近，我們可以看到，各受試者的實際平均音高是估計值（Intercept），而各受試者平均音高的標準差是隨機效應的標準差（Std.Dev）。

#?使用原始數(shù)據(jù)mean##?[1]?193sd##?[1]?63.47#?使用每個子項目的估計截距 mean(subject[1][,'(Intercept)'])##?[1]?193sd##?[1]?62.4#?這也是模型輸出中的總結(jié)summary(res1)

包括固定效應

由于我們預測假設狀態(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，就可能是顯著的。

更多模型信息

我們可以用許多不同類型的信息來評估模型。在比較模型的時候，這些信息可能很有用 一個有用的衡量標準是AIC，即偏差+2?(p+1)，其中p是模型中的參數(shù)數(shù)量（這里，我們將參數(shù)分解，所以1是估計的殘差，p是所有其他參數(shù)，例如，固定效應系數(shù)+估計的隨機效應的方差等）。較低的AIC比較好，因為較高的偏差意味著模型不能很好地擬合數(shù)據(jù)。由于AIC隨著p的增加而增加，所以AIC會因為更多的參數(shù)而受到懲罰。

deviance?=?-2*logLikelihood;?deviance##?[1]?789.5#?用手計算AICp?=?4?#?參數(shù)?=?3（固定效應）+1（隨機效應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

相當接近

現(xiàn)在是M3（xˉ=168.9786）:

220.3196?+?0*(9.7)?+?-1*(54.10244)##?[1]?166.2

還不錯

隨機斜率

在上面的模型中，我們假設禮貌態(tài)度的影響對所有被試都是一樣的，因此，禮貌態(tài)度的系數(shù)是一個。然而，禮貌態(tài)度的影響對于不同的受試者主體可能是不同的；也就是說，可能存在著受試者主體禮貌態(tài)度的相互作用。例如，我們可以預期，有些人在有需要禮儀的場景下更有禮貌，有些人則不那么有禮貌。因此，我們需要的是一個隨機斜率模型，在這個模型中，不僅允許主體有不同的截距，而且還允許它們對禮貌的影響有不同的斜率（即狀態(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，因為我們同時加入了斜率方差和截距與斜率之間的相關關系?？匆幌翧IC值，更復雜的模型的AIC值更高，所以我們想用不太復雜（更簡明）的模型。如果我們看一下估計的斜率，我們可以看到它們在不同的樣本中是非常相似的。所以，我們不需要在模型中加入條件的隨機斜率。

然而，其他人會認為，我們應該保持我們的模型最全面，并且應該包括隨機斜率，即使它們沒有改善模型!

測試顯著性

雖然對是否應該獲得lmer()模型的p值有一些爭論（例如，這個；大多數(shù)爭論圍繞著如何計算dfs），但你可以使用{lmerTest}包獲得df的近似值（以及因此獲得p值）。

獲取P值

summary(res3b)

將模型輸出與SS/Kenward-Roger appox進行比較

anova

anova(res2b)

模型比較

另一方面，有些人認為，用似然比檢驗進行模型比較是檢驗一個參數(shù)是否顯著的更好方法。也就是說，如果在你的模型中加入該參數(shù)能顯著提高模型的擬合度，那么該參數(shù)就應該被納入模型中。

似然比檢驗本質(zhì)上告訴我們，數(shù)據(jù)在更復雜模型下的可能性比在簡單模型下的可能性大多少（這些模型需要嵌套！）：

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）固定效應和（ii）各種隨機效應的正態(tài)分布的變異和協(xié)方差。

• 在ML（最大似然）估計中，我們計算上述（i）和（ii）組中任意選擇的參數(shù)值的數(shù)據(jù)的對數(shù)（似然）（LL）。然后，我們尋找能使L最大化（或最小化-L）的參數(shù)值。這些最佳參數(shù)值被稱為ML參數(shù)估計值。L的最大值（-2倍）被稱為模型的偏差。對于某些目的，如描述數(shù)據(jù)，我們關注ML參數(shù)估計值；對于其他目的，如模型比較，我們關注偏差。在比較固定效應不同的模型時，你應該使用ML，而且你必須包括lmer(, REML=FALSE)。此外，如果你要比較一個lm()和lmer()模型（即測試是否有必要使用任何隨機效應），你也應該使用ML估計。

• 在REML（REstricted ML）估計中，我們的主要興趣是估計隨機效應，而不是固定效應。想象一下，我們通過將上面第(i)組中的固定效應參數(shù)設置為某些合理的值來限制我們的參數(shù)空間。在這個限制的空間里，我們尋找集(ii)中的隨機效應參數(shù)值，使數(shù)據(jù)的LL值最大化；同時注意LL值的最大值。然后多次重復這個過程。然后對固定效應參數(shù)值、隨機效應參數(shù)的估計值和LL的最大值進行平均。這種平均法可以得到REML參數(shù)估計值和REML偏差值。因為這個過程對固定效應參數(shù)的關注度很低，所以它不應該被用來比較固定效應結(jié)構不同的模型。你應該在比較隨機效應不同的模型時使用這個方法。要做到這一點，你應該養(yǎng)成在運行模型比較時包括lmer(, REML=TRUE)，并使用anova(, refit=FALSE)的習慣。

例子

測試隨機效應是否有必要

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分布并不總是一個很好的無效分布的近似值（在測試一些隨機效應時過于保守，而在測試一些固定效應時不夠保守）。

場景效應

不同的場景可能會引起不同的 "音調(diào) "值；因此，在不同的受試者之間，甚至在一個受試者內(nèi)部，在禮貌和非正式的狀態(tài)下，特定場景的音調(diào)可能是相關的。我們可以把這作為一個隨機效應的模型。

ggplot(d,?aes(x=scenario,?y=pitch)

anova

coef(res4)

與受試者的隨機截距相似，現(xiàn)在我們有了每種場景下的平均音高水平。

那么，每個場景的斜率變化呢？

??ggplot(byscenario,?aes(x=condition,?y=pitch))

anova

沒有改善模型。這說明我們的方案在非正式和禮貌狀態(tài)下的差異方面可能是相似的。

混合模型說明

這里有幾件重要的事情要講。你可能會問 "我應該指定哪些隨機斜率？" 甚至是 "隨機斜率到底有沒有必要？"

從概念上講，將隨機斜率和隨機截距一起包括進來是非常有意義的。畢竟，你可以認為人們對實驗操縱的反應不同！同樣，你可以認為人們對實驗操縱的反應不同。同樣，你總是可以認為，實驗操縱的效果對實驗中的所有項目都不一樣。

在上述模型中，我們的整個研究關鍵在于說明有關禮貌的東西。我們對性別差異不感興趣，但它們是很值得控制的。這就是為什么我們對禮貌態(tài)度的影響有隨機斜率（按被試和項目），而不是性別。換句話說，在禮貌態(tài)度對音調(diào)的影響方面，我們只模擬了按主體和按項目的變化。

在線性模型背景下討論的一切都直接適用于混合模型。因此，你還必須擔心共線性和異常值。你還得擔心同方差（方差相等）和潛在的正態(tài)性缺失問題。

獨立性，作為最重要的假設，需要一個特殊的詞。我們轉(zhuǎn)向混合模型而不是僅僅使用線性模型的主要原因之一是為了解決我們數(shù)據(jù)中的非依賴性。然而，混合模型仍然可以違反獨立性。如果你缺少重要的固定或隨機效應。因此，例如，如果我們用一個不包括隨機效應 "主體 "的模型來分析我們的數(shù)據(jù)，那么我們的模型就不會 "知道 "每個主體有多個反應。這就相當于違反了獨立假設。因此，要仔細選擇你的固定效應和隨機效應，解決非獨立性問題。

其他一些說明。

如果你的因變量是...

• 連續(xù)：使用混合效應的線性回歸模型

• 二元：使用混合效應的Logistic回歸模型

函數(shù)lmer用于擬合線性混合模型，函數(shù)glmer用于擬合廣義（非高斯）線性混合模型。

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《R語言用線性混合效應（多水平/層次/嵌套）模型分析聲調(diào)高低與禮貌態(tài)度的關系》。

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

R語言建立和可視化混合效應模型mixed effect model
R語言 線性混合效應模型實戰(zhàn)案例
R語言用潛類別混合效應模型(Latent Class Mixed Model ,LCMM)分析老年癡呆年齡數(shù)據(jù)
R語言貝葉斯廣義線性混合（多層次/水平/嵌套）模型GLMM、邏輯回歸分析教育留級影響因素數(shù)據(jù)R語言估計多元標記的潛過程混合效應模型（lcmm）分析心理測試的認知過程
R語言因子實驗設計nlme擬合非線性混合模型分析有機農(nóng)業(yè)施氮水平
R語言非線性混合效應 NLME模型(固定效應&隨機效應)對抗哮喘藥物茶堿動力學研究
R語言用線性混合效應（多水平/層次/嵌套）模型分析聲調(diào)高低與禮貌態(tài)度的關系
R語言LME4混合效應模型研究教師的受歡迎程度R語言nlme、nlmer、lme4用（非）線性混合模型non-linear mixed model分析藻類數(shù)據(jù)實例
R語言混合線性模型、多層次模型、回歸模型分析學生平均成績GPA和可視化
R語言線性混合效應模型（固定效應&隨機效應）和交互可視化3案例
R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調(diào)查數(shù)據(jù)R語言 線性混合效應模型實戰(zhàn)案例
R語言混合效應邏輯回歸（mixed effects logistic）模型分析肺癌數(shù)據(jù)
R語言如何用潛類別混合效應模型（LCMM）分析抑郁癥狀
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言建立和可視化混合效應模型mixed effect model
R語言LME4混合效應模型研究教師的受歡迎程度
R語言 線性混合效應模型實戰(zhàn)案例
R語言用Rshiny探索lme4廣義線性混合模型（GLMM）和線性混合模型（LMM）
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言如何解決線性混合模型中畸形擬合(Singular fit)的問題
基于R語言的lmer混合線性回歸模型
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
R語言分層線性模型案例
R語言用WinBUGS 軟件對學術能力測驗（SAT）建立分層模型
使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
SPSS中的多層（等級）線性模型Multilevel linear models研究整容手術數(shù)據(jù)
用SPSS估計HLM多層（層次）線性模型模型