混合线性模型 | 常用模型与代码演示
昨天把win10升级到了win11,用着很流畅,推荐大家升级。
昨天群里面有老师问了一个问题,lme4包报错了:
看报错,应该是Rcpp版本过低导致的,我建议老师重新安装一下lme4和Rcpp,如果还不成功,那就回到lib目录,手动删除这两个包,然后再重新安装,毕竟之前写过经验贴:R包安装失败之粗暴解决方法,以及不用砸电脑成功安装R包的方法。
我的电脑lme4没有什么问题,看一下实例数据:
library(lme4)
data("sleepstudy")
dat = sleepstudy
mod1a = lmer(Reaction ~ Days + (1 | Subject), data=dat)
summary(mod1a)
关于混合线性模型,常用模型的拟合方法,之前写过一次总结,这里再放一遍,希望对后来者有所帮助。
这里使用sleepstudy
数据集,看一下免费的R包lme4
和付费包asreml
如何处理不同的混合线性模型,以加深对混合线性模型的理解。
共进行下面几种演示:
- 随机斜率,相同截距(Random slopes with same intercept)
- 随机斜率,随机截距,没有相关性(Random slopes and random intercepts (without correlation))
- 随机斜率,随机截距,有相关性(Random slopes and random intercepts (with correlation)
- 随机斜率,不同截距(Random slopes with a different intercept)
- 其它lme4不能实现的功能
0. 数据描述
❝睡眠剥夺研究中受试者每天的平均反应时间。第0天,受试者有正常的睡眠时间。从那天晚上开始,他们每晚只能睡3个小时。观察结果代表了每天对每个受试者进行的一系列测试的平均反应时间。❞
> 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)
> table(dat$Days)
0 1 2 3 4 5 6 7 8 9
18 18 18 18 18 18 18 18 18 18
> table(dat$Subject)
308 309 310 330 331 332 333 334 335 337 349 350 351 352 369 370 371 372
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
共有18个人,处理天数是0~9天。
这里,Subject为因子,Days为数字类型,Reaction为数字类型。
1. 随机斜率,相同截距
通用的混线性模型,包括固定因子和随机因子。育种中常用的混线性模型,一般是针对于随机因子关系矩阵,而其它领域的一般是针对于不同截距的定义。
- 固定因子:Days
- 随机因子:Subject
「lme4:」
这里,Reaction为y变量,Days为固定因子,随机因子用(1 | Subject)表示。
lme4的基本语法:
library(lme4)
mod1a = lmer(Reaction ~ Days + (1 | Subject), data=dat)
summary(mod1a)
结果:
> summary(mod1a)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
Data: dat
REML criterion at convergence: 1786.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.2257 -0.5529 0.0109 0.5188 4.2506
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1378.2 37.12
Residual 960.5 30.99
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.4051 9.7467 25.79
Days 10.4673 0.8042 13.02
Correlation of Fixed Effects:
(Intr)
Days -0.371
这里,
- 随机因子的方差组分为1378.2,残差的方差组分为960.5。
- 固定因子中,Days为数值协变量,截距的值为251.4,Days的值为10.46
「asreml:」
代码:
library(asreml)
mod1b = asreml(Reaction ~ Days, random = ~ Subject, data=dat)
summary(mod1b)$varcomp
summary(mod1b,coef=T)$coef.fixed
结果:
> mod1b = asreml(Reaction ~ Days, random = ~ Subject, data=dat)
Model fitted using the gamma parameterization.
ASReml 4.1.0 Tue May 11 19:06:29 2021
LogLik Sigma2 DF wall cpu
1 -756.229 1572.711 178 19:06:29 0.0
2 -746.076 1354.228 178 19:06:29 0.0
3 -736.843 1164.250 178 19:06:29 0.0
4 -731.820 1050.393 178 19:06:29 0.0
5 -729.920 986.583 178 19:06:29 0.0
6 -729.670 964.724 178 19:06:29 0.0
7 -729.661 960.621 178 19:06:29 0.0
> summary(mod1b)$varcomp
component std.error z.ratio bound %ch
Subject 1378.4099 504.5367 2.732031 P 0.2
units!R 960.6208 107.0758 8.971412 P 0.0
> summary(mod1b,coef=T)$coef.fixed
solution std error z.ratio
Days 10.46729 0.8042902 13.01432
(Intercept) 251.40510 9.7400376 25.81151
这里,
- 可以看到Subject的方差组分为1378.4,残差的方差组分为960.6,
- Days的值为10.46,截距的值为251.405
计算随机因子的BLUP值:
两者结果一致!
2. 随机斜率,随机截距,没有相关性
这里模型更复杂一点,假定不同的人(项目)有各自的截距,并且他们之间不相关。
❝This is the a model where you assume that the random effect has different intercepts based on the levels of another variable. In addition the || in lme4 assumes that slopes and intercepts have no correlation.❞
「lme4」
mod2a = lmer(Reaction ~ Days + (Days || Subject), data=dat)
summary(mod2a)
结果:
> summary(mod2a)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + ((1 | Subject) + (0 + Days | Subject))
Data: dat
REML criterion at convergence: 1743.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.9626 -0.4625 0.0204 0.4653 5.1860
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 627.57 25.051
Subject.1 Days 35.86 5.988
Residual 653.58 25.565
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 251.405 6.885 36.513
Days 10.467 1.560 6.712
Correlation of Fixed Effects:
(Intr)
Days -0.184
这里,随机因子有了交互效应。
「asreml:」这里,Days是数字变量,Subject/Days,相当于是Subject + Subject:Days,即Subject为随机因子,每个Subject都有一个Days的系数。
mod2b = asreml(Reaction ~ Days, random = ~ Subject/Days, data=dat)
summary(mod2b)$varcomp
结果:
> summary(mod2b)$varcomp
component std.error z.ratio bound %ch
Subject 627.67842 282.67774 2.220474 P 0.3
Subject:Days 35.86583 14.54252 2.466273 P 0.1
units!R 653.70536 76.67803 8.525328 P 0.0
方差组分结果一致:
看一下两者随机因子的效应值:
结果一致。
3. 随机斜率,随机截距,有相关性
❝This is the a model where you assume that the random effect has different intercepts based on the levels of another variable. In addition a single | in lme4 assumes that slopes and intercepts have a correlation to be estimated ❞
「lme4:」
mod3a = lmer(Reaction ~ Days + (Days | Subject), data=dat)
summary(mod3a)
「asreml:」
这里,asreml写法比较复杂,用了str
函数进行定义。
mod3b = asreml(Reaction ~ Days,
random = ~ str(~Subject + Subject:Days, ~us(2):id(Subject)),
data=dat)
summary(mod3b)$varcomp
结果:
结果一致。
看一下随机因子的效应值:
可以看到,结果一致!
4. 随机斜率,不同截距(Random slopes with a different intercept)
模型描述:
❝This is the a model where you assume that the random effect has different intercepts based on the levels of another variable but there’s not a main effect. The 0 in the intercept in lme4 assumes that random slopes interact with an intercept but without a main effect. ❞
「lme4:」
mod4a = lmer(Reaction ~ Days + (0 + Days | Subject), data=dat)
summary(mod4a)
「asreml:」这里,只考虑Subject:Days交互作用,不考虑Subject作为随机因子了。
mod4b = asreml(Reaction ~ Days,
random = ~ Subject:Days,
data=dat)
summary(mod4b)$varcomp
summary(mod4b,coef=T)$coef.fixed
结果比较:
结果一致。
看一下随机因子的效应值:
结果完全一致。
5. asreml能做但是lme4不能做的模型
- 比如diag模型
- 比如us模型
- 比如FA模型
- 比如leg模型
- 比如corgh模型
- ……
相关文章
- 【技术种草】cdn+轻量服务器+hugo=让博客“云原生”一下
- CLB运维&运营最佳实践 ---访问日志大洞察
- vnc方式登陆服务器
- 轻松学排序算法:眼睛直观感受几种常用排序算法
- 十二个经典的大数据项目
- 为什么使用 CDN 内容分发网络?
- 大数据——大数据默认端口号列表
- Weld 1.1.5.Final,JSR-299 的框架
- JavaFX 2012:彻底开源
- 提升as3程序性能的十大要点
- 通过凸面几何学进行独立于边际的在线多类学习
- 利用行动影响的规律性和部分已知的模型进行离线强化学习
- ModelLight:基于模型的交通信号控制的元强化学习
- 浅谈Visual Source Safe项目分支
- 基于先验知识的递归卡尔曼滤波的代理人联合状态和输入估计
- 结合网络结构和非线性恢复来提高声誉评估的性能
- 最佳实践丨云开发CloudBase多环境管理实践
- TimeVAE:用于生成多变量时间序列的变异自动编码器
- 具有线性阈值激活的神经网络:结构和算法
- 内网渗透之横向移动 -- 从域外向域内进行密码喷洒攻击