混合线性模型不同模型拟合的可视化「建议收藏」

大家好，欢迎来到IT知识分享网。

1. 数据

这里使用sleepstudy数据集，看一下免费的R包lme4和付费包asreml如何处理不同的混合线性模型，以加深对混合线性模型的理解。

数据描述：

睡眠剥夺研究中受试者每天的平均反应时间。第0天，受试者有正常的睡眠时间。从那天晚上开始，他们每晚只能睡3个小时，依次进行0~9天。观察结果（y变量）代表了每天对每个受试者进行的一系列测试的平均反应时间。

数据预览：

> head(dat)
  Reaction Days Subject
1 249.5600    0     308
2 258.7047    1     308
3 250.8006    2     308
4 321.4398    3     308
5 356.8519    4     308
6 414.6901    5     308

IT知识分享网

• Reaction:为观测值，遇到刺激的反应时间
• Days：剥夺睡眠的天数
• Subject：实验对象（ID）

原数据可视化：

这里，Subject为每个实验对象（人），做一下折线图，看一下不同人在不同天数的反应时间。

IT知识分享网library(lme4)
data("sleepstudy")
dat = sleepstudy
ggplot(dat,aes(x = Days, y = Reaction, group = Subject, color = Subject)) +
  geom_point() + geom_line()

可以看到，不同的人差别比较大，不同的处理天数差别比较大，但是具体到每个人变化是不同的。

• 有些人在0天时，反应时间比如高（截距），后面随着天数的增加，增加得快，或者增加的慢（斜率）
• 有些人在0天时，反应时间比较短，后面随着天数的增加，增加得快，或者增加得慢
• 截距：intercept，为0天的反应
• 斜率：Slope，为增加的速度

lmer常用模型公式如下：

mod= lmer(data = , formula = y ~ Fixed_Factor + (Random_intercept + Random_Slope | Random_Factor))

• data，为数据集
• y，为观测值，所要分析的性状，因变量
• Fixed_Factor，为固定因子
• ()内为随机因子
  • Random_intercept，为随机截距，即认为不同群体因变量的分布不同（通俗的解释：有的人生下来起点高，是富二代，有的人是一般群众，起点低）
  • Random_Slope，为随机斜率，即认为不同群体受固定因子的影响不同（通俗解释：有的人是学霸，学习能力强，2个小时学会，斜率高；有的人是学渣，2天才能学会，斜率低）
  • Random_Factor，随机因子

参考: https://zhuanlan.zhihu.com/p/63092231

2. 模型1：截距随机，斜率固定

这种模型，认为不同人起点有差异，但是随着剥夺睡眠，他们的变化趋势没有差异（平行的）

这里，随机因子为(1 | Subject)，可以扩展为（1 + 1|Subject）认为斜率是固定的，截距是随机的。

IT知识分享网library(lme4)
mod1a = lmer(Reaction ~ Days + (1 | Subject), data=dat)
summary(mod1a)

library(asreml)
mod1b = asreml(Reaction ~ Days, random = ~ Subject, data=dat)

两者结果一致：

使用asreml软件的predict函数进行预测，查看预测值的分布：

pre = predict(mod1b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)

ggplot(pre,aes(x = Days, y = pre_value, 
               group = Subject, color = Subject)) +
  geom_point() + geom_line()

可以看出，每个人的斜率是一样的，截距不一样。

由结果可以写出拟合的函数，比如fixed的值：

• Days：10.46729，为斜率slope
• Intercept：251.4051，为整体均值，加上各个Subject的效应值，即为各个Subject的斜率

比如 Subject308，它的截距为：251.4051 + 40.7786 = 292.1837，那么它的公式为：

$$y=292.1837+10.46729\ast x$$

如框内所示：

3. 模型2：截距随机，斜率随机，没有相关性

这里，不考虑相关性的写法是||

mod2a =  lmer(Reaction ~ Days + (Days || Subject), data=dat)
summary(mod2a)

mod2b =  asreml(Reaction ~ Days, random = ~ Subject/Days, data=dat)
summary(mod2b)$varcomp
summary(mod2b,coef=T)$coef.fixed

两者结果一致：

使用asreml软件的predict函数进行预测，查看预测值的分布：

pre = predict(mod2b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)

ggplot(pre,aes(x = Days, y = pre_value, 
               group = Subject, color = Subject)) +
  geom_point() + geom_line()

可以看出，不同Subject个体，截距不一样，斜率也不一样。

4. 模型3：截距随机，斜率随机，考虑相关性

mod3a =  lmer(Reaction ~ Days + (Days | Subject), data=dat)
summary(mod3a) # 它没有估计出协方差，asreml可以估计出协方差

mod3b =  asreml(Reaction ~ Days, 
                random = ~ str(~Subject + Subject:Days, ~us(2):id(Subject)),
                data=dat)
summary(mod3b)$varcomp
summary(mod3b,coef=T)$coef.fixed
summary(mod3b,coef=T)$coef.random

结果一致，asreml结果更完整：

使用asreml软件的predict函数进行预测，查看预测值的分布：

pre = predict(mod3b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)

ggplot(pre,aes(x = Days, y = pre_value, 
               group = Subject, color = Subject)) +
  geom_point() + geom_line()

4. 模型3：斜率随机，截距固定

mod4a =  lmer(Reaction ~ Days + (0 + Days | Subject), data=dat)
summary(mod4a)

mod4b =  asreml(Reaction ~ Days, 
                random = ~ Subject:Days,
                data=dat)
summary(mod4b)$varcomp
summary(mod4b,coef=T)$coef.fixed

ranef(mod4a)
summary(mod4b,coef=T)$coef.random

结果一致：

使用asreml软件的predict函数进行预测，查看预测值的分布：

pre = predict(mod4b,classify = "Days:Subject",levels = list(Days = 0:9))
pre = as.data.frame(pre) %>% select(Days =1,Subject=2,pre_value =3)
head(pre)

ggplot(pre,aes(x = Days, y = pre_value, 
               group = Subject, color = Subject)) +
  geom_point() + geom_line()

可以看到，该模型截距都一样，截距固定。

5. lme4包模型比较

模型1：

mod1a = lmer(Reaction ~ Days + (1 | Subject), data=dat)

模型2：

mod2a =  lmer(Reaction ~ Days + (Days || Subject), data=dat)

模型3：

mod3a =  lmer(Reaction ~ Days + (Days | Subject), data=dat)

模型4：

mod4a =  lmer(Reaction ~ Days + (0 + Days | Subject), data=dat)

模型比较：

> anova(mod1a,mod2a,mod3a,mod4a)
refitting model(s) with ML (instead of REML)
Data: dat
Models:
mod1a: Reaction ~ Days + (1 | Subject)
mod4a: Reaction ~ Days + (0 + Days | Subject)
mod2a: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
mod3a: Reaction ~ Days + (Days | Subject)
      npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
mod1a    4 1802.1 1814.8 -897.04   1794.1                          
mod4a    4 1782.1 1794.8 -887.04   1774.1 19.9983  0               
mod2a    5 1762.0 1778.0 -876.00   1752.0 22.0771  1  2.619e-06 ***
mod3a    6 1763.9 1783.1 -875.97   1751.9  0.0639  1     0.8004    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看出，mod2a为最优模型，它与mod3a不显著，与其他模型达到极显著水平。它的AIC和BIC也最低，是最优模型。
即模型中截距随机，斜率随机的模型最优：

6. asreml包模型比较

模型1：

mod1b = asreml(Reaction ~ Days, random = ~ Subject, data=dat)

模型2：

mod2b =  asreml(Reaction ~ Days, random = ~ Subject/Days, data=dat)

模型3：

mod3b =  asreml(Reaction ~ Days, 
                random = ~ str(~Subject + Subject:Days, ~us(2):id(Subject)),
                data=dat)

模型4：

mod4b =  asreml(Reaction ~ Days, 
                random = ~ Subject:Days,
                data=dat)

模型比较：

aic1 = summary(mod1b)$bic
aic2 = summary(mod2b)$bic
aic3 = summary(mod3b)$bic
aic4 = summary(mod4b)$bic

bic1 = summary(mod1b)$bic
bic2 = summary(mod2b)$bic
bic3 = summary(mod3b)$bic
bic4 = summary(mod4b)$bic

re = data.frame(Model = paste0("Model: ",1:4), AIC = c(aic1,aic2,aic3,aic4),
                BIC = c(bic1,bic2,bic3,bic4))
re

结果：

> re
     Model      AIC      BIC
1 Model: 1 1469.687 1469.687
2 Model: 2 1432.073 1432.073
3 Model: 3 1437.213 1437.213
4 Model: 4 1449.746 1449.746

可以看出，模型2最优，它的AIC和BIC结果最低。

注意，asreml和lme4计算AIC和BIC的方法不一样，所以结果有所差异。

使用LRT似然比检验模型：

lrt.asreml(mod1b,mod2b)
lrt.asreml(mod1b,mod3b)
lrt.asreml(mod2b,mod3b)

结果：

> lrt.asreml(mod1b,mod3b)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

            df LR-statistic Pr(Chisq)    
mod3b/mod1b  2       42.837 1.545e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> lrt.asreml(mod2b,mod3b)
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

            df LR-statistic Pr(Chisq)
mod3b/mod2b  1     0.041056    0.4197

模型2最优，模型2和模型3不显著，模型2和模型1达到极显著。

7. 模型2的散点图和拟合图合并

这个模型最优！

欢迎关注我的公众号：育种数据分析之放飞自我。主要分享R语言，Python，育种数据分析，生物统计，数量遗传学，混合线性模型，GWAS和GS相关的知识。