LMM线性混合模型by nlme
笔者邀请您,先思考:
1 如何理解和应用线性混合模型?
介绍
线性混合模型是数据模型中一种重要的类别,它可以分析存在相关关系的数据,模型里面包含固定效应以及随机效应,随机效应描述的是在不同层次的不同水平中,各变量对总体观察变量的贡献。

数据导入
本文选择的数据来自Andrzej Gałecki • Tomasz Burzykowski的nlme的armd数据。
1library(nlme)
2library(nlmeU)
3library(lattice)
4data(armd, package = "nlmeU")
观察数据发现时间作为factor有四个水平,分别是4、12、24、52周。
head(armd)
1 subject treat.f visual0 miss.pat time.f time visual tp
22 1 Active 59 --XX 4wks4 55 1
33 1 Active 59 --XX 12wks 12 45 2
45 2 Active 65 ---- 4wks4 70 1
56 2 Active 65 ---- 12wks 12 65 2
67 2 Active 65 ---- 24wks 24 65 3
78 2 Active 65 ---- 52wks 52 55 4
首先,我们使用lme拟合模型,用lm.form表示固定效应,包含了四部分,分别是测试初值visual0,时间time,治疗方式treat.f,以及治疗方式与时间的相互作用。
1lm.form <-
2 formula(visual ~ visual0 + time + treat.f + treat.f:time)
在构建混合线性模型时,我们需要加入一个随机效应的指定值。~1|subject表示subject各水平分别对应一个随机截距,也就是每一个体对应不同的随机截距。由于个体间存在差异,这样的设定可以有利于构建良好模型。lme使用的拟合方法默认为REML估计(限制极大似然估计)如果需要使用极大似然估计需要在method参数下设置”ML”。结果给出了拟合模型的AIC(赤池信息量准则),BIC(贝叶斯准则)以及logLIk(极大似然估计值)。这些信息可以用于比较不同模型的拟合效果,通常越接近零效果越好。
1(lmm.1 <-
2lme(lm.form,
3 random = ~1|subject, data = armd))
4
5Linear mixed-effects model fit by REML
6 Data: armd
7 AIC BIC logLik
8 6592 6625.3 -3289
9
10Random effects:
11 Formula: ~1 | subject
12(Intercept) Residual
13StdDev: 8.9782 8.6275
14
15Fixed effects: list(lm.form)
16 Value Std.Error DF t-value p-value
17(Intercept) 9.2881 2.68189 631 3.4633 0.0006
18visual0 0.8264 0.04467 231 18.5022 0.0000
19time -0.2122 0.02293 631 -9.2552 0.0000
20treat.fActive -2.4220 1.49997 231 -1.6147 0.1077
21time:treat.fActive -0.0496 0.03356 631 -1.4776 0.1400
22 Correlation:
23 (Intr) visul0 time trt.fA
24visual0-0.920
25time -0.185 -0.003
26treat.fActive -0.295 0.022 0.335
27time:treat.fActive 0.126 0.002 -0.683 -0.476
28
29Standardized Within-Group Residuals:
30 MinQ1 MedQ3 Max
31-4.187505 -0.396925 0.032048 0.551383 2.951321
32
33Number of Observations: 867
34Number of Groups: 234
如果使用常规的线性模型lm,把个体的差异忽略,我们可以得到另一组数据
1(lm.1 <- lm(visual ~ visual0 + time + treat.f + treat.f:time,data = armd))
2
3Call:
4lm(formula = visual ~ visual0 + time + treat.f + treat.f:time,
5data = armd)
6
7Coefficients:
8 (Intercept) visual0time
99.1004 0.8304 -0.2104
10 treat.fActive time:treat.fActive
11 -2.6547 -0.0367
通过计算两个模型的极大似然比检验统计量,对它们进行比较,利用极大似然比检验统计量的差异检验模型中的某些变量是否显著,该检验中,两个模型的差异部分是个体的随机效应。检验统计量属于卡方分布,零假设认为两个模型没有差异,即相差的变量不显著,lower.tail = F 对应的是P[X > x],结果显示拒绝零假设,认为个体的随机效应需要考虑。
1 pchisq(-2*(m1-mm1),df = 1,lower.tail = F)
2'log Lik.' 2.4793e-53 (df=6)
3> (mm1<- logLik(lmm.1))
4'log Lik.' -3289 (df=7)
5> (m1 <- logLik(lm.1))
6'log Lik.' -3407.2 (df=6)
7> pchisq(-2*(m1-mm1),df = 1,lower.tail = F)
8'log Lik.' 2.4793e-53 (df=6)
由于每个病人在不同时间呈现的效果有一定差异,所以在随机效应中有必要把各个体的时间因素单独考虑,在模型lmm.1的基础上,我们使用update更新模型,如果我们只需要在原来模型上稍作修改,update是一个不错的选择。
1lmm.2 <-
2update(lmm.1,
3random = ~1+ time|subject, data = armd)
为了更好的呈现不同个体水平对应的时间随机效应,我们使用qqnorm绘图。

1qqnorm(lmm.2, ~resid(.)|time.f)
针对时间的四个水平,可以分别获得对应的残差QQ图。52wks的qq线跟其它的有点不同,主要源于它的数据量相对其它的少了。
最后看一下lmm.2与lmm.1的差别,结果显示两个模型有显著差异,模型lmm.2更佳。
1anova(lmm.1,lmm.2)
2Model dfAICBIC logLik Test L.Ratio p-value
3lmm.1 1 7 6592.0 6625.3 -3289.0
4lmm.2 2 9 6453.8 6496.7 -3217.9 1 vs 2 142.15 <.0001
内容推荐
数据人网:数据人学习,交流和分享的平台,诚邀您创造和分享数据知识,共建和共享数据智库。
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!