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

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

在本文中，我們將用R語(yǔ)言對(duì)數(shù)據(jù)進(jìn)行線性混合效應(yīng)模型的擬合，然后可視化你的結(jié)果。

線性混合效應(yīng)模型是在有隨機(jī)效應(yīng)時(shí)使用的，隨機(jī)效應(yīng)發(fā)生在對(duì)隨機(jī)抽樣的單位進(jìn)行多次測(cè)量時(shí)。來(lái)自同一自然組的測(cè)量結(jié)果本身并不是獨(dú)立的隨機(jī)樣本。因此，這些單位或群體被假定為從一個(gè)群體的 "人口 "中隨機(jī)抽取的。示例情況包括

• 當(dāng)你劃分并對(duì)各部分進(jìn)行單獨(dú)實(shí)驗(yàn)時(shí)（隨機(jī)組）。

• 當(dāng)你的抽樣設(shè)計(jì)是嵌套的，如橫斷面內(nèi)的四分儀；林地內(nèi)的橫斷面；地區(qū)內(nèi)的林地（橫斷面、林地和地區(qū)都是隨機(jī)組）。

• 當(dāng)你對(duì)相關(guān)個(gè)體進(jìn)行測(cè)量時(shí)（家庭是隨機(jī)組）。

• 當(dāng)你重復(fù)測(cè)量受試者時(shí)（受試者是隨機(jī)組）。

混合效應(yīng)的線性模型在R命令lme4和lmerTest包中實(shí)現(xiàn)。另一個(gè)選擇是使用nmle包中的lme方法。lme4中用于計(jì)算近似自由度的方法比nmle包中的方法更準(zhǔn)確一些，特別是在樣本量不大的時(shí)候。

測(cè)量斑塊長(zhǎng)度

這第一個(gè)數(shù)據(jù)集是從Griffith和Sheldon（2001年，《動(dòng)物行為學(xué)》61：987-993）的一篇論文中提取的，他們?cè)趦赡陜?nèi)對(duì)瑞典哥特蘭島上的30只雄性領(lǐng)頭鶲的白色額斑進(jìn)行了測(cè)量。該斑塊在吸引配偶方面很重要，但其大小每年都有變化。我們?cè)谶@里的目標(biāo)是估計(jì)斑塊長(zhǎng)度（毫米）。

讀取和檢查數(shù)據(jù)

1. 從文件中讀取數(shù)據(jù)。

2. 查看數(shù)據(jù)的前幾行，看是否正確讀取。

3. 創(chuàng)建一個(gè)顯示兩年研究中每只飛鳥(niǎo)的測(cè)量對(duì)圖?？梢試L試制作點(diǎn)陣圖。是否有證據(jù)表明不同年份之間存在著測(cè)量變異性？

構(gòu)建線性混合效應(yīng)模型

1. 對(duì)數(shù)據(jù)進(jìn)行線性混合效應(yīng)模型，將單個(gè)鳥(niǎo)類視為隨機(jī)組。注：對(duì)每只鳥(niǎo)的兩次測(cè)量是在研究的連續(xù)年份進(jìn)行的。為了簡(jiǎn)單起見(jiàn)，在模型中不包括年份。在R中把它轉(zhuǎn)換成一個(gè)字符或因子，這樣它就不會(huì)被當(dāng)作一個(gè)數(shù)字變量。按照下面步驟（2）和（3）所述，用這個(gè)模型重新計(jì)算可重復(fù)性。重復(fù)性的解釋如何改變？

2. 從保存的lmer對(duì)象中提取參數(shù)估計(jì)值（系數(shù)）。檢查隨機(jī)效應(yīng)的輸出。隨機(jī)變異的兩個(gè)來(lái)源是什么？固定效應(yīng)指的是什么？

3. 在輸出中，檢查隨機(jī)效應(yīng)的標(biāo)準(zhǔn)差。應(yīng)該有兩個(gè)標(biāo)準(zhǔn)差：一個(gè)是"（截距）"，一個(gè)是 "殘差"。這是因?yàn)榛旌闲?yīng)模型有兩個(gè)隨機(jī)變異的來(lái)源：鳥(niǎo)類內(nèi)部重復(fù)測(cè)量的差異，以及鳥(niǎo)類之間額斑長(zhǎng)度的真實(shí)差異。這兩個(gè)來(lái)源中的哪一個(gè)對(duì)應(yīng)于"（截距）"，哪一個(gè)對(duì)應(yīng)于 "殘差"？

4. 同時(shí)檢查固定效應(yīng)結(jié)果的輸出。模型公式中唯一的固定效應(yīng)是所有長(zhǎng)度測(cè)量的平均值。它被稱為"（截距）"，但不要與隨機(jī)效應(yīng)的截距相混淆。固定效應(yīng)輸出給了你平均值的估計(jì)值和該估計(jì)值的標(biāo)準(zhǔn)誤差。注意固定效應(yīng)輸出是如何提供均值估計(jì)值的，而隨機(jī)效應(yīng)輸出則提供方差（或標(biāo)準(zhǔn)差）的估計(jì)值。

5. 從擬合模型中提取方差分量，估計(jì)各年斑塊長(zhǎng)度的可重復(fù)性*。

6. 解釋上一步中獲得的重復(fù)性測(cè)量結(jié)果。如果你得到的重復(fù)性小于1.0，那么個(gè)體內(nèi)測(cè)量結(jié)果之間的變化來(lái)源是什么。僅是測(cè)量誤差嗎？

7. 產(chǎn)生一個(gè)殘差與擬合值的圖。注意到有什么問(wèn)題？似乎有一個(gè)輕微的正向趨勢(shì)。這不是一個(gè)錯(cuò)誤，而是最佳線性無(wú)偏預(yù)測(cè)器（BLUPs）"收縮 "的結(jié)果。

分析步驟

讀取并檢查數(shù)據(jù)。

1. 


2. head(fly)

1. # 點(diǎn)陣圖

2. chart(patch ~ bird)

1. # 但顯示成對(duì)數(shù)據(jù)的更好方法是用成對(duì)的交互圖來(lái)顯示

2. plot(res=patch, x = year)

1. # 優(yōu)化版本

2. plot(y = patch, x = factor(year), theme_classic)

擬合一個(gè)線性混合效應(yīng)模型。summary()的輸出將顯示兩個(gè)隨機(jī)變異的來(lái)源：?jiǎn)蝹€(gè)鳥(niǎo)類之間的變異（鳥(niǎo)類截距），以及對(duì)同一鳥(niǎo)類進(jìn)行的重復(fù)測(cè)量之間的變異（殘差）。每個(gè)來(lái)源都有一個(gè)估計(jì)的方差和標(biāo)準(zhǔn)差。固定效應(yīng)只是所有鳥(niǎo)類的平均值--另一個(gè) "截距"。?

1. # 1.混合效應(yīng)模型

2. # 2. 參數(shù)估計(jì)

3. summary(z)

1. # 5. 方差分量

2. VarCorr(z)

1. 


2. # 可重復(fù)性

3. 1.11504^2/(1.11504^2 + 0.59833^2)

## [1] 0.7764342

1. 
   

2. # 7.殘差與擬合值的關(guān)系圖

3. plot(z)

金魚(yú)視覺(jué)

Cronly-Dillon和Muntz(1965; J. Exp. Biol 42: 481-493)用視運(yùn)動(dòng)反應(yīng)來(lái)測(cè)量金魚(yú)的色覺(jué)。在這里，我們將對(duì)數(shù)據(jù)進(jìn)行擬合，包括測(cè)試的全部波長(zhǎng)。5條魚(yú)中的每一條都以隨機(jī)的順序在所有的波長(zhǎng)下被測(cè)試。敏感度的值大表明魚(yú)可以檢測(cè)到低的光強(qiáng)度。視運(yùn)動(dòng)反應(yīng)的一個(gè)重要特點(diǎn)是，魚(yú)不習(xí)慣，在一個(gè)波長(zhǎng)下的視覺(jué)敏感度的測(cè)量不太可能對(duì)后來(lái)在另一個(gè)波長(zhǎng)下的測(cè)量產(chǎn)生影響。

讀取和檢查數(shù)據(jù)

1. 讀取文件中的數(shù)據(jù)，并查看前幾行以確保讀取正確。

2. 使用交互圖來(lái)比較不同光波長(zhǎng)實(shí)驗(yàn)下的個(gè)體魚(yú)的反應(yīng)。

3. 使用什么類型的實(shí)驗(yàn)設(shè)計(jì)？*這將決定在擬合數(shù)據(jù)時(shí)使用的線性混合模型。

構(gòu)建線性混合效應(yīng)模型

1. 對(duì)數(shù)據(jù)擬合一個(gè)線性混合效應(yīng)模型。可以用lmer()來(lái)實(shí)現(xiàn)。發(fā)現(xiàn)“畸形擬合”，“boundary (singular) fit: see ?isSingular
   ”

2. 繪制擬合（預(yù)測(cè)）值**。每條魚(yú)的預(yù)測(cè)值和觀察值之間的差異代表殘差。

3. 你在（1）中做了什么假設(shè)？創(chuàng)建一個(gè)殘差與擬合值的圖，以檢查這些假設(shè)之一。

4. 從保存的lmer對(duì)象中提取參數(shù)估計(jì)值。檢查固定效應(yīng)的結(jié)果。給出的系數(shù)與使用lm分析的分類變量的解釋相同。

5. 檢查隨機(jī)效應(yīng)的輸出。我們的混合效應(yīng)模型中再次出現(xiàn)了兩個(gè)隨機(jī)誤差的來(lái)源。它們是什么？其中哪個(gè)對(duì)應(yīng)于輸出中的"（截距）"，哪個(gè)對(duì)應(yīng)于 "殘差"？注意，在這個(gè)數(shù)據(jù)集中，其中一個(gè)變化源的估計(jì)標(biāo)準(zhǔn)差非常小。這就是畸形擬合信息背后的原因。魚(yú)類之間的方差不太可能真的為零，但是這個(gè)數(shù)據(jù)集非常小，由于抽樣誤差，可能會(huì)出現(xiàn)低方差估計(jì)。

6. 生成基于模型的每個(gè)波長(zhǎng)的平均敏感度的估計(jì)。

7. 各個(gè)波長(zhǎng)之間的差異是否顯著？生成lmer對(duì)象的方差分析表。這里測(cè)試的是什么效應(yīng)，隨機(jī)效應(yīng)還是固定效應(yīng)？解釋方差分析結(jié)果。

*這是一個(gè) "按實(shí)驗(yàn)對(duì)象 "的重復(fù)測(cè)量設(shè)計(jì)，因?yàn)槊織l魚(yú)在每個(gè)實(shí)驗(yàn)下被測(cè)量一次。它本質(zhì)上與隨機(jī)完全區(qū)塊設(shè)計(jì)相同（把每條魚(yú)看作是 "區(qū)塊"）。

*可視化是首選，因?yàn)閿?shù)據(jù)和擬合值都被繪制出來(lái)。請(qǐng)注意魚(yú)與魚(yú)之間的預(yù)測(cè)值是多么的相似。這表明在這項(xiàng)研究中，個(gè)體魚(yú)之間的估計(jì)差異非常小。

***一般來(lái)說(shuō)，在方差分析表中只測(cè)試固定效應(yīng)。使用測(cè)試隨機(jī)效應(yīng)中沒(méi)有方差的無(wú)效假設(shè)是可能的。

分析步驟

讀取并檢查數(shù)據(jù)。

1. x <- read.csv("fish.csv",

2. stringsAsFactors = FALSE)

3. head(x)

擬合一個(gè)線性混合效應(yīng)模型。

該模型假設(shè)所有擬合值的殘差為正態(tài)分布，方差相等。該方法還假設(shè)個(gè)體魚(yú)之間的隨機(jī)截距為正態(tài)分布。該方法還假設(shè)組（魚(yú)）的隨機(jī)抽樣，對(duì)同一魚(yú)的測(cè)量之間沒(méi)有影響。

1. # # 1. 擬合混合效應(yīng)模型。

2. 


## boundary (singular) fit: see ?isSingular

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

2. vis(z)

?

1. # 3.測(cè)試假設(shè)

2. plot(z)

1. # 4. 提取參數(shù)估計(jì)值

2. summary(z)

1. # 6. ?基于模型的平均敏感度估計(jì)

2. means(z)

# 7. ANOVA方差分析

蓍草酚類物質(zhì)的濃度

項(xiàng)目實(shí)驗(yàn)性地調(diào)查了國(guó)家公園的北方森林生態(tài)系統(tǒng)中施肥和食草的影響（Krebs, C.J., Boutin, S. & Boonstra, R., eds (2001a) Ecosystem dynamics of the Boreal Forest.Kluane項(xiàng)目.?牛津大學(xué)出版社，紐約）) ，目前的數(shù)據(jù)來(lái)自于一項(xiàng)關(guān)于植物資源和食草動(dòng)物對(duì)底層植物物種防御性化學(xué)的影響的研究。

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

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

1. 從文件中讀取數(shù)據(jù)。

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

3. 畫(huà)一張圖來(lái)說(shuō)明不同實(shí)驗(yàn)和持續(xù)時(shí)間類別中蓍草中的酚類物質(zhì)的濃度。在每個(gè)實(shí)驗(yàn)和持續(xù)時(shí)間水平的組合中沒(méi)有很多數(shù)據(jù)點(diǎn)，所以按組畫(huà)條形圖可能比按組畫(huà)箱形圖更好。

4. 添加線段來(lái)連接成對(duì)的點(diǎn)。

擬合一個(gè)線性混合效應(yīng)模型

1. 使用的是什么類型的實(shí)驗(yàn)設(shè)計(jì)？*這將決定對(duì)數(shù)據(jù)的線性混合模型的擬合。

2. 在沒(méi)有實(shí)驗(yàn)和持續(xù)時(shí)間之間的交互作用的情況下，對(duì)數(shù)據(jù)進(jìn)行線性混合模型擬合。使用酚類物質(zhì)的對(duì)數(shù)作為因變量，因?yàn)閷?duì)數(shù)轉(zhuǎn)換改善了數(shù)據(jù)與線性模型假設(shè)的擬合。

3. 可視化模型對(duì)數(shù)據(jù)的擬合。按持續(xù)時(shí)間（如果xvar是實(shí)驗(yàn)）或?qū)嶒?yàn)（如果xvar是持續(xù)時(shí)間）分開(kāi)面板。visreg()不會(huì)保留配對(duì)，但會(huì)允許你檢查殘差。

4. 現(xiàn)在重復(fù)模型擬合，但這次包括實(shí)驗(yàn)和持續(xù)時(shí)間之間的相互作用。將模型與數(shù)據(jù)的擬合情況可視化。兩個(gè)模型擬合之間最明顯的區(qū)別是什么，一個(gè)有交互作用，另一個(gè)沒(méi)有？描述包括交互項(xiàng)的模型 "允許 "什么，而沒(méi)有交互項(xiàng)的模型則不允許。判斷，哪個(gè)模型最適合數(shù)據(jù)？

5. 使用診斷圖檢查包括交互項(xiàng)的模型的線性混合模型的一個(gè)關(guān)鍵假設(shè)。

6. 使用擬合模型對(duì)象估計(jì)線性模型的參數(shù)（包括交互作用）。請(qǐng)注意，現(xiàn)在固定效應(yīng)表中有許多系數(shù)。

7. 在上一步的輸出中，你會(huì)看到 "隨機(jī)效應(yīng) "標(biāo)簽下的 "Std.Dev "的兩個(gè)數(shù)量。解釋一下這些數(shù)量指的是什么。

8. 來(lái)估計(jì)所有固定效應(yīng)組合的模型擬合平均值。

9. 生成固定效應(yīng)的方差分析表。哪些項(xiàng)在統(tǒng)計(jì)學(xué)上是顯著的？

10. 默認(rèn)情況下，lmerTest將使用Type 3的平方和來(lái)測(cè)試模型項(xiàng)，而不是按順序（Type 1）。用類型1來(lái)重復(fù)方差分析表。結(jié)果有什么不同嗎？**

*實(shí)驗(yàn)采用了分塊設(shè)計(jì)，即整個(gè)塊被隨機(jī)分配到不同的實(shí)驗(yàn)，然后將第二種實(shí)驗(yàn)（持續(xù)時(shí)間）的不同水平分配到塊的一半。

*應(yīng)該沒(méi)有差別，因?yàn)樵O(shè)計(jì)是完全平衡的。

分析步驟

閱讀并檢查數(shù)據(jù)。

一個(gè)好的策略是對(duì)實(shí)驗(yàn)類別進(jìn)行排序，把對(duì)照組放在前面。這將使線性模型的輸出更加有用。

1. # 1. 讀取數(shù)據(jù)

2. # 2. 檢查

3. head(x)

1. # 3. 分組帶狀圖

2. # 首先，重新排列實(shí)驗(yàn)類別

3. factor(treat,levels=c("cont","exc","fer","bo"))

4. plot(data = x, y = log(phe), x = treat, fill = dura, color = dura)

1. # 4. 在多個(gè)面板上分別繪制成對(duì)的數(shù)據(jù)

2. plot(data = x,y = log(ach, x = dur, fill = dur, col = dur)

擬合一個(gè)線性混合效應(yīng)模型。固定效應(yīng)是 "實(shí)驗(yàn) "和 "持續(xù)時(shí)間"，而 "塊"是隨機(jī)效應(yīng)。擬合交互作用時(shí)，實(shí)驗(yàn)水平之間的差異大小在持續(xù)時(shí)間水平之間會(huì)有所不同。

由于隨機(jī)效應(yīng)也存在（塊），系數(shù)表將顯示兩個(gè)隨機(jī)變化來(lái)源的方差估計(jì)。一個(gè)是擬合模型的殘差的方差。第二個(gè)是（隨機(jī)）塊截距之間的方差。

# 2. 擬合混合效應(yīng)模型-無(wú)交互作用

1. # 3. 可視化

2. vis(z)

1. # 4. 包括交互項(xiàng)和再次視覺(jué)化

2. z.int <- lmer(log(phen.ach) ~ treatment * duration + (1|plot), data=x)

3. 


4. vis(z.int, overlay = TRUE)

1. # 5. 繪制圖表以檢驗(yàn)方差齊性（以及正態(tài)性）

2. plot(z)

1. # 6. 系數(shù)

2. summary(z)

1. # 8. 模型擬合平均值

2. means(z, data = x)

1. # 9. 方差分析表

2. anova(z) # lmerTest中默認(rèn)為3類平方和。

1. # 10. ?改為1類

2. anova(z, type = 1)

最受歡迎的見(jiàn)解

1.基于R語(yǔ)言的lmer混合線性回歸模型

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

3.R語(yǔ)言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

4.R語(yǔ)言線性混合效應(yīng)模型實(shí)戰(zhàn)案例2

5.R語(yǔ)言線性混合效應(yīng)模型實(shí)戰(zhàn)案例

6.線性混合效應(yīng)模型Linear Mixed-Effects Models的部分折疊Gibbs采樣

7.R語(yǔ)言LME4混合效應(yīng)模型研究教師的受歡迎程度

8.R語(yǔ)言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測(cè)GDP增長(zhǎng)

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