林存潔，吳 朦，易丹輝，胡鏡清＊＊

（1.中國(guó)人民大學(xué)應(yīng)用統(tǒng)計(jì)科學(xué)研究中心，北京100872；2.中國(guó)人民大學(xué)統(tǒng)計(jì)學(xué)院，北京 100872；3.中國(guó)中醫(yī)科學(xué)院中醫(yī)基礎(chǔ)理論研究所，北京 100700）

多個(gè)響應(yīng)變量的縱向數(shù)據(jù)聯(lián)合建模方法及應(yīng)用＊

林存潔1,2，吳 朦3，易丹輝1,2，胡鏡清3＊＊

（1.中國(guó)人民大學(xué)應(yīng)用統(tǒng)計(jì)科學(xué)研究中心，北京100872；2.中國(guó)人民大學(xué)統(tǒng)計(jì)學(xué)院，北京 100872；3.中國(guó)中醫(yī)科學(xué)院中醫(yī)基礎(chǔ)理論研究所，北京 100700）

縱向觀測(cè)過程中經(jīng)常會(huì)遇到對(duì)同一個(gè)個(gè)體測(cè)量多個(gè)響應(yīng)變量的情況，由于同一個(gè)個(gè)體的不同響應(yīng)變量之間存在一定的相關(guān)性，因此獨(dú)立分析每一個(gè)響應(yīng)變量將會(huì)損失相關(guān)信息。本文利用聯(lián)合建模的方式對(duì)多個(gè)響應(yīng)變量建立混合效應(yīng)模型，通過不同響應(yīng)變量的隨機(jī)效應(yīng)之間的相關(guān)性刻畫不同響應(yīng)變量之間的關(guān)系，并利用極大似然方法給出模型中參數(shù)的估計(jì)，最后通過中風(fēng)病的實(shí)際數(shù)據(jù)分析來(lái)說明該方法的應(yīng)用。

多元縱向數(shù)據(jù) 線性混合效應(yīng)模型 聯(lián)合模型 隨機(jī)效應(yīng)

在臨床試驗(yàn)、生物醫(yī)學(xué)和社會(huì)科學(xué)等研究領(lǐng)域中，經(jīng)常會(huì)對(duì)不同的個(gè)體在一段時(shí)間內(nèi)進(jìn)行重復(fù)測(cè)量，從而得到縱向觀測(cè)數(shù)據(jù)。雖然不同個(gè)體之間的觀測(cè)彼此獨(dú)立，但是同一個(gè)個(gè)體的不同觀測(cè)之間存在一定的相關(guān)性，因此縱向數(shù)據(jù)分析的一個(gè)重要研究?jī)?nèi)容就是如何考慮個(gè)體內(nèi)的相關(guān)性。其中，最常用的統(tǒng)計(jì)模型之一是線性混合效應(yīng)模型[1]，該模型通過引入隨機(jī)效應(yīng)來(lái)刻畫個(gè)體在不同觀測(cè)點(diǎn)之間的關(guān)系；后來(lái)該模型的思想被拓展到非線性混合效應(yīng)模型[2]以及廣義線性混合效應(yīng)模型[3]等。近年來(lái)這些模型被廣泛應(yīng)用于單個(gè)響應(yīng)變量的縱向數(shù)據(jù)分析中[4,5,6]。

目前關(guān)于縱向數(shù)據(jù)的研究，大部分集中于單個(gè)響應(yīng)變量的分析方法。但是在很多情況下，單個(gè)響應(yīng)變量不能夠全面描述個(gè)體的情況，因此對(duì)同一個(gè)個(gè)體會(huì)測(cè)量?jī)蓚€(gè)或多個(gè)響應(yīng)變量，即得到多個(gè)響應(yīng)變量的縱向觀測(cè)數(shù)據(jù)。例如，在對(duì)中風(fēng)病氣虛血瘀證的中醫(yī)復(fù)雜干預(yù)研究中，需要評(píng)價(jià)中醫(yī)介入治療和單純西醫(yī)治療的療效差異，由于中風(fēng)病的治療方式多樣，常常有“多作用、多靶點(diǎn)”的效果，因此對(duì)于中風(fēng)病的療效評(píng)價(jià)指標(biāo)有很多種。從中醫(yī)角度出發(fā)，有1986年制定的《中風(fēng)病中醫(yī)診斷、療效評(píng)定標(biāo)準(zhǔn)》、1994年制定的《中風(fēng)病辨證診斷標(biāo)準(zhǔn)（試行）》、1996年制定的《中風(fēng)病診斷與療效評(píng)定標(biāo)準(zhǔn)（試行）》以及2011年制定的《缺血性中風(fēng)證候要素診斷量表》。另一方面，現(xiàn)代醫(yī)學(xué)對(duì)腦卒中的評(píng)價(jià)指標(biāo)又有神經(jīng)功能缺失評(píng)定指標(biāo)（NIHSS）、功能預(yù)后指標(biāo)（Brunnstrom）等。對(duì)于同一個(gè)患者而言，不同的指標(biāo)反映不同方面的信息，因此不同的指標(biāo)之間必然存在一定的相關(guān)信息。如何準(zhǔn)確刻畫和利用不同響應(yīng)變量之間的相關(guān)信息來(lái)提高估計(jì)的準(zhǔn)確性，是多個(gè)響應(yīng)變量的縱向數(shù)據(jù)分析研究的重點(diǎn)。

近年來(lái)，關(guān)于多個(gè)響應(yīng)變量的縱向數(shù)據(jù)分析方法也成為研究的熱點(diǎn)。一種方式是將多個(gè)縱向觀測(cè)響應(yīng)變量降維成一個(gè)響應(yīng)指標(biāo)變量，然后利用單個(gè)響應(yīng)變量的縱向數(shù)據(jù)分析方法進(jìn)行分析。該方法容易理解和實(shí)現(xiàn)，但是需要選擇合適的降維方法，而最常用的主成分分析的方法不能完全提取原始數(shù)據(jù)的信息，另外降維之后的指標(biāo)不一定有很好的解釋意義。另外一種方式是通過引入潛在變量的方式將多個(gè)響應(yīng)變量聯(lián)系起來(lái)[7,8,9,10]，并利用EM算法實(shí)現(xiàn)模型的估計(jì)，該方法需要基于一些關(guān)于潛在變量的假設(shè)，但是這些假設(shè)是難以驗(yàn)證的。第三種方式則是基于混合效應(yīng)模型的聯(lián)合建模方法[11,12,13,14,15]，該方法通過對(duì)每個(gè)響應(yīng)變量的縱向觀測(cè)過程假設(shè)一個(gè)隨機(jī)效應(yīng)，進(jìn)而利用不同隨機(jī)效應(yīng)之間的聯(lián)合分布刻畫不同響應(yīng)過程之間的關(guān)系。本文將采用一般混合效應(yīng)聯(lián)合建模方法對(duì)多個(gè)響應(yīng)變量的縱向觀測(cè)數(shù)據(jù)進(jìn)行分析，并給出估計(jì)方法。

本文的第二部分介紹了混合效應(yīng)的聯(lián)合模型和估計(jì)方法；在第三部分將方法應(yīng)用于中風(fēng)病的療效評(píng)價(jià)；第四部分是本文的總結(jié)和討論。

1 多個(gè)響應(yīng)變量的聯(lián)合模型及其估計(jì)

我們可以通過向量拉直的方式將模型（1）轉(zhuǎn)化成一般的縱向數(shù)據(jù)混合效應(yīng)模型。令為個(gè)體i的K個(gè)響應(yīng)指標(biāo)的觀測(cè)向量，為niK×1維的向量，為pK×1維固定效應(yīng)向量，bi=為qK×1維隨機(jī)效應(yīng)向量，定義=diag(Xi,…,Xi)為 niK×pK 的矩陣，=diag(Zi,…,Zi)為niK×qK的矩陣，則模型（1）等價(jià)于下列線性混合效應(yīng)模型：

對(duì)于一般的線性混合效應(yīng)模型的估計(jì)，文獻(xiàn)中已經(jīng)有很多方法，如Laird&Ware（1982）介紹了極大似然估計(jì)方法，并給出EM算法，同時(shí)給出了帶有限制的極大似然估計(jì)（RMLE）和相應(yīng)的計(jì)算方法。對(duì)于模型（2），我們應(yīng)用類似的估計(jì)方法。具體來(lái)說，因?yàn)閅i～，因此我們很容易寫出對(duì)數(shù)似然函數(shù)：

其中c為常數(shù)項(xiàng)。如果協(xié)方差陣Vi是已知的，則可以利用極大似然估計(jì)得到β和bi的估計(jì)：

其中Wi=。但是協(xié)方差陣往往是未知的，即和Ri中含有一些未知參數(shù)（我們用θ表示），那么我們可以關(guān)于 β和θ同時(shí)極大化似然函數(shù)，即得到極大似然估計(jì)。進(jìn)而可以利用經(jīng)驗(yàn)貝葉斯的方法得到的估計(jì)

可以證明該估計(jì)是相合的，而且是漸近正態(tài)的，其估計(jì)的方差可以用如下方式進(jìn)行估計(jì)：而且當(dāng)隨機(jī)效應(yīng)的聯(lián)合分布可以準(zhǔn)確刻畫不同的響應(yīng)變量之間的相關(guān)性時(shí)，聯(lián)合建模的方式可以有效提高估計(jì)的準(zhǔn)確性。另外，當(dāng)K的值比較大的時(shí)候，計(jì)算量會(huì)比較大，可以利用成對(duì)似然估計(jì)的方[16]對(duì)參數(shù)進(jìn)行估計(jì)。

圖1 313例患者的NIHSS評(píng)分和Brunnstrom評(píng)分均值走向圖

2 實(shí)例分析

在本節(jié)中我們將上述方法應(yīng)用到中風(fēng)病氣虛血瘀證中醫(yī)復(fù)雜干預(yù)的研究中。該研究選取313例中風(fēng)?。毙匀毖阅X卒中）患者作為研究示范病例，旨在分析中西醫(yī)聯(lián)合治療的效果。數(shù)據(jù)中記錄了患者的基本信息和中西醫(yī)日常診療信息?；颊弑浑S機(jī)分配到試驗(yàn)組和對(duì)照組，對(duì)照組的患者采用西醫(yī)標(biāo)準(zhǔn)治療，試驗(yàn)組的患者在原來(lái)西醫(yī)標(biāo)準(zhǔn)治療的基礎(chǔ)上，加用丹紅或參麥注射液。所采用的療效評(píng)價(jià)指標(biāo)有美國(guó)國(guó)立衛(wèi)生研究院卒中量表（National Institute of Health Stroke Scale，NIHSS）、Barthel指數(shù)（Barthel Index，BI）、改良Rankin量表（Modified Rankin Scale，MRS），Brunnstrom Scale評(píng)分以及中醫(yī)證候評(píng)分（參照1994年的《中風(fēng)病辨證診斷標(biāo)準(zhǔn)（試行）》）等。我們以NIHSS和Brunnstrom評(píng)分為例，分析不同治療方法的效果。首先對(duì)NIHSS評(píng)分和Brunnstrom評(píng)分作均值圖，如圖1。圖中Treatment=1表示試驗(yàn)組，Treatment=0表示對(duì)照組。從該圖可以看出，觀測(cè)初始點(diǎn)時(shí)試驗(yàn)組的兩種評(píng)分都略高于對(duì)照組的評(píng)分，隨著治療時(shí)間的推移兩種評(píng)分都在下降，但是最后一次觀測(cè)時(shí)間點(diǎn)的對(duì)照組的評(píng)分卻略高于試驗(yàn)組的評(píng)分，尤其是Brunnstrom評(píng)分，最后的差異更大。這說明中西醫(yī)結(jié)合治療方案的治療效果可能會(huì)更好。

為了更加準(zhǔn)確地分析治療方法的效果，利用聯(lián)合建模方法對(duì)NIHSS評(píng)分和Brunnstrom評(píng)分進(jìn)行綜合分析，建立如下聯(lián)合模型：

這里考慮的協(xié)變量Xi有時(shí)間、治療方案分組（0=對(duì)照組，1=實(shí)驗(yàn)組）、時(shí)間和治療方案分組的交互項(xiàng)、性別（0=男性，1=女性）、年齡、BMI指數(shù)、合并疾?。ㄟ@里指合并三高和冠心病的情況，0=未合并其他疾病，1=入組時(shí)已合并其他疾病）、中風(fēng)的初發(fā)與復(fù)發(fā)情況（0=初發(fā)，1=復(fù)發(fā)）、氣虛血瘀證評(píng)分是否≥7分（0=否，1=是），考慮的隨機(jī)效應(yīng)的協(xié)變量Zi為時(shí)間變量。具體估計(jì)結(jié)果見表1。從表1可以看出，對(duì)于NIHSS和Brunnstrom評(píng)分而言，時(shí)間的系數(shù)估計(jì)值分別為-1.517和-0.996，且估計(jì)的p值都很小（＜0.05），這說明隨著治療時(shí)間的推移，對(duì)照組患者的NIHSS評(píng)分和Brunnstrom評(píng)分都會(huì)顯著地降低，即西醫(yī)治療方案能夠緩解中風(fēng)病情。治療方案分組的估計(jì)系數(shù)分別為0.03和1.124，p值為0.926和0.061，這表示在觀測(cè)初始點(diǎn)時(shí)試驗(yàn)組患者的NIHSS評(píng)分均值和Brunnstrom評(píng)分的均值比對(duì)照組患者的NIHSS評(píng)分均值和Brunnstrom評(píng)分的均值要高，但差異并不顯著。分組*時(shí)間的系數(shù)估計(jì)值分別為-0.018和-0.416，這說明隨著時(shí)間的變化，試驗(yàn)組患者的NIHSS評(píng)分和Brunnstrom評(píng)分分別比對(duì)照組的NIHSS評(píng)分和Brunnstrom評(píng)分降低地更快，但是NIHSS評(píng)分的差異不顯著（p值為0.873），而Brunnstrom評(píng)分的差異顯著（p值為0.010）。這說明中西醫(yī)結(jié)合治療方案能夠更快的緩解病情，這些結(jié)論和圖1的結(jié)論也相吻合。除此之外，其它變量的系數(shù)估計(jì)結(jié)果可以看出，對(duì)于NIHSS評(píng)分而言，在其他影響因素相同的情況下，女性患者比男性患者的NIHSS評(píng)分要高，年齡越大的患者NIHSS評(píng)分也會(huì)越高，肥胖的患者評(píng)分也會(huì)越高，但是這些差異不顯著；合并疾病是一個(gè)顯著的影響因素，即有其他合并疾病的患者比未合并其他疾病的患者NIHSS評(píng)分的均值要低；另外復(fù)發(fā)患者以及氣虛血瘀證嚴(yán)重的患者其NIHSS評(píng)分也會(huì)越高，不過差異并不顯著。對(duì)于Brunnstrom評(píng)分，顯著的影響因素是年齡，即老年患者的評(píng)分要顯著的高于其它患者的評(píng)分，而其他因素都不顯著。

表1 NIHSS與Brunnstrom聯(lián)合建模估計(jì)結(jié)果（313例）

表2 隨機(jī)效應(yīng)的相關(guān)性和估計(jì)

為了檢驗(yàn)NIHSS評(píng)分和Brunnstrom評(píng)分之間的相關(guān)性，我們檢驗(yàn)隨機(jī)效應(yīng)的顯著性以及不同響應(yīng)過程之間隨機(jī)效應(yīng)的相關(guān)性。表2展示了混合效應(yīng)模型中隨機(jī)效應(yīng)的方差估計(jì)和相關(guān)系數(shù)估計(jì)，從該表的結(jié)果可以看出NIHSS評(píng)分的隨機(jī)效應(yīng)和Brunnstrom評(píng)分的隨機(jī)效應(yīng)都是顯著不為0（置信上下限不包含0），其中和表示隨機(jī)效應(yīng)的截距項(xiàng)，和表示時(shí)間的隨機(jī)效應(yīng)。而且NIHSS評(píng)分的時(shí)間隨機(jī)效應(yīng)和Brunnstrom評(píng)分的隨機(jī)效應(yīng)具有顯著的線性相關(guān)性，NIHSS評(píng)分的截距項(xiàng)隨機(jī)效應(yīng)和Brunnstrom評(píng)分的截距項(xiàng)隨機(jī)效應(yīng)也具有顯著的線性相關(guān)性。圖2展示了隨機(jī)效應(yīng)的估計(jì)散點(diǎn)圖，從圖中也得出相同的結(jié)論。這也說明了NIHSS評(píng)分和Brunnstrom評(píng)分之間確實(shí)存在一定的相關(guān)性。

3 小結(jié)

圖2 NIHSS評(píng)分和Brunnstrom評(píng)分的隨機(jī)效應(yīng)圖

本文針對(duì)多個(gè)響應(yīng)變量的縱向觀測(cè)數(shù)據(jù)的建模方法進(jìn)行了討論。由于同一個(gè)個(gè)體的不同響應(yīng)變量之間存在一定的相關(guān)性，如果分別分析每個(gè)響應(yīng)變量而忽略不同響應(yīng)變量之間的關(guān)系，可能會(huì)得到錯(cuò)誤的結(jié)論。因此，我們提出利用聯(lián)合建模的方式對(duì)此類數(shù)據(jù)進(jìn)行分析，通過對(duì)每個(gè)響應(yīng)變量建立一般線性混合效應(yīng)模型，并利用隨機(jī)效應(yīng)之間的相關(guān)性刻畫不同的響應(yīng)變量之間的關(guān)系，進(jìn)而充分利用數(shù)據(jù)中的信息提高估計(jì)的精確性。為了估計(jì)模型中的參數(shù)，將數(shù)據(jù)拉直得到一般的線性混合效應(yīng)模型，然后利用極大似然方法給出模型中參數(shù)的估計(jì)。最后通過中風(fēng)病的數(shù)據(jù)說明該方法的應(yīng)用。如果單獨(dú)對(duì)NIHSS評(píng)分建立線性混合效應(yīng)模型會(huì)發(fā)現(xiàn)，顯著的變量有時(shí)間、年齡和合并疾病，如果單獨(dú)對(duì)Brunnstrom評(píng)分建立線性混合效應(yīng)模型，顯著的變量則是時(shí)間和年齡。但是通過NIHSS評(píng)分和Brunnstrom評(píng)分的聯(lián)合建模發(fā)現(xiàn)分組*時(shí)間是Brunnstrom評(píng)分的顯著變量，這也說明試驗(yàn)組的患者病情緩解地更快。因此借助于聯(lián)合建模的方式可以更加準(zhǔn)確的評(píng)價(jià)不同治療方式的效果。除了考慮兩個(gè)響應(yīng)變量的聯(lián)合建模過程，我們還可以考慮NIHSS評(píng)分、Brunnstrom評(píng)分和其它評(píng)價(jià)指標(biāo)的聯(lián)合建模，如Barthel指數(shù)、中醫(yī)證候評(píng)分等。如果響應(yīng)變量的個(gè)數(shù)比較多，我們可以考慮利用成對(duì)似然方法進(jìn)行估計(jì)。

1 Laird N,Ware J.Random-effects models for longitudinal data.Biometrics,1982,38:963-974.

2 Davidian M,Giltinan D.Nonlinear models for repeated measurement data.New York:Chapman&Hall,1995.

3 Breslow N,Clayton D.Approximate inference in generalized linear mixed models.JAm Stat Assoc.1993,88:9-25.

4 Fitzmaurice G,Laird N,Ware J.Applied Longitudinal Analysis.Wiley,New York,2004.

5 Hedeker D.,Gibbons R.D.Longitudinal Data Analysis.Wiley,Hoboken,New Jersey,2006.

6 Diggle P J,Heagerty P,Liang K Y,et al.Analysis of Longitudinal Data.New York:Clarendon Press,2002.

7 Chan D.The conceptualization and analysis of changeover time:An integrative approach incorporating longitudinal mean and covariance structures analysis(LMACS)and multiple indicator latent growth modeling(MLGM).Organiz Res Meth.1998,1:421-483.

8 Roy J,Lin X.Latent variable models for longitudinal data with multiple continuousoutcomes.Biometrics,2000,56:1047-1054.

9 Blozis S.A second-order structured latent curve model for longitudinal data.In:van Montfort K,Oud H and Satorra A(eds)Longitudinal models in the behavioural and related sciences.Mahwah,New Jersey:Lawrence Erlbaum Associates,2006:189-214.

10 Harring J.A nonlinear mixed effects model for latent variables.J Educ Behav Stat,2009,34:293-318.

11 Shah A,Laird N,Schoenfeld D A random-effects model for multiple characteristics with possibly missing data.J Am Stat Assoc,1997,92:775-779.

12 Schafer JL,Yucel RM.Computational strategies for multivariate linear mixed effects modelswith missing values.JComput Graph Statist,2002,11,437-457.

13 Fieuws S,Verbeke G.Joint modelling of multivariate longitudinal profiles:Pitfalls of the random-effects approach.Stat Med,2004,23:3093-3104.

14 Roy A.Estimating correlation coefficient between two variables with repeated observations using mixed effects model.Biom J.2006,48:286-301.

15 Wang W,Fan T.Estimation in multivariate t linear mixed models for multiplelongitudinal data.Stat Sinica,2011,21:1857-1880.

16 Fieuws S.,Verbeke G.Pairwise fitting of mixed models for the joint modeling of multivariate longitudinal profiels.Biometrics.2006,62:424-431.

Joint Modeling of Multivariate Longitudinal Data and Its Application

Lin Cunjie1,2,Wu Meng3,Yi Danhui1,2,Hu Jingqing3
(1.Center for Applied Statistics,Renmin University of China,Beijing 100872,China;2.School of Statistics,Renmin University of China,Beijing 100872,China;3.Instituteof Basic Theory for Chinese Medicine,China Academy of Chinese Medical Sciences,Beijing 100700,China)

Multiple outcomes measured repeatedly for the same subject are common in longitudinal observation.If we use the approach by analyzing each outcome separately,it may lead to wrong conclusions due to the failure of accounting for joint evolution of different outcomes.To adequately capture the interdependence among multiple outcomes,we proposed a joint modeling for multivariate longitudinal data by constructing a linear mixed-effects model for each outcome and accommodating the relationship among multiple outcomes through correlation in random effects.Maximum likelihood method was adopted to estimate parameters in this model.The application of this method was demonstrated through the analysisof strokedata.

Multiplelongitudinal data,linear mixed-effectsmodel,joint modeling,randomeffects

10.11842/wst.2017.09.006

R33

A

2017-05-23

修回日期：2017-07-25

＊ 中國(guó)人民大學(xué)2017年度中央高校建設(shè)世界一流大學(xué)（學(xué)科）和特色發(fā)展引導(dǎo)專項(xiàng)資金（297217000021），負(fù)責(zé)人:易丹輝；教育部人文社會(huì)科學(xué)重點(diǎn)研究基地重大項(xiàng)目（16JJD910002）：基于大數(shù)據(jù)的精準(zhǔn)醫(yī)學(xué)生物統(tǒng)計(jì)分析方法及其應(yīng)用研究，負(fù)責(zé)人：易丹輝；國(guó)家中醫(yī)藥管理局中醫(yī)行業(yè)科研專項(xiàng)（201207005）：30種疾病中醫(yī)臨床評(píng)價(jià)規(guī)范與復(fù)雜干預(yù)評(píng)價(jià)共同路徑研究，負(fù)責(zé)人：胡鏡清。

＊＊ 通訊作者：胡鏡清，研究員，博士生導(dǎo)師，主要研究方向：適應(yīng)中醫(yī)藥理論構(gòu)筑與診療模式的臨床研究方法研究。

（責(zé)任編輯：張娜娜，責(zé)任譯審：王 晶）