【視頻】線性混合效應模型(LMM,Linear Mixed Models)和R語言實現案例

原文鏈接：http://tecdat.cn/?p=23050

原文出處：拓端數據部落公眾號

如果您熟悉線性模型，意識到它們的局限，那么您應該學習線性混合模型mixed-model。本視頻中，我們討論了線性混合模型並在R軟件中進行應用。

視頻：線性混合效應模型(LMM,Linear Mixed Models)和R語言實現

線性混合效應模型(LMM,Linear Mixed Models)和R語言實現案例

 

時長12:13

---

什么是混合效應建模，為什么要使用？

統計分析中許多問題的傳統方法是擬合線性模型，通常使用最小二乘估計。與所有統計方法一樣，最小二乘估計需要做出某些數學假設：數據符合正態分布的並且彼此獨立。

線性統計模型的一個常見示例是多元線性回歸模型：

其中Y被稱為因變量，X是自變量，β是要預測的未知參數，而ϵ是隨機誤差向量。

對於線性回歸模型，我們需要假設誤差是正態分布的並且彼此獨立。自然，嚴重違反這些假設將導致統計模型幾乎沒有用處。

然而，在實際情況中，例如當我們對同一個人重復測量因變量智力分數時，智力分數通常是相關的，因此需要一個模型來解釋這種相關性。

有時因變量顯然不是正態分布的。當我們試圖預測二元因變量時，例如成功/失敗或生存/死亡，誤差只能取兩個值，因此不是正態分布的。但可能通過諸如泊松之類的分布很好地建模。邏輯回歸和泊松回歸分別是在這些情況下使用的模型，並且都是廣義線性模型的特例。

這就是為什么要開發混合模型來處理如此混亂的數據，即使我們的樣本量較小、結構化數據和許多協變量都可以擬合。

線性混合模型

處理相關數據的傳統分析技術是重復測量方差分析和混合模型。相關數據的線性混合模型可以表述為（以回歸模型格式）：

其中  x變量代表固定效應，而  z變量代表隨機效應。

與通常擬合最小二乘的傳統線性模型不同，線性混合模型要么擬合最大似然，要么擬合 REML，限制最大似然。REML 是最大似然的一種變體，通常在變異性估計中具有較小的偏差。

混合模型非常適合聚類數據、重復測量和層次模型。雖然基於經典 ANOVA 的方法可以很好地處理某些特殊情況（例如來自沒有缺失數據的平衡設計的重復測量 ANOVA），但混合模型對於處理更復雜的情況至關重要，包括缺失數據、按不同時間段測量的個體等。

混合模型還可以幫助我們避免假重復的統計錯誤，這是統計推斷中的誤差來源，我們將數據視為獨立的，而實際上並非如此。這導致我們誇大了樣本的大小，從而誇大了自由度和p-值，這可能導致錯誤地得出實際不存在的統計顯着性結論（即 I 類錯誤）。假重復通常發生在具有層次結構的觀察性研究或具有不同空間和/或時間尺度的設計實驗中。

隨機效應和固定效應

噪聲，在統計文獻中被稱為“隨機效應”。指定這些來源決定了我們測量中的相關結構。

在最簡單的線性模型中，我們認為可變性源於測量誤差，因此與其他任何因素無關。但通常是不切實際的。

考慮工業過程控制中的一個問題：測試制造的瓶蓋直徑的變化。我們想研究時間的固定效應：之前與之后。瓶蓋是由幾台機器生產的。很明顯，機器內部和機器之間的直徑存在差異。考慮到來自許多機器的瓶蓋樣本，我們可以通過去除每台機器的平均值來實現測量的標准化。這意味着我們把機器當作固定效應，減去它們，並認為機器內部的變異性是唯一的變異源。減去機器效應后，就去掉了機器間變異性的信息。

另外，在推斷時間固定效應時，我們可以將機器間的變異性視為另一個不確定性的來源。在這種情況下，就不會減去機器效應，而是在LMM框架中把它當作一個隨機效應。

LMM的相關概念

• LMM 涉及到很多基礎概念，因此它有許多名稱：

• 方差分量：因為如示例所示，方差有不止一個來源。

• 分層模型或多級分析：因為我們可以將抽樣視為分層的——首先對類別進行抽樣，然后對其因變量進行抽樣。

• 重復測量：因為我們對每個樣本進行多次測量。

廣義線性混合模型GLMM

廣義線性混合模型相對線性混合模型更加靈活性，即我們可以為因變量假設除正態分布之外的許多族。

廣義線性混合模型的一般形式是

其中  s是固定效應的數量。r是隨機效應的數量。βj是固定效應xij 的參數。bik是隨機效應的參數，而zik是隨機效應的水平。鏈接函數  g(μi)=η用來表示，這樣  y=g(μi) . 因此，混合模型與廣義線性混合模型的結合，形成廣義線性混合模型。

GLMM的鏈接函數

廣義線性混合模型與線性混合模型 之間的不同之處在於因變量可以來自除正態分布之外的不同分布。此外，不是直接對因變量建模，而是應用一些鏈接函數，例如對於二元結果，我們使用Logistic鏈接函數和Logistic的概率密度函數。這些是

對於計數結果，我們使用對數鏈接函數和poisson的概率質量函數，或PMF。請注意，我們稱之為概率質量函數而不是概率密度函數，因為支持是離散的（即對於正整數）。這些是

通過為因變量選擇適當分布族並與線性預測因子相聯系，可以更准確地對具有計數或比例的因變量設計進行建模。隨機效應不再被忽視，而是被估計出來，並且可以對新的數據進行推斷。

R語言對數據進行線性混合效應模型的擬合與可視化

在本文中，我們將用R語言對數據進行線性混合效應模型的擬合，然后可視化你的結果。

線性混合效應模型是在有隨機效應時使用的，隨機效應發生在對隨機抽樣的單位進行多次測量時。來自同一自然組的測量結果本身並不是獨立的隨機樣本。因此，這些單位或群體被假定為從一個群體的 "人口 "中隨機抽取的。示例情況包括

• 當你划分並對各部分進行單獨實驗時（隨機組）。

• 當你的抽樣設計是嵌套的，如橫斷面內的四分儀；林地內的橫斷面；地區內的林地（橫斷面、林地和地區都是隨機組）。

• 當你對相關個體進行測量時（家庭是隨機組）。

• 當你重復測量受試者時（受試者是隨機組）。

混合效應的線性模型在R命令lme4和lmerTest包中實現。另一個選擇是使用nmle包中的lme方法。lme4中用於計算近似自由度的方法比nmle包中的方法更准確一些，特別是在樣本量不大的時候。

測量斑塊長度

這第一個數據集是從Griffith和Sheldon（2001年，《動物行為學》61：987-993）的一篇論文中提取的，他們在兩年內對瑞典哥特蘭島上的30只雄性領頭鶲的白色額斑進行了測量。該斑塊在吸引配偶方面很重要，但其大小每年都有變化。我們在這里的目標是估計斑塊長度（毫米）。

讀取和檢查數據

• 從文件中讀取數據。

• 查看數據的前幾行，看是否正確讀取。

• 創建一個顯示兩年研究中每只飛鳥的測量對圖。可以嘗試制作點陣圖。是否有證據表明不同年份之間存在着測量變異性？

構建線性混合效應模型

• 對數據進行線性混合效應模型，將單個鳥類視為隨機組。注：對每只鳥的兩次測量是在研究的連續年份進行的。為了簡單起見，在模型中不包括年份。在R中把它轉換成一個字符或因子，這樣它就不會被當作一個數字變量。按照下面步驟（2）和（3）所述，用這個模型重新計算可重復性。重復性的解釋如何改變？

• 從保存的lmer對象中提取參數估計值（系數）。檢查隨機效應的輸出。隨機變異的兩個來源是什么？固定效應指的是什么？

• 在輸出中，檢查隨機效應的標准差。應該有兩個標准差：一個是"（截距）"，一個是 "殘差"。這是因為混合效應模型有兩個隨機變異的來源：鳥類內部重復測量的差異，以及鳥類之間額斑長度的真實差異。這兩個來源中的哪一個對應於"（截距）"，哪一個對應於 "殘差"？

• 同時檢查固定效應結果的輸出。模型公式中唯一的固定效應是所有長度測量的平均值。它被稱為"（截距）"，但不要與隨機效應的截距相混淆。固定效應輸出給了你平均值的估計值和該估計值的標准誤差。注意固定效應輸出是如何提供均值估計值的，而隨機效應輸出則提供方差（或標准差）的估計值。

• 從擬合模型中提取方差分量，估計各年斑塊長度的可重復性*。

• 解釋上一步中獲得的重復性測量結果。如果你得到的重復性小於1.0，那么個體內測量結果之間的變化來源是什么。僅是測量誤差嗎？

• 產生一個殘差與擬合值的圖。注意到有什么問題？似乎有一個輕微的正向趨勢。這不是一個錯誤，而是最佳線性無偏預測器（BLUPs）"收縮 "的結果。

分析步驟

讀取並檢查數據。

head(fly)

1.  
   # 點陣圖
2.  
   chart(patch ~ bird)

1.  
   # 但顯示成對數據的更好方法是用成對的交互圖來顯示
2.  
   plot(res=patch, x = year)

1.  
   # 優化版本
2.  
   plot(y = patch, x = factor(year), theme_classic)

擬合一個線性混合效應模型。summary()的輸出將顯示兩個隨機變異的來源：單個鳥類之間的變異（鳥類截距），以及對同一鳥類進行的重復測量之間的變異（殘差）。每個來源都有一個估計的方差和標准差。固定效應只是所有鳥類的平均值--另一個 "截距"。

1.  
   # 1.混合效應模型
2.  
   # 2. 參數估計
3.  
   summary(z)

1.  
   # 5. 方差分量
2.  
   VarCorr(z)

1.  
   # 可重復性
2.  
   1.11504^2/(1.11504^2 + 0.59833^2)

## \[1\] 0.7764342

1.  
   # 7.殘差與擬合值的關系圖
2.  
   plot(z)

---

金魚視覺

Cronly-Dillon和Muntz(1965; J. Exp. Biol 42: 481-493)用視運動反應來測量金魚的色覺。在這里，我們將對數據進行擬合，包括測試的全部波長。5條魚中的每一條都以隨機的順序在所有的波長下被測試。敏感度的值大表明魚可以檢測到低的光強度。視運動反應的一個重要特點是，魚不習慣，在一個波長下的視覺敏感度的測量不太可能對后來在另一個波長下的測量產生影響。

讀取和檢查數據

• 讀取文件中的數據，並查看前幾行以確保讀取正確。

• 使用交互圖來比較不同光波長實驗下的個體魚的反應。

• 使用什么類型的實驗設計？*這將決定在擬合數據時使用的線性混合模型。

構建線性混合效應模型

• 對數據擬合一個線性混合效應模型。可以用lmer()來實現。發現“畸形擬合”，“boundary (singular) fit: see ?isSingular
  ”

• 繪制擬合（預測）值**。每條魚的預測值和觀察值之間的差異代表殘差。

• 你在（1）中做了什么假設？創建一個殘差與擬合值的圖，以檢查這些假設之一。

• 從保存的lmer對象中提取參數估計值。檢查固定效應的結果。給出的系數與使用lm分析的分類變量的解釋相同。

• 檢查隨機效應的輸出。我們的混合效應模型中再次出現了兩個隨機誤差的來源。它們是什么？其中哪個對應於輸出中的"（截距）"，哪個對應於 "殘差"？注意，在這個數據集中，其中一個變化源的估計標准差非常小。這就是畸形擬合信息背后的原因。魚類之間的方差不太可能真的為零，但是這個數據集非常小，由於抽樣誤差，可能會出現低方差估計。

• 生成基於模型的每個波長的平均敏感度的估計。

• 各個波長之間的差異是否顯著？生成lmer對象的方差分析表。這里測試的是什么效應，隨機效應還是固定效應？解釋方差分析結果。

*這是一個 "按實驗對象 "的重復測量設計，因為每條魚在每個實驗下被測量一次。它本質上與隨機完全區塊設計相同（把每條魚看作是 "區塊"）。

*可視化是首選，因為數據和擬合值都被繪制出來。請注意魚與魚之間的預測值是多么的相似。這表明在這項研究中，個體魚之間的估計差異非常小。

*一般來說，在方差分析表中只測試固定效應。使用測試隨機效應中沒有方差的無效假設是可能的。

分析步驟

讀取並檢查數據。

1.  
   x <- read.csv("fish.csv", 
2.  
           stringsAsFactors = FALSE)
3.  
   head(x)

擬合一個線性混合效應模型。

該模型假設所有擬合值的殘差為正態分布，方差相等。該方法還假設個體魚之間的隨機截距為正態分布。該方法還假設組（魚）的隨機抽樣，對同一魚的測量之間沒有影響。

# # 1. 擬合混合效應模型。

## boundary (singular) fit: see ?isSingular

1.  
   # 2. 這就為每條魚分別繪制了擬合值。
2.  
   vis(z)

1.  
   # 3.測試假設
2.  
   plot(z)

1.  
   # 4. 提取參數估計值
2.  
   summary(z)

1.  
   # 6.  基於模型的平均敏感度估計 
2.  
   means(z)

# 7. ANOVA方差分析

---

蓍草酚類物質的濃度

項目實驗性地調查了國家公園的北方森林生態系統中施肥和食草的影響（Krebs, C.J., Boutin, S. & Boonstra, R., eds (2001a) Ecosystem dynamics of the Boreal Forest.Kluane項目. 牛津大學出版社，紐約）) ，目前的數據來自於一項關於植物資源和食草動物對底層植物物種防御性化學的影響的研究。

16個5x5米的小區中的每一個都被隨機分配到四個實驗之一。1）用柵欄圍起來排除食草動物；2）用N-P-K肥料施肥；3）用柵欄和施肥；4）未實驗的對照。然后，16塊地中的每一塊被分成兩塊。每塊地的一側（隨機選擇）在20年的研究中持續接受實驗。每塊地的另一半在頭十年接受實驗，之后讓它恢復到未實驗的狀態。這里要分析的數據記錄了歐蓍草（Achillea millefolium）中酚類物質的濃度（對植物防御化合物的粗略測量），歐蓍草是地塊中常見的草本植物。測量單位是每克干重毫克丹寧酸當量。

可視化數據

• 從文件中讀取數據。

• 檢查前幾行的數據。實驗是作為一個有四個層次的單一變量給出的（而不是作為兩個變量，圍牆和肥料，用2x2因子設計的模型）。持續時間表示半塊土地是否接受了整整20年的實驗，或者是否在10年后停止實驗。變量 "ch "是蓍草中酚類物質的濃度。

• 畫一張圖來說明不同實驗和持續時間類別中蓍草中的酚類物質的濃度。在每個實驗和持續時間水平的組合中沒有很多數據點，所以按組畫條形圖可能比按組畫箱形圖更好。

• 添加線段來連接成對的點。

擬合一個線性混合效應模型

• 使用的是什么類型的實驗設計？*這將決定對數據的線性混合模型的擬合。

• 在沒有實驗和持續時間之間的交互作用的情況下，對數據進行線性混合模型擬合。使用酚類物質的對數作為因變量，因為對數轉換改善了數據與線性模型假設的擬合。

• 可視化模型對數據的擬合。按持續時間（如果xvar是實驗）或實驗（如果xvar是持續時間）分開面板。visreg()不會保留配對，但會允許你檢查殘差。

• 現在重復模型擬合，但這次包括實驗和持續時間之間的相互作用。將模型與數據的擬合情況可視化。兩個模型擬合之間最明顯的區別是什么，一個有交互作用，另一個沒有？描述包括交互項的模型 "允許 "什么，而沒有交互項的模型則不允許。判斷，哪個模型最適合數據？

• 使用診斷圖檢查包括交互項的模型的線性混合模型的一個關鍵假設。

• 使用擬合模型對象估計線性模型的參數（包括交互作用）。請注意，現在固定效應表中有許多系數。

• 在上一步的輸出中，你會看到 "隨機效應 "標簽下的 "Std.Dev "的兩個數量。解釋一下這些數量指的是什么。

• 來估計所有固定效應組合的模型擬合平均值。

• 生成固定效應的方差分析表。哪些項在統計學上是顯著的？

• 默認情況下，lmerTest將使用Type 3的平方和來測試模型項，而不是按順序（Type 1）。用類型1來重復方差分析表。結果有什么不同嗎？**

*實驗采用了分塊設計，即整個塊被隨機分配到不同的實驗，然后將第二種實驗（持續時間）的不同水平分配到塊的一半。

*應該沒有差別，因為設計是完全平衡的。

分析步驟

閱讀並檢查數據。

一個好的策略是對實驗類別進行排序，把對照組放在前面。這將使線性模型的輸出更加有用。

1.  
   # 1. 讀取數據
2.  
   # 2. 檢查
3.  
   head(x)

1.  
   # 3. 分組帶狀圖
2.  
   # 首先，重新排列實驗類別
3.  
   factor(treat,levels=c("cont","exc","fer","bo"))
4.  
   plot(data = x, y = log(phe), x = trea)

1.  
   # 4. 在多個面板上分別繪制成對的數據
2.  
   plot(data = x,y = log(ach, x = dur, fill = dur, col = dur)

擬合一個線性混合效應模型。固定效應是 "實驗 "和 "持續時間"，而 "塊"是隨機效應。擬合交互作用時，實驗水平之間的差異大小在持續時間水平之間會有所不同。

由於隨機效應也存在（塊），系數表將顯示兩個隨機變化來源的方差估計。一個是擬合模型的殘差的方差。第二個是（隨機）塊截距之間的方差。

# 2. 擬合混合效應模型-無交互作用

1.  
   # 3. 可視化
2.  
   vis(z)

1.  
   # 4. 包括交互項和再次視覺化
2.  
    
3.  
   vis(z.int, overlay = TRUE)

1.  
   # 5. 繪制圖表以檢驗方差齊性（以及正態性）
2.  
   plot(z)

1.  
   # 6. 系數
2.  
   summary(z)

1.  
   # 8. 模型擬合平均值
2.  
   means(z, data = x)

1.  
   # 9. 方差分析表
2.  
   anova(z) # lmerTest中默認為3類平方和。

1.  
   # 10.  改為1類
2.  
   anova(z, type = 1)

---

本文摘選《R語言線性混合效應模型（固定效應&隨機效應）和交互可視化3案例》，點擊“閱讀原文”獲取全文完整資料。

---

點擊標題查閱往期內容

非線性混合效應 NLME模型對抗哮喘葯物茶鹼動力學研究

生態學模擬對廣義線性混合模型GLMM進行功率（功效、效能、效力）分析power analysis環境監測數據

有限混合模型聚類FMM、廣義線性回歸模型GLM混合應用分析威士忌市場和研究專利申請數據

如何用潛類別混合效應模型(Latent Class Mixed Model ,LCMM)分析老年痴呆年齡數據

R語言用lme4多層次（混合效應）廣義線性模型（GLM），邏輯回歸分析教育留級調查數據

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

R語言混合效應邏輯回歸（mixed effects logistic）模型分析肺癌數據

R語言如何用潛類別混合效應模型（LCMM）分析抑郁症狀

R語言基於copula的貝葉斯分層混合模型的診斷准確性研究

R語言建立和可視化混合效應模型mixed effect model

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

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

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研究整容手術數據

用SPSS估計HLM多層（層次）線性模型模型