基于縱向和生存時間結局的聯合模型方法及應用＊

李淞淋，魏 戌，易丹輝

（1.農業部信息中心 北京 100125；2.中國人民大學應用統計科學研究中心 北京 100872；3.中國中醫科學院望京醫院 北京 100102）

事物經常是多維度動態發展，這就要求研究者采用全面的、動態的觀點來看待問題，即在兼顧多指標協同關系的情況下，進行動態綜合評估[1,2]。而當評價指標中既有時間資料，又有縱向指標，就需要構建聯合模型進行分析。例如在AIDS研究既關注每立方毫米血漿中CD4細胞的數量，也記錄分析患病時間[3,4]。聯合模型通過構建縱向評價指標和時間資料的聯合分布函數，在考慮指標間關系的基礎上，采用極大似然估計方法進行估計，既可以實現對兩類評價指標的聯合分析，也可以對指標之間相互關系的強度和方向進行衡量[2,5,6]。

值得注意的是，傳統的聯合模型雖然可以綜合分析時間資料和縱向結局指標，但在研究事件發生時間時采用了一個重要卻較為苛刻的前提假設——所有的研究單位都必須“死亡”或發生某一被感興趣的特定事件。也就是說，如果觀測時間可以無限延長，則結局事件必然發生，即。但此假設在某些實際研究中并不適用。首先，大多數研究都有一定的時間限制，t→∞的假設并不現實；其次，某些特定事件不是必然的。例如金融危機爆發時，一些企業不會破產或倒閉[7,8]。已有研究發現，當某件事情的發生率低于10%或5%時，)=1這一前提假設將不再合適；同時，從技術層面看，當發生率低于10%或5%時，繼續使用傳統生存分析方法也無法得到穩健的參數估計結果。

為此，本文在比較分析傳統聯合模型的適用性與優缺點基礎上，探索突破傳統聯合模型的局限性，嘗試對估計方法和統計計算方法進行改進，并應用于缺血性中風早期干預治療的實際案例。

1 傳統聯合模型討論

1.1 適用條件

傳統聯合分析模型意在考慮指標間潛在關系的基礎上，對縱向指標和事件發生時間進行聯合分析，尤其適用于如下3種情況：①關注結局事件的發生風險，但希望將縱向評價指標作為協變量引入生存分析模型來解釋縱向評價指標對結局事件的影響效果；②關注縱向指標走勢，但由于結局事件的發生導致受試患者無法繼續被調查，從而致使結局事件發生以后無法繼續觀測縱向評價指標；③關注縱向評價指標和事件結局發生時間是否有相關關系，則首先假定兩個指標相關，搭建兩者之間的相關關系，構建聯合模型，然后再在對相關關系進行檢驗，最終確定這兩類結局指標之間的關系有無與大小。

1.2 傳統聯合模型形式

傳統聯合模型（Joint Model）由兩部分構成：縱向數據部分和時間資料。假設對m個個體在時間段[0,τ)上進行觀測。第i個個體提供了一系列（可能有部分缺失）的縱向測量指標{yij,j=1…ni}，對應的觀測時間為{sij,j=1,…,ni}，此外還有某一特定結局事件發生前的“生存時間”ti（或許會存在刪失的情況）。定義Ti為觀測到的第i(i=1,…,m)個個體的失效時間（Failure Time），等于真實事件的發生時間（True Event Time）T*i和刪失時間（Censoring Time）Ci中偏小的那一個，即Ti=min(T*i,Ci)。另外，定義事件發生的標識為 δi=I(Ti*≤ Ci)，其中 I(?)是指示函數，當 Ti*≤ Ci時函數取值為1，否則取值為0。于是可將觀測到的事件時間的數據寫成集合{(Ti,δi),i=1,…,m}的形式。傳統聯合模型（Joint Models）的一般形式為：

其中，X2i為解釋變量，β2為與其對應的回歸系數；W2i(t)的形式類似于縱向數據模型中的W1i(sij)，包括個體隨機效應部分和脆弱項（Frailty）。于是聯合分析主要體現在：一是具有一些相同的解釋變量；二是W1i(sij)和W2i(t)是度量同一受試個體的隨機效應。

1.3 傳統聯合模型的優勢和不足

傳統聯合模型利用隨機效應和脆弱項來衡量個體差異并解釋兩個評價指標之間的相關關系，基于Wulfsohn和Tsiatis在1997年提出的經典假設——“給定隨機效應的情況下，縱向指標與時間資料條件獨立”，構建兩個結局指標的聯合分布函數，然后計算似然函數并求解參數估計值。但是，傳統聯合模型要求隨機效應必須服從正態分布，且事件發生率隨著時間延長而趨向于1。這些假設條件極大地限制了模型的適用環境。

2 帶混合Cure的擴展聯合模型

2.1 模型擴展的依據和思路

帶混合Cure的擴展聯合模型是在傳統聯合模型基礎上，放寬假設條件，允許隨機效應服從指數族分布即可，且不要求結局事件必須發生，從而使模型在小概率事件的環境中也可以得到科學的結果。模型擴展的依據，一是傳統聯合模型，二是混合Cure模型。

混合Cure模型最早起源于醫學研究領域。Ziegel等[9]研究指出，隨著技術手段的提高和醫療條件的改善，長期存活個體（Long-Term Survivor）即被治愈的病人數量偏多，死亡等結局事件發生率較小，成為小概率事件。為了分析這些小概率發生的事件，多位學者[10-15]提出了混合Cure模型。

混合Cure模型假設總人群（Whole Population）可以分為g個亞組，i=1,2,…,g，Zi,Xi是第i個個體的兩個協變量向量，Zi中除了有元素1外，還享有Xi中的部分元素。每個亞組在全人群中所占比重為πi(Zi)，每個亞組內人群的“死亡”模式相同，即密度函數為fi(t |Xi)，累積函數為Fi(t |Xi)。

在混合治愈模型中，最常見的情況是總人群可以分為兩個亞組：一組是根本不可能發生結局事件的免疫（或治愈）人群組；另一組則是將研究時間無限延長的話，必將發生結局事件的人群。值得注意的是，由于研究時間有限，不可能無限期延長，因此后一組人群雖然具有發生結局事件的可能性，但在有限的研究時期內，我們只能觀察到其中部分人群實際發生了結局事件。假定，在流行病調查研究中，由于個體抵抗力不同，有一些人是不會感染某種疾病的，即“非易感人群”或“治愈人群”，π1(Zi)=1-p(Zi)；另一類則是“易感人群”或“未治愈人群”，π2(Zi)=p(Zi)。于是，全體人群的分布函數FT(t |Zi,Xi)和生存函數ST(t |Zi,Xi)可分別寫成公式（3-3）和公式（3-4）。

和

其中，F(t|Xi)和S(t|Xi)分別是“未治愈人群”的分布函數和生存函數，

2.2 擴展聯合模型的形式和前提假設

2.2.1 擴展聯合模型的一般形式

假設總人群（Whole Population）中有N個個體，i=1,2,…,N，分為g個亞組，j=1,2,…,g，但這g個亞組不是預先被分好的，而是在診療或干預的過程中自然形成的。引入標識變量rij，若研究結束后第i個個體進入了第 j組，則有rij=1，否則rij=0。記第 j個亞組包含了Nj個個體，則

研究過程中，對這N個個體在時間段[0 ,τ)上進行觀測，個體i被觀測兩套協變量向量Xi={Xi和 Zi分別包含了q1和q2個協變量。對第i個個體在觀測時間序列上進行ni次觀測得到一組（可能有部分缺失）縱向測量指標得分，構成縱向結局指標的取值向量Yik={yi1,yi2,…,yik;k=1,2,…,ni}；此外，研究還記錄了個體i發生某一感興趣結局事件的時間Ti（或許會存在刪失的情況）。研究假定，個體i發生此結局事件后就將離開此研究，不再進行后續的測量。

定義Ti為觀測到的第i個個體的特定事件的發生時間，它等于事件的真實發生時間（True Event Time）和刪失時間（Censoring Time）C中偏小的那個，即iTi=min()。另外，定義事件發生的標識為 δi=I(Ti*≤ Ci)，其中 I(?)表示指示函數，當Ti*≤Ci時，δi=1，表示即研究期間內觀察到個體i發生了結局事件，否則δi=0。于是結局事件的發生情況可寫成集合{(Ti,δi);i=1,…,N}，而所有結局指標的取值效率可以表示成集合{(Yi,Ti,δi,rij);i=1,…,N;j=1,…,g}。

考慮到經過良好的、有效的措施的干預，研究總體的治愈率或免疫力會提高，因此，不良結局事件的發生率很低，limt→∞F(t)=1的假設將不再適用，需要為每個個體落入第 j個亞組的概率 pj=E(Nj/N)進行估計，

于是，“帶混合Cure的擴展聯合模型”由縱向數據子模型、人群分類子模型和時間資料子模型等三部分構成，一般形式如公式（2）。

縱向數據子模型是對研究期間內的縱向評價指標走勢進行擬合分析，其中ui(sik)為固定效應部分，W1ik(sik)為隨機效應部分，εik為測量誤差。若縱向數據子模型為線性混合效應模型，有，即此子模型為?？v向數據子模型暫時假定給定其他結局評價指標，縱向評價指標的條件分布為正態分布，將來還可以擴展到指數分布族。

亞組分類子模型和時間資料子模型類似于混合Cure模型，但都被進行了一定程度擴展。亞組分類子模型為多分類的廣義線性混合效應模型，引入了隨機效應，從而允許各亞組比例在整個研究過程中發生變化，實現對事件發生率或疾病治愈率的動態分析。

時間資料子模型是對第i個個體在研究時期內每個時間點上的結局事件發生風險進行分析，第i個個體屬于第 j組的概率為 pjk。該子模型引進了脆弱項W3i(t)=f2(b)，允許個體差異的存在，同時脆弱項的存在將縱向結局指標yi和事件發生風險h(ti|X3i;rij=1)聯系起來。例如表示：對于第 j組內的每個個體而言，個體效應b0每增加一個單位，其發生該研究事件的風險提高約倍。

需要指出的是，X1ik、X2jk和X3i是解釋變量X的子集，它們可以不相等，可以有相同的部分；Z1ik是解釋變量向量Z的子集；β1、β2、β3、b分別為 X1ik、X2jk、X3i和 Z1ik的待估系數，參數向量 γ1=(γ10,γ11,…,γ1m)和 γ2=(γ20,γ21,...,γ2m)則分別是脆弱效應 W2jk(sjk)和W3i(t)中b的待估系數向量。這是與此前關于聯合分析縱向-生存-治愈模型的研究[16]不同的地方，因此需要使用修正的半參數極大似然估計方法進行估計，使用廣義半參數似然比檢驗進行評估。

2.2.2 聯系縱向評價指標與結局事件發生概率和發生時間的方法

擴展模型搭建縱向評價指標、結局事件發生率和發生時間的聯系如下：①三個子模型有相同的協變量和參數，構建固定效應部分的關系；②縱向數據子模型的隨機效應系數b是亞組分類子模型和時間資料子模型中脆弱項的協變量，將個體效應與結局事件發生率和發生時間聯系起來；③參數向量 γ1=(γ10,γ11,…,γ1m)和 γ2=(γ20,γ21,…,γ2m)的大小可以衡量縱向評價指標與結局事件發生率和發生時間的關系強度的大小和方向。

2.2.3 聯合模型的前提假設

擴展模型需要滿足一定的前提條件：①縱向評價指標和結局事件發生率及發生時間是對同一群體在同一個研究中進行的測量；②三個結局指標（縱向結局指標、結局事件發生率、事件發生時間）之間可能存在一定的關系；③在給定隨機效應或個體差異的情況下，縱向結局指標和結局事件發生率及發生時間是條件獨立的。

3 應用分析

3.1 問題說明

本研究依托于一項“863”課題，目的是綜合評價中醫綜合治療方案和西醫治療方案針對缺血性中風早期干預治療效果。該課題得到的多家醫療單位的支持，采用中央隨機的方法對來自多家醫療中心的1 052例缺血性中風早期病患進行隨機對照試驗，隨機進入中醫組和西醫組的比例是2∶1。

根據有關醫學知識，確定了評價治療中風的藥物一年內的療效指標有NIHSS量表得分和死亡風險及發生時間。美國國立衛生研究院卒中量表（National Institute of Health Stoke Scale，NIHSS）是Thmos等為了急性腦卒中的治療研究于1989年設計的一個包含了15個項目的神經功能檢查量表。截至目前，該量表包含了每個腦動脈病變可能出現的神經系統檢查項目、精神狀態檢查項目、感覺機能、瞳孔反應和足底反射項目等。量表得分越高說明疾病狀態越嚴重。此量表使用簡便、重測信度高、內容一致性好，現已被廣泛用于腦卒中研究的定期測量記錄。此外，腦卒中研究還特別關注死亡這一客觀的結局指標，如果死亡率或死亡風險高，則意味著治療方案的效果非常差。然而，值得特別注意的是：截至目前，已經有很多關于腦卒中治療的研究發現，NIHSS量表得分與死亡風險并不是孤立存在的，兩個結局評價指標之間可能存在一定聯系。有關研究[1，17,18]發現：伴隨著NIHSS評分值的增高，死亡風險性明顯增加。

由上所述可知，目前尚無確鑿證據證明“NIHSS評分”與“因腦卒中致死”存在必然的因果關系，但是二者之間確實存在明顯的、不可被忽視的關系。同時，由于本研究只收集到1年的數據，患者的病死率較低（低于10%），為了得到更加準確的估計結果，決定采用帶混合Cure的聯合分析擴展模型進行研究。

3.2 數據說明

3.2.1 變量介紹

該研究分為住院治療期和出院預后觀察期兩部分。第一部分時長為21天，分別在0天（入院的當前或前一天）、7天、14天和21天時使用NIHSS量表進行測量，即NIHSS得分為本研究中的縱向評價指標。此外，對1 052個患者1年內的生存狀況進行跟蹤調查，及時、準確的記錄每個患者的生存狀況、死亡情況及死亡事件。下面將本研究所關注的一些變量名稱及其符號進行說明：“r_group”表示隨機分組（r_group=0表示西醫組，r_group=1表示中醫組）；“obstime”表示NIHSS量表的測量時間；“r_group*obstime”表示治療組與治療時間的交互效應，也就是考察每一個隨機分組的NIHSS量表得分隨觀測時間的變化趨勢。評價指標中，“Y”為 NIHSS量表得分，“t_death”為生存時間，“death”表示在研究過程中受試者是否死亡（death=1表示死亡，death=0表示刪失，即在1年內尚未觀測到患者死亡）。因為在建模前先對個別指標的極少量缺失值進行處理。于是，在后續研究中，所用數據中不存在數據缺失的情況。

3.2.2 基期比較

除了研究關注的結局指標外，此隨機對照試驗對入組人群的基本情況和可能影響療效評價結果的混雜因素，例如年齡、性別、病程、主要合并病和主要并發癥等都進行了隨機。經過對這些指標的統計學檢驗（表1-5），發現這些指標在兩組中的分布沒有顯著差異，即試驗的隨機效果很好。

表1 兩組間年齡比較（

表1 兩組間年齡比較（

n P值組別中西醫結合組西醫組年齡（歲）62.88±10.13 63.18±10.48 T值701 351 0.44 0.66

表2 兩組間性別比較

表3 兩組間病程比較

表3 兩組間病程比較

n T值P值組別中西醫結合組西醫組701 351病程（小時）49.48±54.62 52.77±59.65 0.892 0.373

表4 兩組間主要合并病比較

表5 兩組間主要并發病比較

圖1 兩組的NIHSS量表得分在4個時間點上變化趨勢圖

圖2 研究人群1年內的生存函數圖

繪制并分析兩組的NIHSS量表得分在0、7、14、21天這4個時間點上的發展曲線（圖1）以及研究人群在1年中的生存函數曲線（圖2）。由圖1可見，兩組的NIHSS量表得分存在差異；由圖2可以看出，1年內因腦卒中致死的人數不足總人群的10%。因此，需要采用本研究提出的帶混合Cure的聯合分析擴展模型進行分析。

3.2.3 模型構建

第一步，討論評價指標的分布特征，確定擴展聯合模型中每個子模型的形式。因為沒有充分的證據可以拒絕“給定隨機效應的情況下，縱向評價指標條件地服從正態分布”這一原假設，為此縱向數據子模型繼續采用線性混合效應模型；從死亡率不高于10%來看，這些患者經過治療，應該是出現了“治愈”和“延長生存時間”兩種結局事件，這也與已有研究的結論[1,17,18]一致。但為了進一步確定事件發生時間子模型的形式，使用圖示法數值法擬合和擬合優度評價的方法進行分析，最終確定生存時間服從Weibull分布。

第二步，變量選擇。本文使用后退法進行變量選擇，標準為參數的t檢驗、半參數似然比檢驗和模型整體擬合效果的半參數似然函數值.

第三步，構建模型。根據以上分析，確定模型形式為公式（7）。

其中，每個變量和參數的含義見公式（6）的相關介紹。對于個體i而言，治愈率為1-pi，未治愈的話則具有生存函數S(ti|r=1)。

死亡概率為

個體i對似然函數的貢獻為{Yi,Ti,δi,bi}的對數聯合似然函數為 Li(θ)。

Li(θ)=LLi(θ)+LSi(θ)+LMi(θ)。其中，

3.2.4 模型估計和檢驗

基于SAS9.2的NLMIXED模塊，使用半參數極大似然估計方法進行模型估計（表6）。

構 建 半 參 數 似 然 比 統 計 量 為 lrtn(β13,β21,β31)=17 422-9 933.67=7 488.33，根據自由度為3的卡方分布 χ2(3)，計算P值為＜0.000 1。顯著拒絕原假設 β13=0,β21=0,β31=0，即中醫組療效顯著優于西醫組。

3.2.5 模型結果分析

從表6所得每個待估參數的半參數估計結果，可知：

（1）在縱向分析子模型中，β10=7.229 2說明入組時（0天）西醫組的NIHSS平均值是7.229 2，顯著非0；β11=-0.382 2說明中醫組中患者在入組時的NIHSS得分均值為7.229 2-0.382 2=6.847 0，但檢驗得兩組人群入組時的NIHSS得分無顯著差異；β11=-0.143 2說明在21天的治療過程中，西醫組NIHSS日均下降約0.143 2，顯著非0；β13=-0.000 3說明在21天的治療過程中，中醫組NIHSS日均下降0.143 2+0.000 3=0.143 5，但檢驗結果顯示與西醫組無顯著的統計學差異。

（2）基于亞組分類子模型，計算未被完全治愈的患者所占比例為

治愈率為

β20=8.058 1說明西醫組的平均治愈率為說明中醫組的平均治愈率為，經檢驗顯著優于西醫組，約為其exp(3.168 6)≈23.77倍；γ10=1.492 0＞0且顯著非零，說明入組時NIHSS得分越高、病情越重的患者，被治愈的可能性越小，0天時NIHSS得分每提高1分，治愈率降低0.775 1；γ10=2.503 4＞0且顯著非零，說明治療過程中NIHSS分數下降越明顯，即對治療方案有明顯反應的患者，被治愈的可能性越高，治療過程中NIHSS日均降幅每增加1，治愈率提高exp(2.503 4)≈11.224 0倍，若NIHSS日均降幅每增加0.1，治愈率提高約exp(0.250 3)≈ 28.45%。

表6 擴展后聯合模型的半參數極大似然估計結果匯總表

（3）時間資料子模型中，β30=8.021 9意味著西醫組平均生存函數為exp{-tγ·exp(-8.021 9)}；β31=0.629 5不顯著非零，即中醫組平均生存函數exp{-tγ·exp(-8.651 4)}，與西醫組無顯著差異；γ20=-0.230 7＜0說明患者基線的NIHSS得分越高，即入組時腦卒中病情越重，患者的生存函數越??；同理可以發現，治療過程中患者NIHSS降低速度越慢，生存函數越小。

（4）從半參數似然比檢驗結果看，統計量lrtn(β13,β21,β31)≈7 488.3，結合（2）、（3）可知中醫組療效顯著優于西醫組，主要體現在前者對高病人治愈率和延長生存時間方面顯著優于西醫組。

4 討論與結論

4.1 擴展后的聯合模型更適用于臨床治療的綜合評價

從理論角度，傳統聯合模型聯合分析縱向和生存時間結局，并研究兩類結局指標之間的關系；擴展的聯合模型吸收傳統聯合分析方法和混合Cure模型的優點，不再繼續假設隨機效應和脆弱項服從正態分布，僅要求隨機效應的期望為零，并使用半參數估計方法和合適的算法，有效擴大適用范圍。從推廣應用角度看，擴展后聯合模型更適用于臨床治療方案的綜合評價，具體體現在：①放寬了臨床療效評價的數據類型：目前一些多靶點的臨床試驗具有多維度、多數據類型的評價指標，擴展后聯合模型不假定數據分布，非常適用于此類型試驗的分析；②提供了醫學倫理要求下，科學評價急、重癥治療方案的模型技術：傳統醫學統計模型需要依托大樣本數據，受醫學倫理要求，無法應用于癌癥、中風等急、重癥疾病臨床治療的效果評價中，而創新后模型打破了樣本數據量的限制，提高了分析精度，增強了臨床推廣應用性；③可綜合分析臨床療效和安全性，獲得準確的試驗預期：當前不少臨床試驗以AD事件發生率作為安全性指標，擴展后的聯合模型可分析提高療效所需承擔的風險，還可探索治療方案的副作用。

4.2 入組NIHSS得分與治愈率、治愈時間有顯著關系。

用擴展的聯合模型分析缺血性中風早期干預的療效，有效地實現了死亡率極低地情況下，對NIHSS量表得分與結局事件發生率和發生時間的聯合分析，研究發現：①入組時NIHSS得分越高、病情越重的患者，被治愈的可能性越低；治療過程中NIHSS分數降速越快，即對治療方案有明顯反應的患者，被治愈的可能性越高；中西醫結合治療方案的治愈率約為是西醫組的exp(3.168 6)≈23.77倍；②患者入組時NIHSS得分越高，即入組時腦卒中病情越重，患者的生存函數越小，即治愈率越低，治愈所需時間越長；治療過程中患者的NIHSS降低地越慢，生存函數越小。