縱向數(shù)據(jù)與生存數(shù)據(jù)聯(lián)合模型中多變點識別問題

沈 佳 坤，　宋 立 新，　孫 秀 峰，　馮 寶 軍

( 1.大連理工大學 管理與經(jīng)濟學部, 遼寧 大連　116024；2.大連理工大學 數(shù)學科學學院, 遼寧 大連　116024 )

?

縱向數(shù)據(jù)與生存數(shù)據(jù)聯(lián)合模型中多變點識別問題

沈 佳 坤1,2，宋 立 新*2，孫 秀 峰1，馮 寶 軍1

( 1.大連理工大學 管理與經(jīng)濟學部, 遼寧 大連116024；2.大連理工大學 數(shù)學科學學院, 遼寧 大連116024 )

提出了共享協(xié)變量和隨機效應的縱向響應中含有多個變點識別的線性混合效應(LME)模型和加速失效時間(AFT)模型的聯(lián)合模型，并通過Gauss-Hermite近似解決極大似然函數(shù)中的復雜積分以得到參數(shù)的估計．通過模擬研究驗證了該方法的有效性，并將其應用于原發(fā)性膽汁性肝硬化(PBC)病變過程，研究發(fā)現(xiàn)：PBC患者的血清膽紅素只在初期治療階段有所降低，兩個月之后迅速開始反彈，直到3.5 a后增速才有所放緩，說明治療方法仍需改進．

多變點；線性混合效應模型；加速失效時間模型；聯(lián)合推斷；極大似然

0　引　言

縱向數(shù)據(jù)既包含橫截面數(shù)據(jù)中多個體間的對比差異，又包含時間序列數(shù)據(jù)中對個體重復測量的變化趨勢，能較好地反映研究對象的動態(tài)變化特征，常被用來評定誘導危害人們身體健康的風險因素的影響、描述患者的生活改善以及評估治療的效率．此外，生存分析揭示了患者從開始治療到疾病復發(fā)或患者死亡等重要臨床事件的時間變化規(guī)律，是醫(yī)學統(tǒng)計研究中的另一個重要內(nèi)容，已取得了很多成果．然而，當縱向響應回歸模型與失效時間回歸模型之間存在某種相關關系時，如體現(xiàn)在一些共同的潛在過程或協(xié)變量上，分別獨立地進行縱向數(shù)據(jù)和生存分析研究會導致結(jié)果產(chǎn)生偏差．綜上，為了將所有影響研究結(jié)果的信息都利用起來得到無偏估計結(jié)果，建立縱向響應和生存時間的聯(lián)合模型勢在必行，并已得到學界的普遍認同[1-2]．

在實際研究中，患者的病情隨時間的發(fā)展往往不是單純的線性關系，而是呈現(xiàn)出階段性的變化特點．為此，學者們一直在努力尋找更簡潔、合理的非線性模型來描述患者病情變化，其中變點回歸以其簡潔的形式、較好的擬合度且具有在變點處連續(xù)的良好性質(zhì)而深受學者們的喜愛，被較多地應用于病情變化波動大的流行病學和癌癥研究中，來識別疾病發(fā)展過程中醫(yī)學指標變化規(guī)律發(fā)生改變的時間點[3-6]．最近，Huang等[7]提出了兩個變點的分段線性模型，并通過貝葉斯的方法得到變點的分布，來描述HIV病毒載量隨時間的變化情況．但由于貝葉斯方法過于依賴給定的先驗分布，不能準確反映新情況下響應變量的自然變化規(guī)律，導致了研究結(jié)果的局限性．

本研究將極大似然估計法引入分段線性模型，以克服貝葉斯推斷固有的缺點，服務于臨床管理中病人用藥的指導服務及醫(yī)療決策中肝移植時機的選擇等．

1　模型和方法

1.1聯(lián)合模型

考慮一個樣本量為N的樣本，個體間相互獨立．令Yij為個體i在時間tij的響應變量，j=1,…,ni且i=1,…,N．本文建立關于縱向數(shù)據(jù)的含有k個變點的線性混合效應(LME)模型(1)，其中k為任意正整數(shù)，并將k個變點在時間軸上的位置作為未知參數(shù)：τ=(τ1τ2…τk)．

Δi1(tij-τ1)++Δi2(tij-τ2)++…+

Δik(tij-τk)++eij

(1)

其中Δi1=(μ3+ui3)-(μ2+ui2)，Δi2=(μ4+ui4)-(μ3+ui3)，…，Δik=(μk+2+uik+2)-(μk+1+uik+1)；τ1<τ2，…，<τk，j=1，…，ni，i=1，…，N，k≥1．

令Ti是個體i的生存時間，在臨床實際中，它往往與個體i的縱向數(shù)據(jù)模型(1)有關．為體現(xiàn)這種關聯(lián)性，研究中假設生存數(shù)據(jù)模型與縱向數(shù)據(jù)模型(1)共享協(xié)變量向量Zi和隨機效應向量ui．特別的，建立混合效應加速失效時間(AFT)模型(2)：

(2)

1.2變點識別

在確定縱向數(shù)據(jù)模型變點個數(shù)階段，本文提出一系列模型的假設檢驗，利用似然比檢驗(LRT)的方法確定LME模型(1)中的變點個數(shù)，使得聯(lián)合模型對樣本的擬合度達到相對最優(yōu)．

然后，本文通過LRT方法確定變點個數(shù)，尋找使聯(lián)合模型具有最優(yōu)擬合度的相應縱向數(shù)據(jù)分段LME模型．第k組LRT方法中檢驗統(tǒng)計量LR(k)的定義如下式所示：

LR(k)=

-2[l(k)0(θ^(k)0)-l(k)1(θ^1(k))]～·

(3)

θ^(k)0

θ^(k)1

1.3聯(lián)合似然推斷

(4)

fTi，δi|Zi，ui(ti，δi|Zi，ui;γ，σ2)fui(ui;G)]dui=

fTi，δi|Zi，ui(ti，δi|Zi，ui;γ，σ2)fui(ui;G)dui=

(5)

式中：K為正交節(jié)點的個數(shù)；us=(us1…us4)T，是節(jié)點值，相應的權重為ws=(ws1…ws4)T．當被積函數(shù)可以被寫作exp(-uTu)l(u)的形式時，合理設置節(jié)點值和相應權重可以得到積分較為準確的估計，其中l(wèi)(u)為階數(shù)小于等于K-1的關于u的多項式，即當節(jié)點數(shù)K足夠大時，近似方程(5)可以無限接近積分的準確值．從而，本文建立的聯(lián)合模型的似然推斷在積分近似意義下是可解的．

2　數(shù)值模擬

2.1模擬方案

不失一般性，本文考慮含有一個隨機效應項和兩個變點的LME模型：

Yij=(μ0+ui)+μ1Zi+μ2tij+(μ3-μ2)×

(tij-τ1)++(μ4-μ3)(tij-τ2)++eij

(6)

其中τ1<τ2，j=1，…，ni，i=1，…，N．

相應的生存分析，建立混合效應AFT模型(7)．生存時間Ti滿足1．3節(jié)中的右刪失假設．

logTi=γ0+γ1Zi+γ2ui+i，

(7)

下面通過一系列數(shù)值模擬研究來驗證上一章中提出的聯(lián)合模型及估計方法的有效性．本節(jié)的數(shù)值模擬基于共享的Z和u的聯(lián)合模型(6)和(7)．假設樣本容量為N=150或N=300，個體間相互獨立，每個個體的重復觀測次數(shù)為ni=6或ni=10．令υ2=σ2=0.12，μ=(0.2-0.3-0.80.90.2)T，且G=0.452．

當ni=6時，為獲得個體i的重復觀測時間，本研究從均勻分布U(0，6)無放回隨機抽取的6個時間點得到次序統(tǒng)計量ti，并建立生存時間的右刪失示性函數(shù)δi=I(Ti≤ti6)，其中ti6是上述抽取均勻分布時間點的最大次序統(tǒng)計量，假設變點τ=(13)，且參數(shù)γ=(1.20.50.2)T；當ni=10時，同理，本文從均勻分布U(0，10)無放回隨機抽取的10個時間點得到次序統(tǒng)計量ti，并建立生存時間的右刪失示性函數(shù)δi=I(Ti≤ti10)，其中ti10是上述抽取均勻分布時間點的最大次序統(tǒng)計量，假設變點τ=(48)，參數(shù)γ=(1．80.50.2)T．同時，假設Zi服從二項分布B(1，0.5)．

從而形成4組縱向數(shù)據(jù)和生存數(shù)據(jù)的數(shù)值模擬方案：

(a)個體數(shù)為150，每個個體重復觀測次數(shù)為6；

(b)個體數(shù)為150，每個個體重復觀測次數(shù)為10；

(c)個體數(shù)為300，每個個體重復觀測次數(shù)為6；

(d)個體數(shù)為300，每個個體重復觀測次數(shù)為10．

2.2數(shù)值結(jié)果

一般的，本研究選擇最大觀測時間的1/3以及2/3分位數(shù)作為極大似然方法估計兩個變點位置的初值，以上4組縱向數(shù)據(jù)和生存數(shù)據(jù)分別作100次數(shù)值模擬．將4種模擬方案下變點(τ1τ2)的估計結(jié)果列示于表1．

表1中的偏差和均方誤差都比原值小4個數(shù)量級以上，表明在該4種數(shù)值模擬方案情形下，極大似然方法可以幾乎準確估計到變點的位置．

表2列示了4種數(shù)值模擬情形下聯(lián)合模型的參數(shù)估計優(yōu)度相關結(jié)果．

通過表2中的普遍低于參數(shù)值本身1～3個數(shù)量級不等的偏差和均方誤差可以驗證參數(shù)估計的無偏性和有效性．此外，注意到偏差和均方誤差隨著重復測量次數(shù)ni由6到10或樣本量N由150到300的增加而降低．

圖1更直觀地呈現(xiàn)了4種方案的數(shù)值模擬，各100次實驗的估計曲線結(jié)果，其中實線表示真實曲線，虛線表示估計曲線．

由圖1中的曲線可以明顯看出重復測量次數(shù)ni由6(圖(a)和(c))到10(圖(b)和(d))，估計曲線與真實曲線之間的差異減?。粯颖玖縉由150(圖(a)和(b))到300(圖(c)和(d))，估計曲線與真實曲線之間的差異明顯減小．

表1　兩個變點位置的數(shù)值模擬結(jié)果

表2　聯(lián)合模型數(shù)值模擬的參數(shù)估計結(jié)果

　　(a)方案(a)

(b) 方案(b)

(c) 方案(c)

(d) 方案(d)

圖14種方案模擬研究曲線

Fig.1The simulation curves of four strategies

3　原發(fā)性膽汁性肝硬化(PBC)實例分析

本文使用的PBC數(shù)據(jù)集來自于美國梅約(Mayo)醫(yī)學研究中心1974年到1984年采集的患者資料[9]，包括312名隨機選取的獨立患者，其中158名患者使用D-青霉胺(D-penicillamine)治療(占51%)，其他154名使用安慰劑(placebo)治療(占49%)．該數(shù)據(jù)集包括患者的基本信息，如年齡、性別等及一些跟蹤測量的生物病理指標，如血清膽紅素(serum bilirubin)、皮膚血管畸形、肝臟腫大等[10]．本文研究的興趣在于反映PBC癥進程的血清膽紅素在不同治療方式下的變化情況，以及相應的生存時間(即從開始治療到死亡或肝移植的時間)的變化特征．

上述美國梅約醫(yī)學研究中心PBC數(shù)據(jù)組，跟蹤記錄的312名隨機選取的患者的情況，包含20個變量，各1 945個觀測值，可以在R語言的JM包中直接獲得．本研究主要關注如下觀測變量：Y代表血清膽紅素serBilir(mg/dL)，是檢測原發(fā)性膽汁性肝硬化程度的指標；t定義為登記日到每次觀測年數(shù)year；Z是用藥情況drug，用藥選擇為D-青霉胺或安慰劑；T定義為登記日到換肝、死亡或研究分析時間的年數(shù)years；δ是生存時間的右刪失示性函數(shù)status2，取值為1時表示患者換肝或死亡，取值為0時表示患者治愈．

本文對這組實例數(shù)據(jù)的分析將分三步進行：(1)確定縱向數(shù)據(jù)模型的變點個數(shù)；(2)估計變點在時間軸上的位置；(3)選擇建立聯(lián)合模型得到相應的參數(shù)估計．

首先，確定縱向數(shù)據(jù)多變點LME模型中變點的個數(shù)．對于樣本量為N的患者間相互獨立的樣本，令Yij為患者i在時間tij的血清膽紅素的測量值．通過1．2節(jié)中假設檢驗規(guī)則，通過原假設和備擇假設中的縱向數(shù)據(jù)模型，建立一系列關于PBC數(shù)據(jù)組的聯(lián)合模型JM0，JM1，…，JMk，檢驗結(jié)果如表3所示．

由表3知，在超過99.99%的置信度水平上含有1個變點的聯(lián)合模型JM1對PBC 數(shù)據(jù)的擬合度優(yōu)于沒有變點的聯(lián)合模型JM0；在超過99.99% 的置信度水平上含有2個變點的聯(lián)合模型JM2對數(shù)據(jù)的擬合度優(yōu)于含有1個變點的聯(lián)合模型JM1；含有3個變點的聯(lián)合模型JM3對數(shù)據(jù)的擬合度不顯著優(yōu)于含有2個變點的聯(lián)合模型JM2．從而含有2個變點的LME模型和相應混合效應AFT模型構成的聯(lián)合模型JM2在本研究中對PBC數(shù)據(jù)擬合最優(yōu)，且AIC和BIC結(jié)果也都支持了該結(jié)論．

表3　PBC聯(lián)合模型似然比檢驗結(jié)果

然后，本文選用如下縱向數(shù)據(jù)模型(8)和生存時間模型(9)通過共享協(xié)變量Zi和隨機效應向量ui構成聯(lián)合模型描述PBC數(shù)據(jù)．

logYij=(μ0+ui1)+μ1Zi+(μ2+ui2)tij+Δi1(tij-τ1)++Δi2(tij-τ2)++eij

(8)

(9)

其中Zi是有兩個取值的等級協(xié)變量，Ti是患者i從登記日到死亡、肝移植或研究分析時間的年數(shù)，相應的生存時間的右刪失示性函數(shù)為δi，參數(shù)向量γ2與隨機效應向量ui的維數(shù)一致．

同上一章模擬研究一樣，本研究選觀測時間變量t最大值的1/3和2/3分位數(shù)即0.988 4和3.983 7，作為極大似然優(yōu)化方法估計兩個變點位置的初值．通過計算可以得到兩個變點在時間軸上的位置(τ1τ2)為(0.141 53.542 2)，即為登記日起的第2個月和第3.5 a．

最后，通過Gauss-Hermite近似意義下的極大似然方法得到聯(lián)合模型中其他參數(shù)的估計值及標準誤，如表4所示．

表4　PBC聯(lián)合模型參數(shù)估計

由表4可知，用藥Z的系數(shù)為-0.306 9，說明D-青霉胺確實對血清膽紅素的升高有抑制作用．通過時間進程的參數(shù)估計及上述兩個變點的估計值可知，從開始用藥到第2個月為第一階段，該階段的血清膽紅素有一個短暫而明顯的降低，增長率為-0.694 7，標準誤較大為0.212 4，并且相應隨機效應的方差估計值也較大為4.416 4，說明由于用藥的不同，這一階段不同個體的血清膽紅素變化差異加大；接著很快開始反彈，血清膽紅素呈現(xiàn)正向增長，增長率為0.172 8，標準誤較小為0.019 0，并且相應隨機效應的方差估計值也較小為0.078 9；直到第3.5 a的時候，血清膽紅素增速變緩，增長率變?yōu)?.130 5直到觀測結(jié)束，標準誤進一步減小為0.013 3，相應隨機效應的方差估計值進一步減小為0.021 6，說明隨著時間的增加，患者抗藥性的增強，個體間的差異逐漸減?。?/p>

圖2顯示了12個隨機選擇的PBC患者血清膽紅素觀測值及擬合曲線．

圖212個隨機選擇的PBC患者血清膽紅素觀測值及擬合曲線

Fig.2TheobservedvaluesandsimulationcurvesofserBilirfor12randomlyselectedPBCpatients

圖2直觀地顯示了治療對降低患者的血清膽紅素只有一開始產(chǎn)生了較好的效果，持續(xù)大概兩個月的時間之后開始反彈，直到3.5 a的時候血清膽紅素的增速開始放緩，但仍保持增長趨勢，說明治療方法需要改進．

4　結(jié)　語

本研究首先對分段LME模型中合適的變點個數(shù)進行選擇，使得聯(lián)合模型對縱向及生存數(shù)據(jù)具有最優(yōu)擬合度，進而建立了共享協(xié)變量和隨機效應的含有多個變點的縱向數(shù)據(jù)LME模型和生存數(shù)據(jù)AFT模型的聯(lián)合模型．在此基礎上，本研究得到了變點在時間軸上的位置估計值，并通過Gauss-Hermite近似解決了聯(lián)合似然函數(shù)中的復雜積分以得到模型中其他參數(shù)的估計，較好地解釋了協(xié)變量及隨機效應對縱向響應變量和生存時間的影響．最后將該模型運用到PBC數(shù)據(jù)中，研究發(fā)現(xiàn)：PBC患者的血清膽紅素只在初期治療階段有所降低，兩個月后迅速開始反彈，直到3.5 a后增速才有所放緩，說明治療方法仍需改進．

[1]Tsiatis A A, Degruttola V, Wulfsohn M S. Modeling the relationship of survival to longitudinal data measured with error. Applications to survival and CD4 counts in patients with AIDS [J]. Journal of the American Statistical Association, 1995, 90(429):27-37.

[2]Wulfsohn M S, Tsiatis A A. A joint model for survival and longitudinal data measured with error[J]. Biometrics, 1997, 53(1):330-339.

[3]Hinkley D V. Inference about the change-point in a sequence of random variables[J]. Biometrika, 1970, 57(1):1-17.

[4]Smith A F M, Cook D G. Straight lines with a change-point:a Bayesian analysis of some renal transplant data [J]. Applied Statistics, 1980, 29(2):180-189.

[5]Kim H M, Lagakos S W. Assessing drug compliance using longitudinal marker data, with application to AIDS[J]. Statistics in Medicine, 1995, 13(19-20):2141-2153.

[6]Kim H J, Fay M P, Feuer E J,etal. Permutation tests for joinpoint regression with applications to cancer rates[J]. Statistics in Medicine, 2000, 19(3):335-351.

[7]Huang Y, Dagne G A, Park J G. Segmental modeling of changing immunologic response for CD4 data with skewness, missingness and dropout [J]. Journal of Applied Statistics, 2013, 40(10):2244-2258.

[8]Evans M, Swartz T. Approximating Integrals via Monte Carlo and Deterministic Methods [M]. Oxford:Oxford University Press, 2000.

[9]Murtaugh P A, Dickson E R, van Dam G M,etal. Primary biliary cirrhosis:Prediction of short-term survival based on repeated patient visits [J]. Hepatology, 1994, 20(1):126-134.

[10]Shapiro J M, Smith H, Schaffner F. Serum bilirubin:a prognostic factor in primary biliary cirrhosis[J]. Gut, 1979, 20(2):137-140.

Multiple change points identification in joint modeling of longitudinal and survival data

SHENJia-kun1,2,SONGLi-xin*2,SUNXiu-feng1,FENGBao-jun1

( 1.Faculty of Management and Economics, Dalian University of Technology, Dalian 116024, China；2.School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China )

A joint model with multiple change points identifying in longitudinal response process is proposed, which combines a linear mixed-effect (LME) model and an accelerated failure time (AFT) model with respect to shared covariates and random effects. All the parameters are estimated by the maximum likelihood function through the Gauss-Hermite approximation to deal with the intractable integrals in it. The effect of the method is elucidated through simulation studies and a real data application about primary biliary cirrhosis (PBC). It is shown that serum bilirubin level declines only at the beginning of treatment and lasts two months, then quickly rebounds and doesn′t slow down until 3.5 years later, which indicates that the treatment methods still need to be improved.

multiple change points; linear mixed-effect (LME) model; accelerated failure time (AFT) model; joint inference; maximum likelihood

1000-8608(2016)05-0539-07

2016-01-15；

2016-07-08．

國家社會科學基金資助項目(16BGL060)；國家自然科學基金資助項目(11371077).

沈佳坤(1991-)，女，博士生，E-mail:shenjiakun@mail.dlut.edu.cn；宋立新*(1966-)，男，教授，博士生導師，E-mail:lxsong@dlut.edu.cn；馮寶軍(1966-)，男，教授，博士生導師，E-mail:fbj066@sina.com.

O212

A

10.7511/dllgxb201605015