★ 鄧飛注：混合線性模型，我以前寫過幾篇學習筆記（線性混合模型系列 1~5），但是偏向基礎知識，看過欒老師的博客后，感覺思路清晰了很多。經(jīng)過欒老師許可，放到公眾號上，我對原來的代碼進行了一些小的更改（有些包的函數(shù)變化了）。非常認可欒老師的觀點：“學習混合線性模型，可以對比一般線性模型，一邊學習理論，一邊動手實踐”。點擊閱讀原文，查看欒老師的博客。另外，欒老師用ggplot2對數(shù)據(jù)和模型進行的可視化分析，非常贊。
” 原文包括兩部分，這里我合二為一。
孿生老師簡介：
以下為正文
學習線性混合效應模型（Linear Mixed Effects Model，LMM）最好的方法，是一邊學習理論，一邊動手實踐，這樣印象最為深刻。本文參考了Bodo Winter博士的教程Linear models and linear mixed effects models in R教程1教程2的結(jié)構(gòu)。
★ 教程1 pdf鏈接：http://www./tutorial/bw_LME_tutorial1.pdf
教程2 pdf鏈接：http://www./tutorial/bw_LME_tutorial2.pdf
” 本文中，為了便于理解，使用的數(shù)據(jù)集是來自csv文件shrimp中的對蝦育種數(shù)據(jù)（對原始數(shù)據(jù)已經(jīng)進行了變換）。
★ 鄧飛注：原始數(shù)據(jù)下載鏈接，https://luansheng./post/datasets/shrimp.csv
” 推薦使用Rstudio來運行R，依賴的R包有：
★ 鄧飛注：這些包，data.table是讀取寫入數(shù)據(jù)，ggplot2作圖包，lme4混合線性包，sjPlot模型作圖包，emmeans計算預測均值，lmerTest是固定因子和隨機因子顯著性檢驗。安裝包：install.packages("data.table") 命令按照 data.table 包，以此類推，這些包都是在CRAN中的，直接使用install.packages進行安裝。
” 1 讀取數(shù)據(jù)文件 fread()函數(shù)來自data.table包，特點是對大文件讀取速度特別快。建議使用data.table作為數(shù)據(jù)處理的主力包。
其中的參數(shù)sep表示文件是用逗號分割，header表示數(shù)據(jù)文件中第一行是列名，stringsAsFactors表示讀入的數(shù)據(jù)，對于字符類型，是否自動處理為因子類型。為了方便后邊模型處理，這里設置為因子類型。
str()函數(shù)可以對數(shù)據(jù)集有一個匯總。從中看可以看出，AnimalID-個體編號，SireID-父本編號，DamID-母本編號，F(xiàn)amilyID-家系編號，SexID-性別，TankID-測試池號等字符類型，已經(jīng)被設置為因子類型了。M1BW-入池前體重，M2BW-收獲體重和M2Age-收獲時日齡均為數(shù)字變量。
注意，下邊代碼，shrimp.csv文件保存在了datasets文件夾下。該文件夾與R腳本文件在同一目錄下。
library(data.table)數(shù)據(jù)文件準備完畢，我們正式開始學習線性混合效應模型。
2 引子：線性混合效應模型可以做什么 我們從一個簡單的問題開始。對蝦是有性別的，分雌雄。如果你對對蝦沒有任何了解，你可能會想知道，雄蝦和雌蝦的體重差別大嗎？我們測定了4282尾蝦的體重。
不同性別的比較：
我們首先直觀的看一下雌雄蝦體重分布點圖。
library(ggplot2)
★ 圖中，可以看出，雌性蝦的平均體重 大于 雄性蝦的平均體重
” shrimp.sex.m2bw <- shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID)]不同池塘的比較：
分析一下數(shù)據(jù)，你會發(fā)現(xiàn)，雌雄蝦分布在2個測試池中，且來自105個家系，并且每個家系的入池體重是存在差別的。而且你會發(fā)現(xiàn)，兩個池子的養(yǎng)殖管理水平存在較大差別。
ggplot(shrimp,aes(x=TankID,y=M2BW,color=TankID))+geom_boxplot()+圖2 ★ 圖中可以看出，池塘1的平均產(chǎn)量要大于池塘2的平均產(chǎn)量
” 我們需要去除測試池、家系對雌雄蝦體重的影響，準確估計雌雄蝦體重。接下來就要用到線性模型了。
★ 鄧飛注：影響體重有很多因素，包括性別，體重，家系等因素，如何判斷哪一個蝦的體重真的好，需要使用模型進行分析。
” 3 線性混合效應模型簡介 模型1
表示一尾蝦的體重由性別和隨機誤差決定。其中Sex作為固定效應，ε作為隨機效應。后者表示所有影響體重的不可測量的效應總和，是隨機和不可控制的。
從數(shù)據(jù)中我們發(fā)現(xiàn)，一尾蝦的體重還受它所在的測試池和所在家系的影響。因此，這兩個效應也需要放到模型中。模型進一步變?yōu)椋?/p>
模型2
模型3
新加入的兩個變量，Tank和Family，如果都作為固定效應。那么上述模型稱為線性模型。如果Family作為隨機效應，那么模型3稱為線性混合效應模型（固定效應+隨機效應）。
固定因子 VS 隨機因子 這里碰到的一個棘手問題是，模型中一個效應到底是作為固定效應，還是隨機效應？準確的說，應該是與研究目的相關。
SAS for Mixed models (Second edition) 手冊中對固定效應的定義為：“An effect is called fixed if the levels in the study represent all possible levels of the factor, or at least all levels about which inference is to be made”?？珊唵蔚乩斫鉃椤霸撔乃兴皆趯嶒炄后w中都已經(jīng)出現(xiàn)”。譬如在本數(shù)據(jù)集中，性別只有雌、雄兩個水平，因此模型中性別一般作為固定效應。再比如，測試投喂5種飼料對對蝦體重的影響。由于目的很明確，只是評估這5種飼料的差異，因此飼料應作為固定效應。
隨機效應的定義為：“Factor effects are random if they are used in the study to represent only a sample (ideally, a random sample) of a larger set of potential levels”?？珊唵蔚乩斫鉃椤霸囼炄后w出現(xiàn)的該效應的水平只是一個很大水平數(shù)中的隨機抽樣。
固定效應和隨機效應的差別在哪里？“一個效應作為固定效應，還是隨機效應，應該依據(jù)研究的目的而定”?！癮 factor is considered random if its levels plausibly represent a larger population with a probability distribution”。如果我們分析一個效應的目的是為了研究它所在的一個具有概率分布的大群體的情況，那么這個效應應該作為隨機效應。隨機效應有兩個特點，a) 它是大群體中的一個樣本，b) 它具有概率分布。
譬如在shrimp數(shù)據(jù)集中，我們當前的目的是分析雌雄兩個性別的體重差異，那么105個家系就是很大家系中的一個小樣本，因此作為隨機效應更為合適。
4 線型混合效應模型R實戰(zhàn)分析 4.1 簡單線性模型 lm()是R自帶的函數(shù)。summary()函數(shù)輸出shrimp.lm的結(jié)果。
shrimp.lm <- lm(M2BW~SexID,shrimp)現(xiàn)在解釋上邊的結(jié)果。首先來看Multiple R-squared: 0.2704，它表示模型對總體方差的解釋能力。具體意思可以解釋為，總體方差中的27.04%，可以由這個模型來解釋。Adjusted R-squared是對Multiple R-squared的矯正，主要是考慮了固定效應。固定效應越多，該值越低。
下一個概念是非常的重要，那就是p值。P值是否小于0.05或者0.001，已經(jīng)成為文章結(jié)果是否可靠，是否能夠發(fā)表的一個重要標志。但是，p值一定程度上就像SCI論文的影響因子，有點濫用的味道。Nature上關于p值的故事:中文版本；英文版本。
★ 鄧飛注：中文版文章的已經(jīng)掛了，英文版的文章：https://www./news/scientific-method-statistical-errors-1.14700
” 我們還本溯源，本質(zhì)上p值是一個條件概率，它表示無效假設為真時的概率。那什么是無效假設呢？在上文中，無效假設是“對于體重性狀，雌雄蝦間沒有差異”，也就是說雌雄蝦體重是相等的。在本文中，p值為<0.001 ，意味著如果”雌雄間沒有差異，那么數(shù)據(jù)基本上不可能是這樣“，因為雌雄間沒有差異的概率太低！這反過來表明，性別影響了對蝦的體重，雌雄蝦體重是有差別的，也就是說統(tǒng)計上是顯著的。
需要注意的另外一個問題是，模型所有效應的顯著性（最底部）與系數(shù)列表單個效應的顯著性還是存在差別的。上文中，兩個結(jié)果是一致的，主要是因為模型中只包括1個固定效應。如果有更多效應，這兩個值就不再相等。
F值是我們應該關注的另外一個參數(shù)，表示模型是否顯著的一個重要參數(shù)。F值可以簡單理解為處理方差與誤差方差的比值，譬如在上文中，可以理解為性別間體重方差與殘差方差的比值，這個值越大，那么表示雌雄間體重差異越大。需要注意，在上文中F值是1928，與兩個自由度有關系（性別2-1；誤差4281-1）。
接下來重點討論系數(shù)列表。你會看到SexIDMale，你可能會問，SexID有兩個水平，F(xiàn)emale去哪里了？Estimates這一列表示的固定效應值到底是什么意思？
需要注意，系數(shù)列表中最后一列p值，表示估計值偏離0的程度。
首先看一下系數(shù)列表中的(Intercept) 項，估計值是34.376，是不是感覺很熟悉？它是Female體重的均值。在本文前邊我們估計了雌雄體重的均值。
再看一下SexIDMale的估計值，是-6.213。Intercept估計值+SexIDMale估計值=雄蝦體重值28.162。為什么要用雌蝦體重值作為截距？為什么雄蝦體重固定效應估計值要表示為與截距，也就是雌蝦體重值的差值？
第一個問題，之所以選擇雌蝦而不是雄蝦作為截距，只是因為Female和Male這兩個水平Female根據(jù)字母排序在Male前邊。
第二個問題，之所以把雄蝦體重固定效應值表示為與截距之差，是為了與協(xié)變量等統(tǒng)一起來。為了解釋，大家先看圖。
sex.bw <- shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID)]
★ 鄧飛注：以前我只知道，固定因子的第一個Level為強制為0，原來原因是這樣的。
” 我們可以根據(jù)模型1，預測雌雄蝦體重：
我們隨機選取一尾蝦，如果是雄蝦，那么斜率為-6.2137345，雄蝦的預測體重為：
我們隨機選取一尾蝦，如果是雌蝦，因為截距本身設置為雌蝦體重均值，因此斜率為0，雌蝦的預測體重為：
我們繼續(xù)來看一個稍微復雜的例子，在模型中加入TankID效應
shrimp.lm.sex.tank <- lm(M2BW~SexID+TankID,shrimp)Call:從上邊模型中可以看出，截距是35.8560193。我們把數(shù)據(jù)匯總一下：
> shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID,TankID)]我們發(fā)現(xiàn)，模型實際上是把Female.T1這個組合均值設置為截距（群體均值）。Male到Female.T1的斜率為-6.2081881，T2到Female.T1的斜率為-2.9644499；同理，F(xiàn)emale到Female.T1的斜率為0，T1到Female.T1的斜率也為0。
那么，根據(jù)模型2,可以預測雌雄蝦體重：
如果隨機從T1池撈取一尾雄蝦，它的預測體重為：
這個數(shù)值跟Male.T1的均值非常接近。
如果隨機從T2池撈取一尾雄蝦，它的預測體重為
這個數(shù)值跟Male.T2的均值非常接近。
如果隨機從T2池撈取一尾雌蝦，它的預測體重為
這個數(shù)值跟Female.T2的均值非常接近。
4.2 帶協(xié)變量的線性模型 什么是協(xié)變量（covariate)？wikipedia中是這樣解釋協(xié)變量的：In statistics, a covariate is a variable that is possibly predictive of the outcome under study. A covariate may be of direct interest or it may be a confounding or interacting variable. The alternative terms explanatory variable, independent variable, or predictor, are used in a regression analysis.
舉個簡單例子，我們從每個家系選擇選擇30尾個體，尾部用VIE顏色標記后，混合養(yǎng)殖在一起，比較每個家系的生長性能。需要考慮的一個問題就是，混合養(yǎng)殖前每個家系的大小，是不一樣的。一般，混養(yǎng)前體重越大的家系，混養(yǎng)結(jié)束時體重也會越大。我們在分析混養(yǎng)后體重時，需要考慮混養(yǎng)前每個家系體重對它的影響。這個混養(yǎng)前家系體重，就是一個協(xié)變量，在回歸分析中，也會稱為自變量，解釋變量。
計算每個家系的混養(yǎng)后平均收獲體重（M2BW）
> shrimp.taggingbw.per.family <- shrimp[,lapply(.SD,mean,na.rm=TRUE),by=.(FamilyID),.SDcols=c("M1BW","M2BW")]繪制出混養(yǎng)前家系初始體重（M1BW）和混養(yǎng)后體重的散點圖（M2BW）：
ggplot(data=shrimp.taggingbw.per.family,aes(x=M1BW,y=M2BW))+
分析下邊2個模型：
模型4
> shrimp.lm.m1bw <- lm(M2BW ~ M1BW,shrimp)上邊實際上是簡單的回歸分析，系數(shù)表中Intercept為24.4789861，表示當M1BW為0時，混養(yǎng)體重M2BW為24.4789861g。實際上初始體重不可能為0。因此這個截距沒有實際意義。
M1BW的系數(shù)估計值為0.9484804，表示初始體重每增加1g，混養(yǎng)體重會增加0.9484804g。
模型5
> shrimp.lm.sex.m1bw <- lm(M2BW~SexID+SexID:M1BW,shrimp)這個模型稍微復雜一些。從系數(shù)列表中，針對雌性兩個性別，給出了不同的回歸系數(shù)。這主要是由于雌雄生長速度的差異造成的，后期雌蝦生長速要快于雄蝦。
ggplot(data=shrimp,aes(x=M1BW,y=M2BW,color=SexID))+geom_point()圖5 我們根據(jù)模型5來預測一尾雄蝦的體重,假定它的混養(yǎng)前體重為5g，那么它的混養(yǎng)后體重為26.0340769g。
4.3 包括交互效應的線性模型 模型已經(jīng)變的越來越復雜了。我們進一步在模型5中加入養(yǎng)殖池TankID效應。
模型6
> # 模型6從Adjusted R-squared參數(shù)（0.3798813 vs 0.2701994可以看出，模型6 對數(shù)據(jù)的解釋能力，相比模型1，已經(jīng)提高了。也就是說根據(jù)模型6預測體重，準確性比模型1要高。
什么是交互效應？Interaction effects represent the combined effects of factors on the dependent measure. When an interaction effect is present, the impact of one factor depends on the level of the other factor。兩個效應間相互影響，促進或者制約的關系。
譬如，我們數(shù)據(jù)中有Sex和Tank兩個固定效應，那么我們可能會想 雌蝦會不會在特別偏愛某種環(huán)境，譬如在T1池中長得比T2池中大，但是雄蝦可能會在T2池中長的比T1池大？因此我們可以在模型中試著加入交互效應。
模型7
> # 模型7從系數(shù)表中，可以看出，相比兩個主效應（SexIDMale：-5.563925、TankIDT2：-3.0719274），Sex和Tank的交互效應非常?。?.2507485），p > 0.05，統(tǒng)計檢驗與0相比，達不到顯著水平。Adjusted R-squared: 0.3798447也表明，相比模型6，模型7的解釋能力并沒有提高。
我們圖形展示一下模型7固定效應：
library(sjPlot)
4.4 包括隨機效應的線性混合效應模型 請加載一個新的數(shù)據(jù)集shrimpex.csv，其中有一個PopID字段，包括Pop1到Pop4共計4個水平，表示shrimp數(shù)據(jù)由四個群體組成。現(xiàn)在考慮這樣一個問題：四個群體間收獲體重是否存在差異。
首先加載數(shù)據(jù)文件。畫出四個群體收獲體重的箱形圖，加上jitter點。
ggplot(data=shrimp,aes(x=PopID,y=M2BW,color=PopID)) +從上圖中，大致可以看出，群體間是存在差異的。
進一步分析數(shù)據(jù)，你會發(fā)現(xiàn)每個群體由多個家系組成，見下圖。
ggplot(data=shrimp,aes(x=PopID,y=M2BW,fill=FamilyID)) +
理論上，我們從每個群體抽樣時，抽樣個體是代表該群體的隨機樣本。但是，一個群體內(nèi)的個體往往存在親緣關系，譬如（全同胞、半同胞個體）。因此抽樣個體存在兩個層次：每個群體包括多個家系，每個家系包括數(shù)量不等的個體。
從上圖中可以看出，每個群體內(nèi)的不同家系間是存在差異的。
每一個家系內(nèi)的個體，遺傳自同一對親本，相互間相似性更強。不同家系個體的體重均值是不一樣的。
這實際上違背了樣本觀察值的獨立性原則，同一個家系內(nèi)的全同胞個體的體重值實際上是由他們親本所決定。
針對這種情況，我們把家系效應作為隨機效應加入模型中。這相當于，給每個家系設置一個基線，類似于不同的家系有不同的平均體重，也稱作隨機截距（random intercept）。
在模型中加入家系隨機效應，那么觀察值的非獨立性問題就解決了。
為了說明家系結(jié)構(gòu)對分析結(jié)果的影響，故意在每個群體中設置了一個均值特別高的家系。在實際測試數(shù)據(jù)中，這種現(xiàn)象也會經(jīng)常出現(xiàn)。如果我們分析時不考慮群體內(nèi)的家系結(jié)構(gòu)，那么家系方差會被累加到殘差方差中。
如果采取方差分析的方法，你也會發(fā)現(xiàn)，忽略家系結(jié)構(gòu)，群體的均方值可能會非常大，被嚴重高估。
根據(jù)教程1對于固定效應和隨機效應的討論，由于我們的目的是要分析四個群體間的差異，獲得每個群體的性能，因此群體更適合做固定效應。每個群體是由多個家系組成的，這些家系只是大量家系的一個隨機抽樣，因此更加適合作為隨機效應。
下邊我們通過兩個模型實例，來看一下家系結(jié)構(gòu)對分析結(jié)果的影響。
模型8
模型8不考慮家系結(jié)構(gòu)，Pop、Sex和Tank為固定效應，Sex:M1BW為協(xié)變量。
分析結(jié)果如下：
library(lmerTest)模型9 Pop：Family為隨機效應
Pop，Sex和Tank為固定效應
Sex：M1BW為協(xié)變量
在模型中加入隨機效應，需要使用lme4包中的lmer函數(shù)。下邊代碼中的(1|PopID:FamilyID)，表示針對不同的家系，單獨估計其隨機截距（random intercept）。其中1表示隨機截距。
分析結(jié)果如下：
> library(lme4)> anova(shrimp.lm.9)把模型8的Residual standard error的平方，與模型9 Random Effects部分對比，你會發(fā)現(xiàn)，如果不考慮家系結(jié)構(gòu)，殘差方差明顯被高估,估計值為18.4559804 ?？紤]家系結(jié)構(gòu)后，殘差方差為12.5267429, 明顯變小， 從殘差中分離出了大部分的家系方差。由家系隨機效應產(chǎn)生的方差，估計值為5.9328232，占表型方差的32%。
從anova方差分析(如果對基礎概念不了解，可參考這篇blog：自由度和方差分析)的角度看，加入家系隨機效應后，群體固定效應（PopID）盡管仍然也達到了顯著水平，但是均方和F值明顯變小。這表明存在這樣一種風險，如果考慮群體內(nèi)的家系結(jié)構(gòu)，本來兩個群體的差異可能達不到顯著水平，但是如果忽視了這種家系結(jié)構(gòu)，兩個群體間的差異統(tǒng)計上會表現(xiàn)為顯著水平，從而誤判群體間的實際性能差別。
我們看一下，基于模型9（不包括家系的隨機效應），預測四個群體家系的性能，如下圖所示：你會發(fā)現(xiàn)，每個群體中特別大的家系效應，已經(jīng)被剔除掉了。
ps：擬合值反應的是包括所有固定和隨機效應的結(jié)果，lmer中通過fitted()函數(shù)獲得該值。預測值，是可以設定不包括隨機效應的，lmer中通過predict()函數(shù)獲得該值。
shrimp.lm.9.predict <- predict(shrimp.lm.9,re.form=NA) #擬合值4.5 隨機效應的顯著性檢驗 模型中加入隨機效應后，如何檢驗其顯著性？
一般是通過似然比率檢驗（Likelihood ratio test, LRT)來實現(xiàn)。lmerTest包中的rand()函數(shù)提供該功能。
首先我們來看一下似然比率檢驗的結(jié)果：
library(lmerTest)> ranova(shrimp.lm.9)4.6 獲得每個群體的性能 調(diào)用emmeans包中的函數(shù)，計算四個群體的估計邊際均值(estimated marginal means)，或者說最小二乘均值(least-squares means)。根據(jù)邊際均值，我們可以對群體的性能進行排序和比較。
關于emmeans包，請參考日志最小二乘均值的估計模型。盡管該日志介紹的是lsmeans包，但用法跟emmeans包都是一樣的。而且根據(jù)作者介紹，在不久的將來，emmeans包要替代lsmeans包。
注意，安裝emmeans還需要pbkrtest包，這個包沒有自動安裝，需要手動安裝。
> library(emmeans)從上邊結(jié)果中查找emmean列，可以看到Pop1群體的邊際均值最大，這表明四個群體中該群體性能最好。
從圖1中可以大體看出，雌蝦體重比雄蝦高。然后，我們實際計算雌雄蝦體重均值，發(fā)現(xiàn)雌蝦的確比雄蝦重。
因為Female和Male均為因子變量，因此在x軸上可以將Female標準化為0，Male與Female的間距為1，二者體重差值為-6.2137345，那么斜線的斜率可以認為等于-6.213。這樣就可以很方便的求解或者說預測體重。
很明顯，二者表現(xiàn)為正相關，相關系數(shù)為0.39。因此在對比雌雄體重差異時，需要考慮加入M1BW作為協(xié)變量，對M2BW進行校正。
不考慮除殘差外的隨機效應，目前模型6是最優(yōu)模型。我們根據(jù)模型6，可以回答最初的問題，雌雄體重間差異顯著。接下來，我們會考慮在模型中加入隨機效應，進入線性混合效應模型部分。
這里遇到了一個問題，在評價群體間的差異時，是否需要考慮每個群體內(nèi)的家系結(jié)構(gòu)？
