縱向數據與生存數據聯合模型中多變點識別問題

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

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

?

縱向數據與生存數據聯合模型中多變點識別問題

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

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

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

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

0　引　言

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

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

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

1　模型和方法

1.1聯合模型

考慮一個樣本量為N的樣本，個體間相互獨立．令Yij為個體i在時間tij的響應變量，j=1,…,ni且i=1,…,N．本文建立關于縱向數據的含有k個變點的線性混合效應(LME)模型(1)，其中k為任意正整數，并將k個變點在時間軸上的位置作為未知參數：τ=(τ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的縱向數據模型(1)有關．為體現這種關聯性，研究中假設生存數據模型與縱向數據模型(1)共享協變量向量Zi和隨機效應向量ui．特別的，建立混合效應加速失效時間(AFT)模型(2)：

(2)

1.2變點識別

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

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

LR(k)=

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

(3)

θ^(k)0

θ^(k)1

1.3聯合似然推斷

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

2　數值模擬

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節中的右刪失假設．

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

(7)

下面通過一系列數值模擬研究來驗證上一章中提出的聯合模型及估計方法的有效性．本節的數值模擬基于共享的Z和u的聯合模型(6)和(7)．假設樣本容量為N=150或N=300，個體間相互獨立，每個個體的重復觀測次數為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個時間點得到次序統計量ti，并建立生存時間的右刪失示性函數δi=I(Ti≤ti6)，其中ti6是上述抽取均勻分布時間點的最大次序統計量，假設變點τ=(13)，且參數γ=(1.20.50.2)T；當ni=10時，同理，本文從均勻分布U(0，10)無放回隨機抽取的10個時間點得到次序統計量ti，并建立生存時間的右刪失示性函數δi=I(Ti≤ti10)，其中ti10是上述抽取均勻分布時間點的最大次序統計量，假設變點τ=(48)，參數γ=(1．80.50.2)T．同時，假設Zi服從二項分布B(1，0.5)．

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

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

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

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

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

2.2數值結果

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

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

表2列示了4種數值模擬情形下聯合模型的參數估計優度相關結果．

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

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

由圖1中的曲線可以明顯看出重復測量次數ni由6(圖(a)和(c))到10(圖(b)和(d))，估計曲線與真實曲線之間的差異減小；樣本量N由150(圖(a)和(b))到300(圖(c)和(d))，估計曲線與真實曲線之間的差異明顯減小．

表1　兩個變點位置的數值模擬結果

表2　聯合模型數值模擬的參數估計結果

　　(a)方案(a)

(b) 方案(b)

(c) 方案(c)

(d) 方案(d)

圖14種方案模擬研究曲線

Fig.1The simulation curves of four strategies

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

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

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

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

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

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

表3　PBC聯合模型似然比檢驗結果

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

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

(8)

(9)

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

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

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

表4　PBC聯合模型參數估計

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

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

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

Fig.2TheobservedvaluesandsimulationcurvesofserBilirfor12randomlyselectedPBCpatients

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

4　結　語

本研究首先對分段LME模型中合適的變點個數進行選擇，使得聯合模型對縱向及生存數據具有最優擬合度，進而建立了共享協變量和隨機效應的含有多個變點的縱向數據LME模型和生存數據AFT模型的聯合模型．在此基礎上，本研究得到了變點在時間軸上的位置估計值，并通過Gauss-Hermite近似解決了聯合似然函數中的復雜積分以得到模型中其他參數的估計，較好地解釋了協變量及隨機效應對縱向響應變量和生存時間的影響．最后將該模型運用到PBC數據中，研究發現：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