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