线性混合效应模型教程1
Sheng Luan
2017年11月19日
4.4 包括随机效应的线性混合效应模型
请加载一个新的数据集shrimpex.csv,其中有一个PopID字段，包括Pop1到Pop4共计4个水平，表示shrimp数据由四个群体组成。现在考虑这样一个问题：四个群体间收获体重是否存在差异。
首先加载数据文件。画出四个群体收获体重的箱形图，加上jitter点。
shrimp <- fread(input = "shrimpex.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>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")
如上图所示，如果我们分析时不考虑群体内的家系结构，那么家系方差会被累加到残差方差中，导致残差方差值被高估。
\[M2BW = Pop + Sex + Tank + Sex:M1BW\]
模型8不考虑家系结构，Pop、Sex和Tank为固定效应，Sex:M1BW为协变量。
\[M2BW = Pop + Sex + Tank + Sex:M1BW + Pop:Family\]
模型9考虑家系结构，Pop:Family为随机效应，Pop、Sex和Tank为固定效应，Sex:M1BW为协变量。
模型8的分析结果
shrimp.lm.8 <- lm(M2BW ~ 1 + PopID + SexID + TankID + SexID:M1BW,shrimp)
summary(shrimp.lm.8)##
## 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在模型中加入随机效应，需要使用lme4包中的lmer函数。
模型9的分析结果
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 ['lmerMod']
## 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 t value
## (Intercept) 35.7188 1.0605 33.68
## PopIDPop2 -1.6596 0.6884 -2.41
## PopIDPop3 -3.7742 0.6733 -5.61
## PopIDPop4 -5.8638 0.7359 -7.97
## SexIDMale -5.6599 0.5901 -9.59
## TankIDT2 -2.9491 0.1097 -26.88
## SexIDFemale:M1BW 0.3757 0.1204 3.12
## SexIDMale:M1BW 0.2904 0.1231 2.36
##
## 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把模型8的Residual standard error的平方，与模型9 Random Effects部分对比，你会发现，如果不考虑家系结构，残差方差明显被高估,估计值为18.4559804 。考虑家系结构后，残差方差为12.5267428, 明显变小， 从残差中分离出了大部分的家系方差。
4.5 获得每个群体的性能
根据教程1对于固定效应和随机效应的讨论，由于我们的目的是要分析四个群体间的差异，获得每个群体的性能，因此群体更适合做固定效应。每个群体是由多个家系组成的，这些家系只是大量家系的一个随机抽样，因此更加适合作为随机效应。
对于模型9，把Pop转换为固定效应，重新进行分析：
shrimp.lm.10 <- lmer(M2BW ~ 1 + PopID +SexID + TankID + SexID:M1BW + (1|PopID:FamilyID),shrimp)
summary(shrimp.lm.10)## Linear mixed model fit by REML ['lmerMod']
## 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 t value
## (Intercept) 35.7188 1.0605 33.68
## PopIDPop2 -1.6596 0.6884 -2.41
## PopIDPop3 -3.7742 0.6733 -5.61
## PopIDPop4 -5.8638 0.7359 -7.97
## SexIDMale -5.6599 0.5901 -9.59
## TankIDT2 -2.9491 0.1097 -26.88
## SexIDFemale:M1BW 0.3757 0.1204 3.12
## SexIDMale:M1BW 0.2904 0.1231 2.36
##
## 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然后，调用emmeans包中的函数，计算四个群体的估计边际均值(estimated marginal means)，或者说最小二乘均值(least-squares means)。根据边际均值，我们可以对群体的性能进行排序和比较。
关于emmeans包，请参考日志最小二乘均值的估计模型。尽管该日志介绍的是lsmeans，但用法跟emmeans包都是一样的。而且根据作者介绍，在将来，emmeans包要替代lsmeans包。
注意，安装emmeans还需要pbkrtest包，这个包没有自动安装，需要手动安装。
require(emmeans)## Loading required package: emmeansshrimp.lm10.rgl <- ref_grid(shrimp.lm.10)## Loading required namespace: pbkrtest## 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.emmeans(shrimp.lm10.rgl,"PopID")## PopID emmean SE df asymp.LCL asymp.UCL
## Pop1 33.83782 0.4853477 NA 32.88656 34.78909
## Pop2 32.17826 0.4903717 NA 31.21715 33.13937
## Pop3 30.06359 0.4659862 NA 29.15028 30.97691
## Pop4 27.97398 0.5373764 NA 26.92074 29.02722
##
## Results are averaged over the levels of: SexID, TankID
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95从上边结果中，可以看到Pop1群体的性能最好。