R语言之纵向数据分析：多级线性增长模型2

这篇文章已经是纵向数据分析系列的第三篇了。之前，我们介绍了什么是纵向数据，我们如何把长型的数据集转换成纵向的数据集，并通过建立相应的多级模型进行分析。显然，仅仅介绍基本的多级模型并不足以对虚拟随机对照实验数据集进行分析。这篇文章则延续之前的分析，并介绍一种不错的方法来处理多级模型里的治疗效果问题。

创建一个纵向数据集并把它转成长型格式

和往常一样，我们还是先创建一个纵向数据集，然后把它转成长型的格式。

1. library(MASS)
3. dat.tx.a<-mvrnorm(n=250,mu=c(30,20,28),
4. Sigma=matrix(c(25.0,17.5,12.3,
5. 17.5,25.0,17.5,
6. 12.3,17.5,25.0),nrow=3,byrow=TRUE))
8. dat.tx.b<-mvrnorm(n=250,mu=c(30,20,22),
9. Sigma=matrix(c(25.0,17.5,12.3,
10. 17.5,25.0,17.5,
11. 12.3,17.5,25.0),nrow=3,byrow=TRUE))
13. dat<-data.frame(rbind(dat.tx.a,dat.tx.b))
14. names(dat)<-c('measure.1','measure.2','measure.3')
16. dat<-data.frame(subject.id=factor(1:500),tx=rep(c('A','B'),each=250),dat)
18. rm(dat.tx.a,dat.tx.b)
20. dat<-reshape(dat,varying=c('measure.1','measure.2','measure.3'),
21. idvar='subject.id',direction='long')

多级增长模型分析治疗效果

这是一个RCT数据集，也就说明，这里潜在存两个以疗效分组的差异。上一次，我们没有考虑异质性，而仅仅对两个组的普通疗效进行了细分。从直觉上看，我们可以添加治疗变量来修复这个模型，从而抓到组间的差异。

1. library(lmerTest)
2. m1<-lmer(measure~time+tx+(1|subject.id),data=dat)

这里，我们任然使用lmerTest包是因为它允许修正模型的检验使用自由度。模型的公式和是和上次的一样，除非我们把治疗变量tx作为独立变量。

1. summary(m1)
2. Linearmixed model fitbyREML t-testsuseSatterthwaiteapproximations to degrees of freedom[
3. lmerMod]
4. Formula:measure~time+tx+(1|subject.id)
5. Data:dat
7. REML criterion at convergence:9700.8
9. Scaledresiduals:
10. Min1QMedian3QMax
11. -3.03326-0.656180.080810.658932.64555
13. Randomeffects:
14. GroupsNameVarianceStd.Dev.
15. subject.id(Intercept)6.8342.614
16. Residual31.9405.652
17. Numberof obs:1500,groups:subject.id,500
19. Fixedeffects:
20. EstimateStd.Errordf t valuePr(>|t|)
21. (Intercept)30.61470.44461494.000068.856<2e-16***
22. time-2.34340.1787999.0000-13.112<2e-16***
23. txB-2.29110.3740498.0000-6.1271.82e-09***
24. ---
25. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
27. CorrelationofFixedEffects:
28. (Intr)time
29. time-0.804
30. txB-0.4210.000

从上述的summary里的信息可以很明显的看到，time和tx这两个变量都相当显著，而接受B治疗的人群在抑郁症得分上比接受A治疗的人群低2.34分。那么，我们是否可以认为，B治疗方案效果更好？绝对不是。模型m1实际上并不合适对纵向数据集里的组间差异进行比较。m1里的治疗效果，其实也就是一段时间内的平均治疗效果。换句话说，治疗方案A和治疗方案B的疗效所导致的分差在-2.34分，其实只不过是1,2,3时间段内抑郁症均分的差。这是RCT，所以，我们所期待的志愿者参与两组治疗方案里时间段1的差要小于1。我们对两组疗效的均值差异感兴趣，而现在的轨迹在时间段2、3里存在差异。
这个概率可能对初学者来说有点复杂，但是其执行过程很简单，我们只需要添加交互变量time和tx。

1. m2<-lmer(measure~time*tx+(1|subject.id),data=dat)

在R，我们可以使用:符号来表明变量间的交互关系。因此，我们可以这样设置：time+tx+time:tx。但是，更简短的写法是：timetx。操作符要R涵盖主要变量和治疗疗效。此时，使用操作符就能简单实现各个效应的展示。请记住，操作符也可用于高维变量交互。比如，ABC就要R包含三个主要的效应，三个两两交互的效应，以及一个三元交互效应。

1. summary(m2)
2. Linearmixed model fitbyREML t-testsuseSatterthwaiteapproximations to degrees of freedom[
3. lmerMod]
4. Formula:measure~time*tx+(1|subject.id)
5. Data:dat
7. REML criterion at convergence:9640.8
9. Scaledresiduals:
10. Min1QMedian3QMax
11. -3.0454-0.67760.10410.66222.4296
13. Randomeffects:
14. GroupsNameVarianceStd.Dev.
15. subject.id(Intercept)7.4482.729
16. Residual30.1005.486
17. Numberof obs:1500,groups:subject.id,500
19. Fixedeffects:
20. EstimateStd.Errordf t valuePr(>|t|)
21. (Intercept)27.88080.55741421.500050.017<2e-16***
22. time-0.97650.2454998.0000-3.9807.40e-05***
23. txB3.17670.78831421.50004.0305.88e-05***
24. time:txB-2.73390.3470998.0000-7.8798.44e-15***
25. ---
26. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
28. CorrelationofFixedEffects:
29. (Intr)time txB
30. time-0.880
31. txB-0.7070.622
32. time:txB0.622-0.707-0.880

现在，我们从m2模型中获得一个三元交互关系。在这个模型当中，time效应指的就是疗效A（这是特定条件下的疗效和疗效时间的交互），而txB效应就是时间点为0时的效应（数据集里一个虚的时间点）。这个术语解释起来有点复杂：系数值反映了每个时间段内两种治疗方案的疗效上的差异，其条件是取平均时间和治疗效果。口头上说，人们可以理解的是，这表明两组治疗方案是如何随着时间的推移而不同。对于每个单位时间内的增长（即，从时间段1到时间段2），在完成平均治疗效应差的分析以后，接受B治疗方案的参与者理论上在抑郁症得分方面比接受A治疗方案的要低2.88分。所以，在时间段2内，接受B治疗方案的参与者得分差是 -2.88 x 2 + 3.43 = -2.33，而在时间段3内，接受B治疗方案的参与者得分差是-2.88 x 3+ 3.43 = -5.21。我们应当时刻注意，在一个给定的时间段里对比疗效差的时候，主要疗效以及疗效变量间的交互。但是，为什么要这样？你不妨设置分析一下，当你忽略了主要治疗效应的时候，你在时间段1内所得到的疗效差是怎样的情形。

选择最后的模型

我可以告诉你，模型m2比m1和m0（上一篇文章所采用的模型）都要好，原因在于，这个数据集是我们自己创建的。但在实际生活中，我们会在选择最好的一个模型的时候会遇到很多困难。其中一种选择的方式就是，按照我们之前所讲的，使用summary()函数，以此算出额外变量的统计显著性。其它一个常用的方法就是使用方差分析（ANOVA)表（仅仅针对嵌套模型），而这种方法适合用于这个模型。

1. m0<-lmer(measure~time+(1|subject.id),data=dat)
2. anova(m0,m1,m2)
4. Data:dat
5. Models:
6. object:measure~time+(1|subject.id)
7. ..1:measure~time+tx+(1|subject.id)
8. ..2:measure~time*tx+(1|subject.id)
9. DfAIC BIC logLik devianceChisqChiDfPr(>Chisq)
10. object49741.99763.1-4866.99733.9
11. ..159707.69734.1-4848.89697.636.33411.663e-09***
12. ..269649.29681.1-4818.69637.260.34517.962e-15***
13. ---
14. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

这里，我们使用anova()对m0、m1和m2模型进行方差分析，原因在于我们把m0嵌套于m1，然后把m1嵌套于m2。如果还有另外的子集，我们就把模型嵌套到那个子集上。换句话说，含有独立变量A、B的模型则嵌套于含有变量A、B、C的模型中。
如果你做过方差分析，那这个方差分析表对于你来说应该是很熟悉了。Df列表明了模型里的自由度，它简单的反映了这个案例的参数估计。例如，自由度为4的模型m0表明了这个模型预测了4个参数（1.修正截距效应系数，2.时间效应，3.随机截距方差，4.残差）。
AIC和BIC是两个常用的拟合指数。AIC和BIC都考虑到两个因素：模型的数据拟合效果有多好，和这个模型复不复杂。AIC和BIC的不同之处在于模型复杂度的罚分上（BIC对模型复杂度的罚分影响更大）。我发现，AIC理解起来更容易，因为它渐进的使用留一交叉检验法（LOOCV)来预测运行效果。我们应当理解交叉检验和预测性能的概念，随后，如果你对此不太理解，你可以先把它跳过。

LogLik列模型参数的似然对数。从根本上说，它是模型里数据拟合效果的指标。其偏差是−2×logLik−2×logLik。为什么我们需要它？因为两个嵌套模型的偏差在原假设的条件下遵循卡方检验（假设两个嵌套模型的偏差一样）。此时，我们可以把它用作模型拟合差的测试。偏差的差则在Chisq列显示。卡方检验的统计相关自由度则在Chi Df（简单来说，就是额外参数的个数）。最后一行的结果显而易见：它算出了测试两个嵌套模型的差的p值。虽然是这样，我们还是要警惕ANOVA的结果，因为在众多的候选模型中，它会导致错误类型1上升极快。

由于模型m2比其它两个候选模型AIC和BIC都要好，而对于方差分析，我们可以认为这个模型是目前为止我们拟合的最好的了。当然，这个模型还有许多要改进的地方。后续，我们可能还可以在模型的性能上进行提升，并作图进行预测。