R語言之縱向數據分析：多級增長模

上一次，我們討論了如何對長型數據轉換成長型的數據，同時還是用了一個隨機創建的對照實驗數據集來對其增長趨勢進行可視化。但是，我們是否能夠進一步的分析并預測結果的增長趨勢與時間之間的關系。
是的，當然可以！我們可以使用多級增長模型（也稱之為層次模型或者混合模型）進行估計。

產生一個水平數據集并把它轉成寬格式
下面，我們先從我之前的一篇文章的實例進行講解：

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

多級增長模型

這里有很多R語言包可以幫助你進行多級分析，其中，我發現lme4包是最好的一個，因為它使用比較簡單，而且建模能力也很強（尤其是輸出二進制結果或者計數結果）。當然，nlme包也是相當不錯的，它可以給連續型結果提供了類似的結果（正態/高斯分布）。

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

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

含有p估計值的多級增長模型

下面大部分的代碼和上面的類似，除非我們要使用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)

其結果很相似，但現在，我們可以得到自由度和p的估計值。所以，我們可以很自信的說普通RCT參與治療的人，現在，隨著時間的推移，他們的抑郁癥得分在下降，其速度為每下降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

計算95%置信區間和預測區間

有時，我們想在單個軌跡的均值進行作圖。如果要展示均值里的一些不確定因素，我們需要使用擬合好的模型，利用擬合值進行計算，算出95%置信區間和95%預測區間。

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

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

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%置信區間和95%預測區間的圖像了。注意，預測區間的圖像要寬于置信區間。也就是說，預測均值的結果比用單個值預測要好。

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

如果你和我一樣，對數據也很敏感，你應該能觀察到圖線的擬合效果并不太好。這里，有兩種辦法可以得到更好的結果，而這個我們在后面將會講到。保持關注。