原文連結：http://tecdat.cn/?p=23681

原文出處：拓端資料部落公眾號

定義

線性混合效應模型與我們已經知道的線性模型有什麼不同？

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

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

例子

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

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

聲調=禮貌態度+性別+ϵ

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

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

隨機效應

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

我們將一些資料為例進行分析。

table(subject)

把資料視覺化。

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

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

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

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

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

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

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

聲調=截距+禮貌+性別+ϵ

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

聲調=截距+禮貌+性別+(1|個體)+ϵ

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

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

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

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

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

2. coef(lmer(pitch ~ (1 | subject))

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

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

1. # 使用原始資料
2. mean

## [1] 193

sd

## [1] 63.47

1. # 使用每個子專案的估計截距
2. mean(subject[1][,'(Intercept)'])

## [1] 193

sd

## [1] 62.4

1. # 這也是模型輸出中的總結
2. summary(res1)

包括固定效應

由於我們預測假設狀態的條件（"非正式 "與 "禮貌態度"）會影響音調（也許在非正式狀態下音調會更高），此外還有受試者的性別（女性的音調可能會更高），讓我們把這些條件納入模型，同時也考慮到每個受試者的隨機截距（讓截距因受試者而異）。

lmer(音調~禮貌+性別+(1|個體))

音調j=截距+截距j+狀態+性別

音調 A=截距+A截距的變動〔InterceptShift〕+狀態+性別

2. ggplot(x=condition, y=mean)

在這裡，我們還將對比狀態和性別，這樣我們就可以看到條件在女性和男性之間的 "平均值 "的影響，以及性別在 "非正式 "和 "不禮貌 "之間的平均值的影響。

可以看到，我們的平均音調是192.88，非正式狀態的音調比禮貌狀態的音調高，b=9.71，t=3.03，女性的音調比男性高，b=54.10，t=5.14。我們可以使用一個經驗法則，如果t大於2，就可能是顯著的。

更多模型資訊

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

deviance = -2*logLikelihood; deviance

## [1] 789.5

1. # 用手計算AIC
2. p = 4 # 引數 = 3（固定效應）+1（隨機效應
3. deviance + 2*(p+1) # 總引數 = 4 + 1 用於估計殘差

## [1] 799.5

AIC(res2)

## [1] 799.5

提取所有係數

2. ggplot(aes(x=factor(subject), y=mean_pitch))

coef(res2)

在這裡，我們可以看到，這個模型為每個受試者產生了一個單獨的截距，此外還有一個引數估計值/斜率的條件和性別，在各受試者中是恆定的。從這裡，我們可以嘗試根據這些係數來估計一個給定樣本的平均音高 。例如，讓我們嘗試用他們估計的截距和作為女性的影響來估計受試者F1的平均數（x¯=232.0357）。

179.3003 + 0*(9.7) + 1*(54.10244)

## [1] 233.4

相當接近

現在是M3（x¯=168.9786）:

220.3196 + 0*(9.7) + -1*(54.10244)

## [1] 166.2

還不錯

隨機斜率

在上面的模型中，我們假設禮貌態度的影響對所有被試都是一樣的，因此，禮貌態度的係數是一個。然而，禮貌態度的影響對於不同的受試者主體可能是不同的；也就是說，可能存在著受試者主體禮貌態度的相互作用。例如，我們可以預期，有些人在有需要禮儀的場景下更有禮貌，有些人則不那麼有禮貌。因此，我們需要的是一個隨機斜率模型，在這個模型中，不僅允許主體有不同的截距，而且還允許它們對禮貌的影響有不同的斜率（即狀態對音調的不同影響）。

讓我們開始將資料視覺化。

根據這幅圖，看起來各受試者的斜率是否很不穩定？

現在加入隨機斜率。

coef(res3)

anova(res2, res3, refit=FALSE)

增加每個受試者的隨機斜率又佔用了2個自由度，並沒有明顯改善模型擬合，χ2(2)=0.024,P=0.99。注意df=2，因為我們同時加入了斜率方差和截距與斜率之間的相關關係。看一下AIC值，更復雜的模型的AIC值更高，所以我們想用不太複雜（更簡明）的模型。如果我們看一下估計的斜率，我們可以看到它們在不同的樣本中是非常相似的。所以，我們不需要在模型中加入條件的隨機斜率。

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

測試顯著性

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

獲取P值

2. summary(res3b)

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

anova

anova(res2b)

模型比較

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

似然比檢驗本質上告訴我們，資料在更復雜模型下的可能性比在簡單模型下的可能性大多少（這些模型需要巢狀！）：

D的分佈大約是χ2，自由度為df2-df1。我們要麼 "手動 "做這個計算，要麼就直接使用anova()函式!

anova(res4, res4b)

2. # 手動計算
3. dev0 <- -2*logLik(res4) # 較簡單的偏差模型
4. devdiff <- (dev0-dev1) # 偏差差值
5. dfdiff # 引數的差異（使用dfs）。

在這裡，我們比較了兩個巢狀模型，一個沒有條件，另一個有條件。通過模型比較，我們得出結論，在我們的模型中加入條件是有必要的，因為它明顯改善了模型的擬合，χ2(1)=8.79,P<0.01。

REML與ML

讓我們從一個統計模型開始，指定（i）固定效應和（ii）各種隨機效應的正態分佈的變異和協方差。

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

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

例子

測試隨機效應是否有必要

1. dev1 <- -2*logLik(res5)
2. dev0 <- -2*logLik(res5b)
3. devdiff <- as.numeric(dev0-dev1); devdiff

1. # 比較AICs
2. AIC(res5b, res5)

AICtab(res5b, res5) #  AIC

用anova直接比較lm和nlme模型

看起來，是的，納入受試者的隨機截距是有道理的，χ2(1) = 19.51, P < 0.001.

在這裡，χ2分佈並不總是一個很好的無效分佈的近似值（在測試一些隨機效應時過於保守，而在測試一些固定效應時不夠保守）。

場景效應

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

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

anova

coef(res4)

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

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

3. ggplot(byscenario, aes(x=condition, y=pitch))

anova

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

混合模型說明

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

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

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

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

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

其他一些說明。

如果你的因變數是...

• 連續：使用混合效應的線性迴歸模型
• 二元：使用混合效應的Logistic迴歸模型

函式lmer用於擬合線性混合模型，函式glmer用於擬合廣義（非高斯）線性混合模型。

最受歡迎的見解

1.基於R語言的lmer混合線性迴歸模型

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

3.R語言線性混合效應模型實戰案例

4.R語言線性混合效應模型實戰案例2

5.R語言線性混合效應模型實戰案例

6.線性混合效應模型Linear Mixed-Effects Models的部分摺疊Gibbs取樣

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

8.R語言中基於混合資料抽樣(MIDAS)迴歸的HAR-RV模型預測GDP增長

9.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM