拓端tecdat|R語言線性混合效應模型（固定效應&隨機效應）和交互可視化3案例

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

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

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

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

• 當你划分並對各部分進行單獨實驗時（隨機組）。
• 當你的抽樣設計是嵌套的，如橫斷面內的四分儀；林地內的橫斷面；地區內的林地（橫斷面、林地和地區都是隨機組）。
• 當你對相關個體進行測量時（家庭是隨機組）。
• 當你重復測量受試者時（受試者是隨機組）。

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

---

測量斑塊長度

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

讀取和檢查數據

1. 從文件中讀取數據。
2. 查看數據的前幾行，看是否正確讀取。
3. 創建一個顯示兩年研究中每只飛鳥的測量對圖。可以嘗試制作點陣圖。是否有證據表明不同年份之間存在着測量變異性？

構建線性混合效應模型

1. 對數據進行線性混合效應模型，將單個鳥類視為隨機組。注：對每只鳥的兩次測量是在研究的連續年份進行的。為了簡單起見，在模型中不包括年份。在R中把它轉換成一個字符或因子，這樣它就不會被當作一個數字變量。按照下面步驟（2）和（3）所述，用這個模型重新計算可重復性。重復性的解釋如何改變？
2. 從保存的lmer對象中提取參數估計值（系數）。檢查隨機效應的輸出。隨機變異的兩個來源是什么？固定效應指的是什么？
3. 在輸出中，檢查隨機效應的標准差。應該有兩個標准差：一個是"（截距）"，一個是 "殘差"。這是因為混合效應模型有兩個隨機變異的來源：鳥類內部重復測量的差異，以及鳥類之間額斑長度的真實差異。這兩個來源中的哪一個對應於"（截距）"，哪一個對應於 "殘差"？
4. 同時檢查固定效應結果的輸出。模型公式中唯一的固定效應是所有長度測量的平均值。它被稱為"（截距）"，但不要與隨機效應的截距相混淆。固定效應輸出給了你平均值的估計值和該估計值的標准誤差。注意固定效應輸出是如何提供均值估計值的，而隨機效應輸出則提供方差（或標准差）的估計值。
5. 從擬合模型中提取方差分量，估計各年斑塊長度的可重復性*。
6. 解釋上一步中獲得的重復性測量結果。如果你得到的重復性小於1.0，那么個體內測量結果之間的變化來源是什么。僅是測量誤差嗎？
7. 產生一個殘差與擬合值的圖。注意到有什么問題？似乎有一個輕微的正向趨勢。這不是一個錯誤，而是最佳線性無偏預測器（BLUPs）"收縮 "的結果。

分析步驟

讀取並檢查數據。

1.  
    
2.  
   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.  
   # 可重復性
3.  
   1.11504^2/(1.11504^2 + 0.59833^2)

## [1] 0.7764342

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

---

金魚視覺

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

讀取和檢查數據

1. 讀取文件中的數據，並查看前幾行以確保讀取正確。
2. 使用交互圖來比較不同光波長實驗下的個體魚的反應。
3. 使用什么類型的實驗設計？*這將決定在擬合數據時使用的線性混合模型。

構建線性混合效應模型

1. 對數據擬合一個線性混合效應模型。可以用lmer()來實現。發現“畸形擬合”，“boundary (singular) fit: see ?isSingular
   ”
2. 繪制擬合（預測）值**。每條魚的預測值和觀察值之間的差異代表殘差。
3. 你在（1）中做了什么假設？創建一個殘差與擬合值的圖，以檢查這些假設之一。
4. 從保存的lmer對象中提取參數估計值。檢查固定效應的結果。給出的系數與使用lm分析的分類變量的解釋相同。
5. 檢查隨機效應的輸出。我們的混合效應模型中再次出現了兩個隨機誤差的來源。它們是什么？其中哪個對應於輸出中的"（截距）"，哪個對應於 "殘差"？注意，在這個數據集中，其中一個變化源的估計標准差非常小。這就是畸形擬合信息背后的原因。魚類之間的方差不太可能真的為零，但是這個數據集非常小，由於抽樣誤差，可能會出現低方差估計。
6. 生成基於模型的每個波長的平均敏感度的估計。
7. 各個波長之間的差異是否顯著？生成lmer對象的方差分析表。這里測試的是什么效應，隨機效應還是固定效應？解釋方差分析結果。

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

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

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

分析步驟

讀取並檢查數據。

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

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

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

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

## 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）中酚類物質的濃度（對植物防御化合物的粗略測量），歐蓍草是地塊中常見的草本植物。測量單位是每克干重毫克丹寧酸當量。

可視化數據

1. 從文件中讀取數據。
2. 檢查前幾行的數據。實驗是作為一個有四個層次的單一變量給出的（而不是作為兩個變量，圍牆和肥料，用2x2因子設計的模型）。持續時間表示半塊土地是否接受了整整20年的實驗，或者是否在10年后停止實驗。變量 "ch "是蓍草中酚類物質的濃度。
3. 畫一張圖來說明不同實驗和持續時間類別中蓍草中的酚類物質的濃度。在每個實驗和持續時間水平的組合中沒有很多數據點，所以按組畫條形圖可能比按組畫箱形圖更好。
4. 添加線段來連接成對的點。

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

1. 使用的是什么類型的實驗設計？*這將決定對數據的線性混合模型的擬合。
2. 在沒有實驗和持續時間之間的交互作用的情況下，對數據進行線性混合模型擬合。使用酚類物質的對數作為因變量，因為對數轉換改善了數據與線性模型假設的擬合。
3. 可視化模型對數據的擬合。按持續時間（如果xvar是實驗）或實驗（如果xvar是持續時間）分開面板。visreg()不會保留配對，但會允許你檢查殘差。
4. 現在重復模型擬合，但這次包括實驗和持續時間之間的相互作用。將模型與數據的擬合情況可視化。兩個模型擬合之間最明顯的區別是什么，一個有交互作用，另一個沒有？描述包括交互項的模型 "允許 "什么，而沒有交互項的模型則不允許。判斷，哪個模型最適合數據？
5. 使用診斷圖檢查包括交互項的模型的線性混合模型的一個關鍵假設。
6. 使用擬合模型對象估計線性模型的參數（包括交互作用）。請注意，現在固定效應表中有許多系數。
7. 在上一步的輸出中，你會看到 "隨機效應 "標簽下的 "Std.Dev "的兩個數量。解釋一下這些數量指的是什么。
8. 來估計所有固定效應組合的模型擬合平均值。
9. 生成固定效應的方差分析表。哪些項在統計學上是顯著的？
10. 默認情況下，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 = treat, fill = dura, color = dura)

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

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

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

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

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

1.  
   # 4. 包括交互項和再次視覺化
2.  
   z.int <- lmer(log(phen.ach) ~ treatment * duration + (1|plot), data=x)
3.  
    
4.  
   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)

---

最受歡迎的見解

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