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

這篇文章已經(jīng)是縱向數(shù)據(jù)分析系列的第三篇了。之前，我們介紹了什么是縱向數(shù)據(jù)，我們?nèi)绾伟验L(zhǎng)型的數(shù)據(jù)集轉(zhuǎn)換成縱向的數(shù)據(jù)集，并通過建立相應(yīng)的多級(jí)模型進(jìn)行分析。顯然，僅僅介紹基本的多級(jí)模型并不足以對(duì)虛擬隨機(jī)對(duì)照實(shí)驗(yàn)數(shù)據(jù)集進(jìn)行分析。這篇文章則延續(xù)之前的分析，并介紹一種不錯(cuò)的方法來處理多級(jí)模型里的治療效果問題。

創(chuàng)建一個(gè)縱向數(shù)據(jù)集并把它轉(zhuǎn)成長(zhǎng)型格式

和往常一樣，我們還是先創(chuàng)建一個(gè)縱向數(shù)據(jù)集，然后把它轉(zhuǎn)成長(zhǎng)型的格式。

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)模型分析治療效果

這是一個(gè)RCT數(shù)據(jù)集，也就說明，這里潛在存兩個(gè)以療效分組的差異。上一次，我們沒有考慮異質(zhì)性，而僅僅對(duì)兩個(gè)組的普通療效進(jìn)行了細(xì)分。從直覺上看，我們可以添加治療變量來修復(fù)這個(gè)模型，從而抓到組間的差異。

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

這里，我們?nèi)稳皇褂胠merTest包是因?yàn)樗试S修正模型的檢驗(yàn)使用自由度。模型的公式和是和上次的一樣，除非我們把治療變量tx作為獨(dú)立變量。

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這兩個(gè)變量都相當(dāng)顯著，而接受B治療的人群在抑郁癥得分上比接受A治療的人群低2.34分。那么，我們是否可以認(rèn)為，B治療方案效果更好？絕對(duì)不是。模型m1實(shí)際上并不合適對(duì)縱向數(shù)據(jù)集里的組間差異進(jìn)行比較。m1里的治療效果，其實(shí)也就是一段時(shí)間內(nèi)的平均治療效果。換句話說，治療方案A和治療方案B的療效所導(dǎo)致的分差在-2.34分，其實(shí)只不過是1,2,3時(shí)間段內(nèi)抑郁癥均分的差。這是RCT，所以，我們所期待的志愿者參與兩組治療方案里時(shí)間段1的差要小于1。我們對(duì)兩組療效的均值差異感興趣，而現(xiàn)在的軌跡在時(shí)間段2、3里存在差異。
這個(gè)概率可能對(duì)初學(xué)者來說有點(diǎn)復(fù)雜，但是其執(zhí)行過程很簡(jiǎn)單，我們只需要添加交互變量time和tx。

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

在R，我們可以使用:符號(hào)來表明變量間的交互關(guān)系。因此，我們可以這樣設(shè)置：time+tx+time:tx。但是，更簡(jiǎn)短的寫法是：timetx。操作符要R涵蓋主要變量和治療療效。此時(shí)，使用操作符就能簡(jiǎn)單實(shí)現(xiàn)各個(gè)效應(yīng)的展示。請(qǐng)記住，操作符也可用于高維變量交互。比如，ABC就要R包含三個(gè)主要的效應(yīng)，三個(gè)兩兩交互的效應(yīng)，以及一個(gè)三元交互效應(yīng)。

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

現(xiàn)在，我們從m2模型中獲得一個(gè)三元交互關(guān)系。在這個(gè)模型當(dāng)中，time效應(yīng)指的就是療效A（這是特定條件下的療效和療效時(shí)間的交互），而txB效應(yīng)就是時(shí)間點(diǎn)為0時(shí)的效應(yīng)（數(shù)據(jù)集里一個(gè)虛的時(shí)間點(diǎn)）。這個(gè)術(shù)語解釋起來有點(diǎn)復(fù)雜：系數(shù)值反映了每個(gè)時(shí)間段內(nèi)兩種治療方案的療效上的差異，其條件是取平均時(shí)間和治療效果?？陬^上說，人們可以理解的是，這表明兩組治療方案是如何隨著時(shí)間的推移而不同。對(duì)于每個(gè)單位時(shí)間內(nèi)的增長(zhǎng)（即，從時(shí)間段1到時(shí)間段2），在完成平均治療效應(yīng)差的分析以后，接受B治療方案的參與者理論上在抑郁癥得分方面比接受A治療方案的要低2.88分。所以，在時(shí)間段2內(nèi)，接受B治療方案的參與者得分差是 -2.88 x 2 + 3.43 = -2.33，而在時(shí)間段3內(nèi)，接受B治療方案的參與者得分差是-2.88 x 3+ 3.43 = -5.21。我們應(yīng)當(dāng)時(shí)刻注意，在一個(gè)給定的時(shí)間段里對(duì)比療效差的時(shí)候，主要療效以及療效變量間的交互。但是，為什么要這樣？你不妨設(shè)置分析一下，當(dāng)你忽略了主要治療效應(yīng)的時(shí)候，你在時(shí)間段1內(nèi)所得到的療效差是怎樣的情形。

選擇最后的模型

我可以告訴你，模型m2比m1和m0（上一篇文章所采用的模型）都要好，原因在于，這個(gè)數(shù)據(jù)集是我們自己創(chuàng)建的。但在實(shí)際生活中，我們會(huì)在選擇最好的一個(gè)模型的時(shí)候會(huì)遇到很多困難。其中一種選擇的方式就是，按照我們之前所講的，使用summary()函數(shù)，以此算出額外變量的統(tǒng)計(jì)顯著性。其它一個(gè)常用的方法就是使用方差分析（ANOVA)表（僅僅針對(duì)嵌套模型），而這種方法適合用于這個(gè)模型。

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()對(duì)m0、m1和m2模型進(jìn)行方差分析，原因在于我們把m0嵌套于m1，然后把m1嵌套于m2。如果還有另外的子集，我們就把模型嵌套到那個(gè)子集上。換句話說，含有獨(dú)立變量A、B的模型則嵌套于含有變量A、B、C的模型中。
如果你做過方差分析，那這個(gè)方差分析表對(duì)于你來說應(yīng)該是很熟悉了。Df列表明了模型里的自由度，它簡(jiǎn)單的反映了這個(gè)案例的參數(shù)估計(jì)。例如，自由度為4的模型m0表明了這個(gè)模型預(yù)測(cè)了4個(gè)參數(shù)（1.修正截距效應(yīng)系數(shù)，2.時(shí)間效應(yīng)，3.隨機(jī)截距方差，4.殘差）。
AIC和BIC是兩個(gè)常用的擬合指數(shù)。AIC和BIC都考慮到兩個(gè)因素：模型的數(shù)據(jù)擬合效果有多好，和這個(gè)模型復(fù)不復(fù)雜。AIC和BIC的不同之處在于模型復(fù)雜度的罰分上（BIC對(duì)模型復(fù)雜度的罰分影響更大）。我發(fā)現(xiàn)，AIC理解起來更容易，因?yàn)樗鼭u進(jìn)的使用留一交叉檢驗(yàn)法（LOOCV)來預(yù)測(cè)運(yùn)行效果。我們應(yīng)當(dāng)理解交叉檢驗(yàn)和預(yù)測(cè)性能的概念，隨后，如果你對(duì)此不太理解，你可以先把它跳過。

LogLik列模型參數(shù)的似然對(duì)數(shù)。從根本上說，它是模型里數(shù)據(jù)擬合效果的指標(biāo)。其偏差是?2×logLik?2×logLik。為什么我們需要它？因?yàn)閮蓚€(gè)嵌套模型的偏差在原假設(shè)的條件下遵循卡方檢驗(yàn)（假設(shè)兩個(gè)嵌套模型的偏差一樣）。此時(shí)，我們可以把它用作模型擬合差的測(cè)試。偏差的差則在Chisq列顯示?？ǚ綑z驗(yàn)的統(tǒng)計(jì)相關(guān)自由度則在Chi Df（簡(jiǎn)單來說，就是額外參數(shù)的個(gè)數(shù)）。最后一行的結(jié)果顯而易見：它算出了測(cè)試兩個(gè)嵌套模型的差的p值。雖然是這樣，我們還是要警惕ANOVA的結(jié)果，因?yàn)樵诒姸嗟暮蜻x模型中，它會(huì)導(dǎo)致錯(cuò)誤類型1上升極快。

由于模型m2比其它兩個(gè)候選模型AIC和BIC都要好，而對(duì)于方差分析，我們可以認(rèn)為這個(gè)模型是目前為止我們擬合的最好的了。當(dāng)然，這個(gè)模型還有許多要改進(jìn)的地方。后續(xù)，我們可能還可以在模型的性能上進(jìn)行提升，并作圖進(jìn)行預(yù)測(cè)。