育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

混合线性模型的可视化

已有 813 次阅读 2021-6-5 17:29 |个人分类:农学统计|系统分类:科研笔记

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
  • Reaction:为观测值,遇到刺激的反应时间
  • Days:剥夺睡眠的天数
  • Subject:实验对象(ID)

「原数据可视化:」

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

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)认为斜率是固定的,截距是随机的。

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,那么它的公式为:

如框内所示:

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相关的知识。




http://blog.sciencenet.cn/blog-2577109-1289879.html

上一篇:混合线性模型的理解 | 学霸和富二代的故事
下一篇:R语言循环第三境界:purrr包map函数!

1 李宏翰

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-9-25 02:03

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部