R語(yǔ)言之縱向數(shù)據(jù)分析：多級(jí)增長(zhǎng)模

上一次，我們討論了如何對(duì)長(zhǎng)型數(shù)據(jù)轉(zhuǎn)換成長(zhǎng)型的數(shù)據(jù)，同時(shí)還是用了一個(gè)隨機(jī)創(chuàng)建的對(duì)照實(shí)驗(yàn)數(shù)據(jù)集來(lái)對(duì)其增長(zhǎng)趨勢(shì)進(jìn)行可視化。但是，我們是否能夠進(jìn)一步的分析并預(yù)測(cè)結(jié)果的增長(zhǎng)趨勢(shì)與時(shí)間之間的關(guān)系。
是的，當(dāng)然可以！我們可以使用多級(jí)增長(zhǎng)模型（也稱之為層次模型或者混合模型）進(jìn)行估計(jì)。

產(chǎn)生一個(gè)水平數(shù)據(jù)集并把它轉(zhuǎn)成寬格式
下面，我們先從我之前的一篇文章的實(shí)例進(jìn)行講解：

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')

多級(jí)增長(zhǎng)模型

這里有很多R語(yǔ)言包可以幫助你進(jìn)行多級(jí)分析，其中，我發(fā)現(xiàn)lme4包是最好的一個(gè)，因?yàn)樗褂帽容^簡(jiǎn)單，而且建模能力也很強(qiáng)（尤其是輸出二進(jìn)制結(jié)果或者計(jì)數(shù)結(jié)果）。當(dāng)然，nlme包也是相當(dāng)不錯(cuò)的，它可以給連續(xù)型結(jié)果提供了類似的結(jié)果（正態(tài)/高斯分布）。

1. install.packages('lme4')
2. library(lme4)
3. m<-lmer(measure~time+(1|subject.id),data=dat)

如果你之前做過(guò)回歸分析，你應(yīng)該對(duì)這樣的語(yǔ)法結(jié)構(gòu)比較熟悉了。通常來(lái)說(shuō)，它就是lm()函數(shù)當(dāng)中含有額外的隨即效應(yīng)公式。
隨即效應(yīng)，如果你對(duì)這個(gè)術(shù)語(yǔ)不熟悉的話，其實(shí)可以這么理解，通常來(lái)說(shuō)，它就是一個(gè)實(shí)驗(yàn)所無(wú)法控制的誤差，即變化。因此，比方來(lái)說(shuō)，一個(gè)志愿者所收到的治療效果就是一種混合的效應(yīng)，因?yàn)?，假設(shè)我們是實(shí)驗(yàn)人員，我們會(huì)決定哪些人接受A治療方案，哪些接受B治療方案。然而，抑郁癥評(píng)分的基線在治療的初始階段會(huì)因人而異，一些人可能會(huì)更加抑郁，一些其實(shí)并沒(méi)有這么憂郁。由于這是無(wú)法控制的，我們會(huì)把它看成是隨即效應(yīng)。
尤其是，抑郁評(píng)分基線的差異可以看作是一個(gè)隨機(jī)區(qū)間（即，不同的志愿者參與不同等級(jí)的治療）。我們也可以在建模的時(shí)候，對(duì)它們的斜率進(jìn)行隨機(jī)設(shè)置：例如，如果我們有理由相信盡管大家接受的治療是一樣的，一些參與治療的人可以收到很好的療效，而其它人則收效甚微。
結(jié)果的隨機(jī)效應(yīng)部分陳述了數(shù)據(jù)的方差結(jié)構(gòu)。在這個(gè)模型中，存在兩種方差結(jié)構(gòu)：殘差（通常用在線性模型）和個(gè)體之間的差異（即，每一個(gè)主體的id）。量化個(gè)體差異程度的一種常用方法就是研究同類相關(guān)系數(shù)（ICC）。我們可能可以從多級(jí)模型那里計(jì)算ICC，而且,這意味著，24.3%的抑郁平分變化可以由個(gè)體差異程度來(lái)解釋。
現(xiàn)在，我們把目光轉(zhuǎn)到修正效應(yīng)。嗯…，那些p值在哪里呢？這，盡管SAS和其它統(tǒng)計(jì)軟件有給多級(jí)模型的修正效應(yīng)計(jì)算提供p值方面的信息，其實(shí)，很多統(tǒng)計(jì)學(xué)家的計(jì)算結(jié)果并不一致。舉個(gè)簡(jiǎn)單的例子，我們對(duì)自由度與這些t檢驗(yàn)的關(guān)聯(lián)程度了解的不深，而且沒(méi)有自由度的話，我們比不知道t檢驗(yàn)的具體分布，因此，我們無(wú)法得到p值方面的信息。SAS和其它軟件都有相應(yīng)的工作區(qū)來(lái)處理估計(jì)值，這時(shí)lme4包開(kāi)發(fā)人員感到不舒服的地方。結(jié)果，lmer包并沒(méi)有刻意的匯報(bào)p值的信息（所以，不要害怕你得不到p值！或許有其它的方法在顯著性的測(cè)量上比我們的模型做的還好）。
這么說(shuō)，如果你絕對(duì)需要p值，我們可以使用基于lme4包所產(chǎn)生的lmerTest包來(lái)估算p值。

含有p估計(jì)值的多級(jí)增長(zhǎng)模型

下面大部分的代碼和上面的類似，除非我們要使用lmerTest包。

1. # install.packages('lme4')
2. # Please note the explanation and limitations:
3. # https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html
4. library(lmerTest)
5. m<-lmer(measure~time+(1|subject.id),data=dat)

其結(jié)果很相似，但現(xiàn)在，我們可以得到自由度和p的估計(jì)值。所以，我們可以很自信的說(shuō)普通RCT參與治療的人，現(xiàn)在，隨著時(shí)間的推移，他們的抑郁癥得分在下降，其速度為每下降1分，下降的量為2.24。

1. summary(m)
2. Linearmixed model fitbyREML t-tests
3. useSatterthwaiteapproximations to
4. degrees of freedom[merModLmerTest]
5. Formula:measure~time+(1|subject.id)
6. Data:dat
8. REML criterion at convergence:9639.6
10. Scaledresiduals:
11. Min1QMedian3QMax
12. -2.45027-0.705960.008320.659512.78759
14. Randomeffects:
15. GroupsNameVarianceStd.Dev.
16. subject.id(Intercept)9.2893.048
17. Residual28.8605.372
18. Numberof obs:1500,groups:subject.id,500
20. Fixedeffects:
21. EstimateStd.Errordf t valuePr(>|t|)
22. (Intercept)29.85080.39151449.400076.25<2e-16***
23. time-2.24200.1699999.0000-13.20<2e-16***
24. ---
25. Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
27. CorrelationofFixedEffects:
28. (Intr)
29. time-0.868

計(jì)算95%置信區(qū)間和預(yù)測(cè)區(qū)間

有時(shí)，我們想在單個(gè)軌跡的均值進(jìn)行作圖。如果要展示均值里的一些不確定因素，我們需要使用擬合好的模型，利用擬合值進(jìn)行計(jì)算，算出95%置信區(qū)間和95%預(yù)測(cè)區(qū)間。

1. # See for details: http://glmm.wikidot.com/faq
2. dat.new<-data.frame(time=1:3)
3. dat.new$measure<-predict(m,dat.new,re.form=NA)
5. m.mat<-model.matrix(terms(m),dat.new)
7. dat.new$var<-diag(m.mat%*%vcov(m)%*%t(m.mat))+VarCorr(m)$subject.id[1]
8. dat.new$pvar<-dat.new$var+sigma(m)^2
9. dat.new$ci.lb<-with(dat.new,measure-1.96*sqrt(var))
10. dat.new$ci.ub<-with(dat.new,measure+1.96*sqrt(var))
11. dat.new$pi.lb<-with(dat.new,measure-1.96*sqrt(pvar))
12. dat.new$pi.ub<-with(dat.new,measure+1.96*sqrt(pvar))

第一行代碼指出我們想要求出均值的一個(gè)點(diǎn)，它們一般來(lái)說(shuō)是在我們這個(gè)案例的前三次預(yù)測(cè)的時(shí)候。第二行代碼使用了predict()函數(shù)來(lái)得到模型的均值，它不考慮條件隨機(jī)效應(yīng)（re.form=NA)。第三第四行計(jì)算了均值的方差，一般來(lái)說(shuō)是矩陣交叉與隨機(jī)效應(yīng)截距相加。第五行計(jì)算了單個(gè)觀測(cè)值的方差，它的方差等于方差均值假設(shè)殘差方差。第六到第九行則按普通方法，并假設(shè)它是正態(tài)分布來(lái)計(jì)算95%置信區(qū)間和預(yù)測(cè)區(qū)間。最后所給的代碼是：

1. dat.new
2. time measurevarpvar ci.lb ci.ub pi.lb pi.ub
3. 127.7242110.8566943.0405421.2661134.1823114.86557440.58285
4. 225.2234210.8245143.0083518.7749031.6719412.36959238.07725
5. 322.7226310.8566943.0405416.2645329.180739.86399335.58127

作均值圖像

最后，我們要作它的95%置信區(qū)間和95%預(yù)測(cè)區(qū)間的圖像了。注意，預(yù)測(cè)區(qū)間的圖像要寬于置信區(qū)間。也就是說(shuō)，預(yù)測(cè)均值的結(jié)果比用單個(gè)值預(yù)測(cè)要好。

1. ggplot(data=dat.new,aes(x=time,y=measure))+
2. geom_line(data=dat,alpha=.02,aes(group=subject.id))+
3. geom_errorbar(width=.02,colour='red',
4. aes(x=time-.02,ymax=ci.ub,ymin=ci.lb))+
5. geom_line(colour='red',linetype='dashed',aes(x=time-.02))+
6. geom_point(size=3.5,colour='red',fill='white',aes(x=time-.02))+
7. geom_errorbar(width=.02,colour='blue',
8. aes(x=time+.02,ymax=pi.ub,ymin=pi.lb))+
9. geom_line(colour='blue',linetype='dashed',aes(x=time+.02))+
10. geom_point(size=3.5,colour='blue',fill='white',aes(x=time+.02))+
11. theme_bw()

如果你和我一樣，對(duì)數(shù)據(jù)也很敏感，你應(yīng)該能觀察到圖線的擬合效果并不太好。這里，有兩種辦法可以得到更好的結(jié)果，而這個(gè)我們?cè)诤竺鎸?huì)講到。保持關(guān)注。