栾生老师 || 线性混合效应模型教程
★邓飞注：混合线性模型，我以前写过几篇学习笔记（线性混合模型系列1~5），但是偏向基础知识，看过栾老师的博客后，感觉思路清晰了很多。经过栾老师许可，放到公众号上，我对原来的代码进行了一些小的更改（有些包的函数变化了）。非常认可栾老师的观点：“学习混合线性模型，可以对比一般线性模型，一边学习理论，一边动手实践”。点击阅读原文，查看栾老师的博客。另外，栾老师用
ggplot2对数据和模型进行的可视化分析，非常赞。 ”
原文包括两部分，这里我合二为一。
孪生老师简介：
以下为正文
学习线性混合效应模型（Linear Mixed Effects Model，LMM）最好的方法，是一边学习理论，一边动手实践，这样印象最为深刻。本文参考了Bodo Winter博士的教程Linear models and linear mixed effects models in R教程1教程2的结构。
★教程1 pdf链接：http://www.bodowinter.com/tutorial/bw_LME_tutorial1.pdf 教程2 pdf链接：http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf ”
本文中，为了便于理解，使用的数据集是来自csv文件shrimp中的对虾育种数据（对原始数据已经进行了变换）。
★邓飞注：原始数据下载链接，https://luansheng.netlify.com/post/datasets/shrimp.csv ”
推荐使用Rstudio来运行R，依赖的R包有：
- data.table
- ggplot2
- lme4
- sjPlot
- emmeans
- lmerTest
★邓飞注：这些包，data.table是读取写入数据，ggplot2作图包，lme4混合线性包，sjPlot模型作图包，emmeans计算预测均值，lmerTest是固定因子和随机因子显著性检验。安装包：install.packages("data.table") 命令按照 data.table 包，以此类推，这些包都是在CRAN中的，直接使用
install.packages进行安装。 ”
1 读取数据文件
fread()函数来自data.table包，特点是对大文件读取速度特别快。建议使用data.table作为数据处理的主力包。
其中的参数sep表示文件是用逗号分割，header表示数据文件中第一行是列名，stringsAsFactors表示读入的数据，对于字符类型，是否自动处理为因子类型。为了方便后边模型处理，这里设置为因子类型。
str()函数可以对数据集有一个汇总。从中看可以看出，AnimalID-个体编号，SireID-父本编号，DamID-母本编号，FamilyID-家系编号，SexID-性别，TankID-测试池号等字符类型，已经被设置为因子类型了。M1BW-入池前体重，M2BW-收获体重和M2Age-收获时日龄均为数字变量。
注意，下边代码，shrimp.csv文件保存在了datasets文件夹下。该文件夹与R脚本文件在同一目录下。
library(data.table)
shrimp <- fread(input = "shrimp.csv",sep = ",",header = TRUE,stringsAsFactors = TRUE)
> str(shrimp)
Classes ‘data.table’ and 'data.frame': 4282 obs. of 10 variables:
$ AnimalID: Factor w/ 4282 levels "13G1000001","13G1000002",..: 3308 3307 2215 1303 3601 2184 2194 2175 1585 2176 ...
$ SireID : Factor w/ 100 levels "12G000K010","12G000K065",..: 81 81 81 81 81 81 81 81 81 81 ...
$ DamID : Factor w/ 91 levels "12G000K052","12G000K097",..: 81 81 81 81 81 81 81 81 81 81 ...
$ PopID : Factor w/ 4 levels "Pop1","Pop2",..: 1 1 1 1 1 1 1 1 1 1 ...
$ FamilyID: Factor w/ 105 levels "13F1306003","13F1306004",..: 6 6 6 6 6 6 6 6 6 6 ...
$ SexID : Factor w/ 2 levels "Female","Male": 2 2 1 1 1 2 1 2 2 1 ...
$ TankID : Factor w/ 2 levels "T1","T2": 1 1 1 1 1 1 1 1 1 1 ...
$ M1BW : num 8.13 8.13 8.13 8.13 8.13 8.55 8.55 8.55 8.55 8.55 ...
$ M2BW : num 29 30.5 33.3 40.1 43 29.1 30.7 30.7 32.5 35.6 ...
$ M2Age : int 219 219 219 219 219 219 219 219 219 219 ...
- attr(*, ".internal.selfref")=<externalptr>
数据文件准备完毕，我们正式开始学习线性混合效应模型。
2 引子：线性混合效应模型可以做什么
我们从一个简单的问题开始。对虾是有性别的，分雌雄。如果你对对虾没有任何了解，你可能会想知道，雄虾和雌虾的体重差别大吗？我们测定了4282尾虾的体重。
不同性别的比较：
我们首先直观的看一下雌雄虾体重分布点图。
library(ggplot2)
ggplot(data=shrimp,aes(x=SexID,y=M2BW,color=SexID))+
geom_boxplot()+geom_dotplot(binaxis = "y",
stackdir = "center",position = "dodge",binwidth = 0.25)
从图1中可以大体看出，雌虾体重比雄虾高。然后，我们实际计算雌雄虾体重均值，发现雌虾的确比雄虾重。
★图中，可以看出，雌性虾的平均体重 大于 雄性虾的平均体重 ”
shrimp.sex.m2bw <- shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID)]
> shrimp.sex.m2bw
SexID V1
1: Male 28.16273
2: Female 34.37646
不同池塘的比较：
分析一下数据，你会发现，雌雄虾分布在2个测试池中，且来自105个家系，并且每个家系的入池体重是存在差别的。而且你会发现，两个池子的养殖管理水平存在较大差别。
ggplot(shrimp,aes(x=TankID,y=M2BW,color=TankID))+geom_boxplot()+
geom_dotplot(binaxis = "y",stackdir = "center",position = "dodge",binwidth = 0.3)
图2
★图中可以看出，池塘1的平均产量要大于池塘2的平均产量 ”
我们需要去除测试池、家系对雌雄虾体重的影响，准确估计雌雄虾体重。接下来就要用到线性模型了。
★邓飞注：影响体重有很多因素，包括性别，体重，家系等因素，如何判断哪一个虾的体重真的好，需要使用模型进行分析。 ”
3 线性混合效应模型简介
模型1
表示一尾虾的体重由性别和随机误差决定。其中Sex作为固定效应，ε作为随机效应。后者表示所有影响体重的不可测量的效应总和，是随机和不可控制的。
从数据中我们发现，一尾虾的体重还受它所在的测试池和所在家系的影响。因此，这两个效应也需要放到模型中。模型进一步变为：
模型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”。可简单地理解为“该效应的所有水平在实验群体中都已经出现”。譬如在本数据集中，性别只有雌、雄两个水平，因此模型中性别一般作为固定效应。再比如，测试投喂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”。可简单地理解为“试验群体出现的该效应的水平只是一个很大水平数中的随机抽样。
固定效应和随机效应的差别在哪里？“一个效应作为固定效应，还是随机效应，应该依据研究的目的而定”。“a factor is considered random if its levels plausibly represent a larger population with a probability distribution”。如果我们分析一个效应的目的是为了研究它所在的一个具有概率分布的大群体的情况，那么这个效应应该作为随机效应。随机效应有两个特点，a) 它是大群体中的一个样本，b) 它具有概率分布。
譬如在shrimp数据集中，我们当前的目的是分析雌雄两个性别的体重差异，那么105个家系就是很大家系中的一个小样本，因此作为随机效应更为合适。
4 线型混合效应模型R实战分析
4.1 简单线性模型
lm()是R自带的函数。summary()函数输出shrimp.lm的结果。
shrimp.lm <- lm(M2BW~SexID,shrimp)
shrimp.lm.sum <- summary(shrimp.lm)
> shrimp.lm.sum
Call:
lm(formula = M2BW ~ SexID, data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-22.9765 -3.1765 -0.0765 3.0235 18.9235
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.3765 0.1082 317.65 <2e-16 ***
SexIDMale -6.2137 0.1560 -39.82 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.101 on 4280 degrees of freedom
Multiple R-squared: 0.2704, Adjusted R-squared: 0.2702
F-statistic: 1586 on 1 and 4280 DF, p-value: < 2.2e-16
现在解释上边的结果。首先来看Multiple R-squared: 0.2704，它表示模型对总体方差的解释能力。具体意思可以解释为，总体方差中的27.04%，可以由这个模型来解释。Adjusted R-squared是对Multiple R-squared的矫正，主要是考虑了固定效应。固定效应越多，该值越低。
下一个概念是非常的重要，那就是p值。P值是否小于0.05或者0.001，已经成为文章结果是否可靠，是否能够发表的一个重要标志。但是，p值一定程度上就像SCI论文的影响因子，有点滥用的味道。Nature上关于p值的故事:中文版本；英文版本。
★邓飞注：中文版文章的已经挂了，英文版的文章：https://www.nature.com/news/scientific-method-statistical-errors-1.14700 ”
我们还本溯源，本质上p值是一个条件概率，它表示无效假设为真时的概率。那什么是无效假设呢？在上文中，无效假设是“对于体重性状，雌雄虾间没有差异”，也就是说雌雄虾体重是相等的。在本文中，p值为<0.001 ，意味着如果”雌雄间没有差异，那么数据基本上不可能是这样“，因为雌雄间没有差异的概率太低！这反过来表明，性别影响了对虾的体重，雌雄虾体重是有差别的，也就是说统计上是显著的。
需要注意的另外一个问题是，模型所有效应的显著性（最底部）与系数列表单个效应的显著性还是存在差别的。上文中，两个结果是一致的，主要是因为模型中只包括1个固定效应。如果有更多效应，这两个值就不再相等。
F值是我们应该关注的另外一个参数，表示模型是否显著的一个重要参数。F值可以简单理解为处理方差与误差方差的比值，譬如在上文中，可以理解为性别间体重方差与残差方差的比值，这个值越大，那么表示雌雄间体重差异越大。需要注意，在上文中F值是1928，与两个自由度有关系（性别2-1；误差4281-1）。
接下来重点讨论系数列表。你会看到SexIDMale，你可能会问，SexID有两个水平，Female去哪里了？Estimates这一列表示的固定效应值到底是什么意思？
需要注意，系数列表中最后一列p值，表示估计值偏离0的程度。
首先看一下系数列表中的(Intercept) 项，估计值是34.376，是不是感觉很熟悉？它是Female体重的均值。在本文前边我们估计了雌雄体重的均值。
再看一下SexIDMale的估计值，是-6.213。Intercept估计值+SexIDMale估计值=雄虾体重值28.162。为什么要用雌虾体重值作为截距？为什么雄虾体重固定效应估计值要表示为与截距，也就是雌虾体重值的差值？
第一个问题，之所以选择雌虾而不是雄虾作为截距，只是因为Female和Male这两个水平Female根据字母排序在Male前边。
第二个问题，之所以把雄虾体重固定效应值表示为与截距之差，是为了与协变量等统一起来。为了解释，大家先看图。
sex.bw <- shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID)]
ggplot(data=shrimp,aes(x=SexID,y=M2BW,color=SexID))+
geom_dotplot(binaxis = "y",stackdir = "center",position = "dodge",binwidth = 0.25)+
geom_vline(xintercept = 1)+
geom_vline(xintercept = 2)+
geom_hline(yintercept = 33.93992)+
geom_hline(yintercept = 27.77438)+
geom_abline(intercept = 40.75,slope=-6.6155)+
geom_text(x=2.4,y=31.5,label="-6.6155")+
annotate("segment",x=2.5,xend=2.5,y=27.77438,yend = 33.93992)
因为Female和Male均为因子变量，因此在x轴上可以将Female标准化为0，Male与Female的间距为1，二者体重差值为-6.2137345，那么斜线的斜率可以认为等于-6.213。这样就可以很方便的求解或者说预测体重。
★邓飞注：以前我只知道，固定因子的第一个Level为强制为0，原来原因是这样的。 ”
我们可以根据模型1，预测雌雄虾体重：
我们随机选取一尾虾，如果是雄虾，那么斜率为-6.2137345，雄虾的预测体重为：
我们随机选取一尾虾，如果是雌虾，因为截距本身设置为雌虾体重均值，因此斜率为0，雌虾的预测体重为：
我们继续来看一个稍微复杂的例子，在模型中加入TankID效应
shrimp.lm.sex.tank <- lm(M2BW~SexID+TankID,shrimp)
shrimp.lm.sex.tank.sum <- summary(shrimp.lm.sex.tank)
shrimp.lm.sex.tank.sum
Call:
lm(formula = M2BW ~ SexID + TankID, data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-21.4916 -3.0478 -0.1697 2.8166 20.4084
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 35.8560 0.1276 281.10 <2e-16 ***
SexIDMale -6.2082 0.1493 -41.58 <2e-16 ***
TankIDT2 -2.9644 0.1492 -19.87 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.882 on 4279 degrees of freedom
Multiple R-squared: 0.332, Adjusted R-squared: 0.3317
F-statistic: 1063 on 2 and 4279 DF, p-value: < 2.2e-16
从上边模型中可以看出，截距是35.8560193。我们把数据汇总一下：
> shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID,TankID)]
SexID TankID V1
1: Male T1 29.62356
2: Female T1 35.87844
3: Male T2 26.70756
4: Female T2 32.86907
我们发现，模型实际上是把Female.T1这个组合均值设置为截距（群体均值）。Male到Female.T1的斜率为-6.2081881，T2到Female.T1的斜率为-2.9644499；同理，Female到Female.T1的斜率为0，T1到Female.T1的斜率也为0。
那么，根据模型2,可以预测雌雄虾体重：
如果随机从T1池捞取一尾雄虾，它的预测体重为：
这个数值跟Male.T1的均值非常接近。
如果随机从T2池捞取一尾雄虾，它的预测体重为
这个数值跟Male.T2的均值非常接近。
如果随机从T2池捞取一尾雌虾，它的预测体重为
这个数值跟Female.T2的均值非常接近。
4.2 带协变量的线性模型
什么是协变量（covariate)？wikipedia中是这样解释协变量的：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颜色标记后，混合养殖在一起，比较每个家系的生长性能。需要考虑的一个问题就是，混合养殖前每个家系的大小，是不一样的。一般，混养前体重越大的家系，混养结束时体重也会越大。我们在分析混养后体重时，需要考虑混养前每个家系体重对它的影响。这个混养前家系体重，就是一个协变量，在回归分析中，也会称为自变量，解释变量。
计算每个家系的混养后平均收获体重（M2BW）
> shrimp.taggingbw.per.family <- shrimp[,lapply(.SD,mean,na.rm=TRUE),by=.(FamilyID),.SDcols=c("M1BW","M2BW")]
> shrimp.taggingbw.per.family
FamilyID M1BW M2BW
1: 13F1306014 8.503478 32.73478
2: 13F1306015 8.614000 36.00000
3: 13F1306016 8.753810 34.89286
4: 13F1306019 7.747561 34.25122
5: 13F1306022 8.931842 33.66842
---
101: 13F1306802 7.274138 27.69655
102: 13F1367013 5.438571 25.70000
103: 13F1367075 5.016250 26.19750
104: 13F1367097 6.937500 26.18125
105: 13F1367106 9.000222 38.87778
绘制出混养前家系初始体重（M1BW）和混养后体重的散点图（M2BW）：
ggplot(data=shrimp.taggingbw.per.family,aes(x=M1BW,y=M2BW))+
geom_point()+
geom_smooth(method = "lm")
很明显，二者表现为正相关，相关系数为0.39。因此在对比雌雄体重差异时，需要考虑加入M1BW作为协变量，对M2BW进行校正。
分析下边2个模型：
模型4
> shrimp.lm.m1bw <- lm(M2BW ~ M1BW,shrimp)
> summary(shrimp.lm.m1bw)
Call:
lm(formula = M2BW ~ M1BW, data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-20.6953 -4.1186 -0.3901 3.8326 21.3708
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.47899 0.47668 51.35 <2e-16 ***
M1BW 0.94848 0.06434 14.74 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.832 on 4239 degrees of freedom
(41 observations deleted due to missingness)
Multiple R-squared: 0.04876, Adjusted R-squared: 0.04854
F-statistic: 217.3 on 1 and 4239 DF, p-value: < 2.2e-16
上边实际上是简单的回归分析，系数表中Intercept为24.4789861，表示当M1BW为0时，混养体重M2BW为24.4789861g。实际上初始体重不可能为0。因此这个截距没有实际意义。
M1BW的系数估计值为0.9484804，表示初始体重每增加1g，混养体重会增加0.9484804g。
模型5
> shrimp.lm.sex.m1bw <- lm(M2BW~SexID+SexID:M1BW,shrimp)
> summary(shrimp.lm.sex.m1bw)
Call:
lm(formula = M2BW ~ SexID + SexID:M1BW, data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-23.719 -3.078 -0.090 2.846 18.417
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 27.28624 0.54727 49.858 < 2e-16 ***
SexIDMale -5.91782 0.80948 -7.311 3.16e-13 ***
SexIDFemale:M1BW 0.97548 0.07385 13.208 < 2e-16 ***
SexIDMale:M1BW 0.93313 0.08053 11.587 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.933 on 4237 degrees of freedom
(41 observations deleted due to missingness)
Multiple R-squared: 0.3196, Adjusted R-squared: 0.3191
F-statistic: 663.4 on 3 and 4237 DF, p-value: < 2.2e-16
这个模型稍微复杂一些。从系数列表中，针对雌性两个性别，给出了不同的回归系数。这主要是由于雌雄生长速度的差异造成的，后期雌虾生长速要快于雄虾。
ggplot(data=shrimp,aes(x=M1BW,y=M2BW,color=SexID))+geom_point()
图5
我们根据模型5来预测一尾雄虾的体重,假定它的混养前体重为5g，那么它的混养后体重为26.0340769g。
4.3 包括交互效应的线性模型
模型已经变的越来越复杂了。我们进一步在模型5中加入养殖池TankID效应。
模型6
> # 模型6
> shrimp.lm.sex.tank.m1bw <- lm(M2BW ~ SexID + TankID + SexID:M1BW,shrimp)
> shrimp.lm.sex.tank.m1bw.sum <- summary(shrimp.lm.sex.tank.m1bw)
> print(shrimp.lm.sex.tank.m1bw.sum)
Call:
lm(formula = M2BW ~ SexID + TankID + SexID:M1BW, data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-22.2635 -2.9218 -0.1317 2.7273 19.9333
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.59445 0.52621 54.340 < 2e-16 ***
SexIDMale -5.43136 0.77289 -7.027 2.44e-12 ***
TankIDT2 -2.95103 0.14467 -20.399 < 2e-16 ***
SexIDFemale:M1BW 0.99876 0.07049 14.168 < 2e-16 ***
SexIDMale:M1BW 0.88930 0.07689 11.566 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.708 on 4236 degrees of freedom
(41 observations deleted due to missingness)
Multiple R-squared: 0.3805, Adjusted R-squared: 0.3799
F-statistic: 650.4 on 4 and 4236 DF, p-value: < 2.2e-16
从Adjusted R-squared参数（0.3798813 vs 0.2701994可以看出，模型6 对数据的解释能力，相比模型1，已经提高了。也就是说根据模型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。两个效应间相互影响，促进或者制约的关系。
譬如，我们数据中有Sex和Tank两个固定效应，那么我们可能会想 雌虾会不会在特别偏爱某种环境，譬如在T1池中长得比T2池中大，但是雄虾可能会在T2池中长的比T1池大？因此我们可以在模型中试着加入交互效应。
模型7
> # 模型7
> shrimp.lm.sex.tank.m1bw.interaction <- lm(M2BW ~ SexID + TankID + SexID:TankID + SexID:M1BW,shrimp)
> shrimp.lm.sex.tank.m1bw.interaction.sum <- summary(shrimp.lm.sex.tank.m1bw.interaction)
> print(shrimp.lm.sex.tank.m1bw.interaction.sum)
Call:
lm(formula = M2BW ~ SexID + TankID + SexID:TankID + SexID:M1BW,
data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-22.2039 -2.8961 -0.1335 2.7157 19.9954
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.64805 0.52985 54.068 < 2e-16 ***
SexIDMale -5.56393 0.78792 -7.062 1.92e-12 ***
TankIDT2 -3.07193 0.20104 -15.280 < 2e-16 ***
SexIDMale:TankIDT2 0.25075 0.28953 0.866 0.387
SexIDFemale:M1BW 0.99972 0.07050 14.180 < 2e-16 ***
SexIDMale:M1BW 0.89123 0.07692 11.586 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.708 on 4235 degrees of freedom
(41 observations deleted due to missingness)
Multiple R-squared: 0.3806, Adjusted R-squared: 0.3798
F-statistic: 520.4 on 5 and 4235 DF, p-value: < 2.2e-16
从系数表中，可以看出，相比两个主效应（SexIDMale：-5.563925、TankIDT2：-3.0719274），Sex和Tank的交互效应非常小（0.2507485），p > 0.05，统计检验与0相比，达不到显著水平。Adjusted R-squared: 0.3798447也表明，相比模型6，模型7的解释能力并没有提高。
我们图形展示一下模型7固定效应：
library(sjPlot)
plot_model(shrimp.lm.sex.tank.m1bw.interaction,show.values = TRUE)
不考虑除残差外的随机效应，目前模型6是最优模型。我们根据模型6，可以回答最初的问题，雌雄体重间差异显著。接下来，我们会考虑在模型中加入随机效应，进入线性混合效应模型部分。
4.4 包括随机效应的线性混合效应模型
请加载一个新的数据集shrimpex.csv，其中有一个PopID字段，包括Pop1到Pop4共计4个水平，表示shrimp数据由四个群体组成。现在考虑这样一个问题：四个群体间收获体重是否存在差异。
首先加载数据文件。画出四个群体收获体重的箱形图，加上jitter点。
ggplot(data=shrimp,aes(x=PopID,y=M2BW,color=PopID)) +
geom_boxplot(outlier.size = 0)+
geom_jitter(alpha=0.3)+
labs(x="群体",y="收获体重")+
theme_gray(base_size = 20)+
theme(legend.position = "none")
从上图中，大致可以看出，群体间是存在差异的。
进一步分析数据，你会发现每个群体由多个家系组成，见下图。
ggplot(data=shrimp,aes(x=PopID,y=M2BW,fill=FamilyID)) +
geom_boxplot(outlier.size = 0)+
labs(x="群体",y="收获体重")+
theme_gray(base_size = 20)+
theme(legend.position = "none")
这里遇到了一个问题，在评价群体间的差异时，是否需要考虑每个群体内的家系结构？
理论上，我们从每个群体抽样时，抽样个体是代表该群体的随机样本。但是，一个群体内的个体往往存在亲缘关系，譬如（全同胞、半同胞个体）。因此抽样个体存在两个层次：每个群体包括多个家系，每个家系包括数量不等的个体。
从上图中可以看出，每个群体内的不同家系间是存在差异的。
每一个家系内的个体，遗传自同一对亲本，相互间相似性更强。不同家系个体的体重均值是不一样的。
这实际上违背了样本观察值的独立性原则，同一个家系内的全同胞个体的体重值实际上是由他们亲本所决定。
针对这种情况，我们把家系效应作为随机效应加入模型中。这相当于，给每个家系设置一个基线，类似于不同的家系有不同的平均体重，也称作随机截距（random intercept）。
在模型中加入家系随机效应，那么观察值的非独立性问题就解决了。
为了说明家系结构对分析结果的影响，故意在每个群体中设置了一个均值特别高的家系。在实际测试数据中，这种现象也会经常出现。如果我们分析时不考虑群体内的家系结构，那么家系方差会被累加到残差方差中。
如果采取方差分析的方法，你也会发现，忽略家系结构，群体的均方值可能会非常大，被严重高估。
根据教程1对于固定效应和随机效应的讨论，由于我们的目的是要分析四个群体间的差异，获得每个群体的性能，因此群体更适合做固定效应。每个群体是由多个家系组成的，这些家系只是大量家系的一个随机抽样，因此更加适合作为随机效应。
下边我们通过两个模型实例，来看一下家系结构对分析结果的影响。
模型8
模型8不考虑家系结构，Pop、Sex和Tank为固定效应，Sex:M1BW为协变量。
分析结果如下：
library(lmerTest)
shrimp.lm.8 <- lm(M2BW ~ 1 + PopID + SexID + TankID + SexID:M1BW,shrimp)
summary(shrimp.lm.8) #加载lmerTest包后，lmer的返回结果，每个固定效应系数带有P值
Call:
lm(formula = M2BW ~ 1 + PopID + SexID + TankID + SexID:M1BW,
data = shrimp)
Residuals:
Min 1Q Median 3Q Max
-20.7242 -2.5787 -0.2092 2.2353 18.9070
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 35.42932 0.53924 65.702 < 2e-16 ***
PopIDPop2 -1.61211 0.18406 -8.759 < 2e-16 ***
PopIDPop3 -3.61817 0.18086 -20.005 < 2e-16 ***
PopIDPop4 -5.76930 0.21237 -27.166 < 2e-16 ***
SexIDMale -5.39346 0.70557 -7.644 2.59e-14 ***
TankIDT2 -2.93073 0.13206 -22.192 < 2e-16 ***
SexIDFemale:M1BW 0.40396 0.06778 5.960 2.73e-09 ***
SexIDMale:M1BW 0.30223 0.07363 4.105 4.13e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.296 on 4233 degrees of freedom
(41 observations deleted due to missingness)
Multiple R-squared: 0.4845, Adjusted R-squared: 0.4837
F-statistic: 568.4 on 7 and 4233 DF, p-value: < 2.2e-16
模型9Pop：Family为随机效应 Pop，Sex和Tank为固定效应 Sex：M1BW为协变量
在模型中加入随机效应，需要使用lme4包中的lmer函数。下边代码中的(1|PopID:FamilyID)，表示针对不同的家系，单独估计其随机截距（random intercept）。其中1表示随机截距。
分析结果如下：
> library(lme4)
> shrimp.lm.9 <- lmer(M2BW ~ 1 + PopID + SexID + TankID + SexID:M1BW + (1|PopID:FamilyID),shrimp)
> summary(shrimp.lm.9)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: M2BW ~ 1 + PopID + SexID + TankID + SexID:M1BW + (1 | PopID:FamilyID)
Data: shrimp
REML criterion at convergence: 23070.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.5304 -0.5854 0.0240 0.6388 3.8977
Random effects:
Groups Name Variance Std.Dev.
PopID:FamilyID (Intercept) 5.933 2.436
Residual 12.527 3.539
Number of obs: 4241, groups: PopID:FamilyID, 105
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 35.7188 1.0605 487.3240 33.683 < 2e-16 ***
PopIDPop2 -1.6596 0.6884 100.5550 -2.411 0.01774 *
PopIDPop3 -3.7742 0.6733 101.7519 -5.606 1.78e-07 ***
PopIDPop4 -5.8638 0.7359 110.5088 -7.968 1.58e-12 ***
SexIDMale -5.6599 0.5901 4148.3019 -9.591 < 2e-16 ***
TankIDT2 -2.9491 0.1097 4140.0606 -26.884 < 2e-16 ***
SexIDFemale:M1BW 0.3757 0.1204 992.9529 3.122 0.00185 **
SexIDMale:M1BW 0.2904 0.1231 1068.2960 2.360 0.01847 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) PpIDP2 PpIDP3 PpIDP4 SxIDMl TnIDT2 SIDF:M
PopIDPop2 -0.384
PopIDPop3 -0.419 0.506
PopIDPop4 -0.519 0.476 0.494
SexIDMale -0.248 0.002 0.007 -0.004
TankIDT2 -0.042 -0.006 -0.004 -0.001 -0.030
SxIDFm:M1BW -0.889 0.077 0.108 0.250 0.291 -0.010
SxIDMl:M1BW -0.711 0.074 0.100 0.246 -0.352 0.010 0.786
> anova(shrimp.lm.9)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
PopID 914.2 304.7 3 103.8 24.3255 5.319e-12 ***
SexID 1152.4 1152.4 1 4148.3 91.9951 < 2.2e-16 ***
TankID 9053.7 9053.7 1 4140.1 722.7464 < 2.2e-16 ***
SexID:M1BW 122.4 61.2 2 1414.2 4.8839 0.007695 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
把模型8的Residual standard error的平方，与模型9 Random Effects部分对比，你会发现，如果不考虑家系结构，残差方差明显被高估,估计值为18.4559804 。考虑家系结构后，残差方差为12.5267429, 明显变小， 从残差中分离出了大部分的家系方差。由家系随机效应产生的方差，估计值为5.9328232，占表型方差的32%。
从anova方差分析(如果对基础概念不了解，可参考这篇blog：自由度和方差分析)的角度看，加入家系随机效应后，群体固定效应（PopID）尽管仍然也达到了显著水平，但是均方和F值明显变小。这表明存在这样一种风险，如果考虑群体内的家系结构，本来两个群体的差异可能达不到显著水平，但是如果忽视了这种家系结构，两个群体间的差异统计上会表现为显著水平，从而误判群体间的实际性能差别。
我们看一下，基于模型9（不包括家系的随机效应），预测四个群体家系的性能，如下图所示：你会发现，每个群体中特别大的家系效应，已经被剔除掉了。
ps：拟合值反应的是包括所有固定和随机效应的结果，lmer中通过fitted()函数获得该值。预测值，是可以设定不包括随机效应的，lmer中通过predict()函数获得该值。
shrimp.lm.9.predict <- predict(shrimp.lm.9,re.form=NA) #拟合值
shrimp.lm.9.predict.dt <- data.table(ObsSeq =as.integer(names(shrimp.lm.9.predict)),PredictedValue=shrimp.lm.9.predict)
shrimp[,":="(ObsSeq=seq(nrow(shrimp)))]
#把拟合值合并到shrimp数据集
shrimp.predicted.value <- merge(shrimp,shrimp.lm.9.predict.dt,by = c("ObsSeq"),all.y = TRUE)
ggplot(data=shrimp.predicted.value,aes(x=PopID,y=PredictedValue,fill=FamilyID)) +
geom_boxplot(outlier.size = 0)+
labs(x="群体",y="收获体重预测值")+
theme_gray(base_size = 20)+
theme(legend.position = "none")
4.5 随机效应的显著性检验
模型中加入随机效应后，如何检验其显著性？
一般是通过似然比率检验（Likelihood ratio test, LRT)来实现。lmerTest包中的rand()函数提供该功能。
首先我们来看一下似然比率检验的结果：
library(lmerTest)
ranova(shrimp.lm.9) # rand弃用，现在是ranova，表示随机因子的显著性检验
> ranova(shrimp.lm.9)
ANOVA-like table for random-effects: Single term deletions
Model:
M2BW ~ PopID + SexID + TankID + (1 | PopID:FamilyID) + SexID:M1BW
npar logLik AIC LRT Df Pr(>Chisq)
<none> 10 -11535 23090
(1 | PopID:FamilyID) 9 -12206 24430 1341.4 1 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
4.6 获得每个群体的性能
调用emmeans包中的函数，计算四个群体的估计边际均值(estimated marginal means)，或者说最小二乘均值(least-squares means)。根据边际均值，我们可以对群体的性能进行排序和比较。
关于emmeans包，请参考日志最小二乘均值的估计模型。尽管该日志介绍的是lsmeans包，但用法跟emmeans包都是一样的。而且根据作者介绍，在不久的将来，emmeans包要替代lsmeans包。
注意，安装emmeans还需要pbkrtest包，这个包没有自动安装，需要手动安装。
> library(emmeans)
> shrimp.lm9.rgl <- ref_grid(shrimp.lm.9)
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, set emm_options(pbkrtest.limit = 4241) or larger,
but be warned that this may result in large computation time and memory use.
Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, set emm_options(lmerTest.limit = 4241) or larger,
but be warned that this may result in large computation time and memory use.
> emmeans(shrimp.lm9.rgl,"PopID")
PopID emmean SE df asymp.LCL asymp.UCL
Pop1 33.8 0.485 Inf 32.9 34.8
Pop2 32.2 0.490 Inf 31.2 33.1
Pop3 30.1 0.466 Inf 29.2 31.0
Pop4 28.0 0.537 Inf 26.9 29.0
Results are averaged over the levels of: SexID, TankID
Degrees-of-freedom method: asymptotic
Confidence level used: 0.95
从上边结果中查找emmean列，可以看到Pop1群体的边际均值最大，这表明四个群体中该群体性能最好。
什么，你意犹未尽？
快点击阅读原文，查看栾老师的个人博客吧，上面很育种，很数据分析，很专业。
ps：微信不支持链接到netlify上，链接附上： https://luansheng.netlify.com/2017/11/19/lmm-tutorial-1/
https://luansheng.netlify.com/2017/11/29/lmm-tutorial-2/#8