学习线型混合效应模型最好的方法，是一边学习理论，一边动手实践。这样印象最为深刻。 本文参考了Linear models and linear mixed effects models in R的结构。 本文中，使用的数据集是来自文件shrimp.csv中的对虾育种数据。

1 读取数据文件

fread函数来自data.table包，特点是对大文件读取速度特别快。建议使用data.table作为数据处理的主力包。

其中的参数sep表示文件是用逗号分割，header表示数据文件中第一行是列名，stringsAsFactors表示读入的数据，对于字符类型，是否自动处理为因子类型。为了方便后边模型处理，这里设置为因子类型。

str()函数可以对数据集有一个汇总。从中看可以看出，AnimalID-个体编号，SireID-父本编号，DamID-母本编号，FamilyID-家系编号，SexID-性别，TankID-测试池号等字符类型，已经被设置为因子类型了。M1BW-入池前体重，M2BW-收获体重和M2Age-收获时日龄均为数字变量。

shrimp <- fread(input = "shrimp.csv",sep = ",",header = TRUE,stringsAsFactors = TRUE)
str(shrimp)
## Classes 'data.table' and 'data.frame':   4282 obs. of  9 variables:
##  $ AnimalID: Factor w/ 4282 levels "13G1000001","13G1000002",..: 1322 1310 1317 1684 2342 2358 3897 3979 3329 3504 ...
##  $ SireID  : Factor w/ 100 levels "12G000K010","12G000K065",..: 12 12 12 12 12 12 76 76 12 12 ...
##  $ DamID   : Factor w/ 91 levels "12G000K052","12G000K097",..: 70 70 70 70 70 70 10 10 70 70 ...
##  $ FamilyID: Factor w/ 105 levels "13F1306003","13F1306004",..: 21 21 21 21 21 21 61 61 21 21 ...
##  $ SexID   : Factor w/ 2 levels "Female","Male": 2 1 1 1 1 1 1 1 2 2 ...
##  $ TankID  : Factor w/ 2 levels "T1","T2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ M1BW    : num  3.5 3.5 3.5 3.5 3.5 3.5 3.58 3.58 3.76 3.76 ...
##  $ M2BW    : num  23 25.6 28.5 33.2 34.2 34.6 24.3 26.7 26.1 26.2 ...
##  $ M2Age   : int  212 212 212 212 212 212 222 222 213 213 ...
##  - attr(*, ".internal.selfref")=<externalptr>

数据文件准备完毕，我们正式开始学习线性混合效应模型。

2 引子：线性混合效应模型可以做什么

我们从一个简单的问题开始。对虾是有性别的，分雌雄。如果你对对虾没有任何了解，你可能会想知道，雄虾和雌虾的体重差别大吗？

我们测定了4282尾虾的体重。

我们首先直观的看一下雌雄虾体重分布点图。

ggplot(data=shrimp,aes(x=SexID,y=M2BW,color=SexID))+geom_boxplot()+geom_dotplot(binaxis = "y",stackdir = "center",position = "dodge",binwidth = 0.25)
图1 雌雄虾体重分布点图

图1 雌雄虾体重分布点图

从图1中可以大体看出，雌虾体重比雄虾高。然后，我们实际计算雌雄虾体重均值，发现雌虾的确比雄虾重。

shrimp.sex.m2bw <- shrimp[,mean(M2BW,na.rm=TRUE),by=.(SexID)]
shrimp.sex.m2bw
##     SexID       V1
## 1:   Male 27.77438
## 2: Female 33.93992

我们不禁要问，雌雄虾的这个差别是真实的吗？这是雌雄虾体重真正的差值吗？这个结果有没有偏差？

分析一下数据，你会发现，雌雄虾分布在2多个测试池中，且来自105个家系，并且每个家系的入池体重是存在差别的。而且你会发现，两个池子的养殖管理水平存在较大差别。

ggplot(shrimp,aes(x=TankID,y=M2BW,color=TankID))+geom_boxplot()+geom_dotplot(binaxis = "y",stackdir = "center",position = "dodge",binwidth = 0.3)

我们需要去除测试池、家系对雌雄虾体重的影响，准确估计雌雄虾体重。接下来就要用到线性模型了。

3 线性混合效应模型简介

\[M2BW = Sex + \varepsilon\] 表示一尾虾的体重由性别和随机误差决定。其中\[Sex\]作为固定效应，\[\varepsilon\]作为随机效应。后者表示所有影响体重的不可测量的效应总和，是随机和不可控制的。

从数据中我们发现，一尾虾的体重还受它所在的测试池和所在家系的影响。因此，这两个效应也需要放到模型中。模型进一步变为：

\[M2BW = Sex + Tank + Family +\varepsilon\] 新加入的两个变量，Tank和Family，如果都作为固定效应。那么上述模型称为线性模型。如果Family作为随机效应，那么上述模型称为线性混合效应模型（固定效应+随机效应）。

这里碰到的一个棘手问题是，模型中一个效应到底是作为固定效应，还是随机效应？准确的说，应该是与研究目的相关。

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)
summary(shrimp.lm)
## 
## Call:
## lm(formula = M2BW ~ SexID, data = shrimp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.5399  -2.8744   0.1256   3.0601  15.2601 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  33.9399     0.0974   348.4   <2e-16 ***
## SexIDMale    -6.1655     0.1404   -43.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.591 on 4280 degrees of freedom
## Multiple R-squared:  0.3105, Adjusted R-squared:  0.3104 
## F-statistic:  1928 on 1 and 4280 DF,  p-value: < 2.2e-16

现在解释上边的结果。首先来看Multiple R-squared: 0.3105，它表示模型对总体方差的解释能力。具体意思可以解释为，总体方差中的31%，可以由这个模型来解释。Adjusted R-squared是对Multiple R-squared的矫正，主要是考虑了固定效应。固定效应越多，该值越低。

下一个概念是非常的重要，那就是P值。