張?zhí)灬?，董圣杰，楊智榮，武珊珊，田金徽，孫鳳

在生物醫(yī)學(xué)領(lǐng)域，有很多研究感興趣的結(jié)局是到某個(gè)事件發(fā)生的時(shí)間，發(fā)生的事件可以是不利結(jié)果（如死亡、腫瘤復(fù)發(fā)等）、有益結(jié)果（如懷孕、出院等）、中性結(jié)果（如停止母乳時(shí)間）[1]，這類數(shù)據(jù)稱為至事件時(shí)間數(shù)據(jù)（time-toevent data），傳統(tǒng)上稱為生存數(shù)據(jù)或生存資料（survival data）、隨訪資料（follow-up data），對(duì)該類數(shù)據(jù)行的統(tǒng)計(jì)分析稱為生存分析，主要考察每個(gè)研究對(duì)象出現(xiàn)某一結(jié)局所經(jīng)歷的時(shí)間，廣泛應(yīng)用于惡性腫瘤、慢性疾病等隨訪研究中事件分析，比如疾病的發(fā)生、復(fù)發(fā)、轉(zhuǎn)移、傷口的愈合、某種癥狀的消失等[2-4]。生存數(shù)據(jù)具有獨(dú)特的特點(diǎn)：①至事件發(fā)生的時(shí)間趨向于偏態(tài)分布，則采用正態(tài)似然不適合；②在隨訪期內(nèi)并非每例患者都出現(xiàn)待觀察的事件，即有截尾觀測(cè)。

在生存數(shù)據(jù)中，結(jié)局事件和發(fā)生時(shí)間同等重要[5]，主要的分析方法為風(fēng)險(xiǎn)或生存函數(shù)[6]。針對(duì)生存數(shù)據(jù)進(jìn)行經(jīng)典的Meta分析和網(wǎng)絡(luò)Meta分析（network meta-analysis，NMA），與連續(xù)型數(shù)據(jù)和二分類數(shù)據(jù)相比更加棘手，既不能簡(jiǎn)單地將時(shí)間作為連續(xù)型數(shù)據(jù)處理，也不能把事件簡(jiǎn)單地作為二分類數(shù)據(jù)處理。生存數(shù)據(jù)NMA，一般可分為一步法和兩步法兩種分析策略，一步法適用個(gè)體參與者數(shù)據(jù)（individual participant data，IPD）[7,8]；兩步法適用于IPD和聚合數(shù)據(jù)（aggregate data，AD），目前多使用兩步法。本文主要梳理基于兩步法生存數(shù)據(jù)NMA的相關(guān)方法和模型，以期為研究人員了解和開(kāi)展相關(guān)工作提供參考。

1 生存數(shù)據(jù)類型

1.1 個(gè)體參與者數(shù)據(jù)理想情況下，系統(tǒng)評(píng)價(jià)員能夠獲得納入NMA每個(gè)原始研究的IPD，收集到原始研究的起點(diǎn)事件和終點(diǎn)事件及相關(guān)人數(shù)、生存時(shí)間、影響因素（協(xié)變量）等。IPD既直接用于一步法NMA建模；也可采用適合的模型計(jì)算出單個(gè)研究的尺度參數(shù)和位置參數(shù)，然后進(jìn)行NMA效應(yīng)量合并。

1.2 聚合數(shù)據(jù)如果不能獲得原始研究的IPD，則需要依賴于研究水平的聚合數(shù)據(jù)如相應(yīng)的概括統(tǒng)計(jì)量（summary statistics）。一般情況下，報(bào)告的生存數(shù)據(jù)效應(yīng)指標(biāo)多為生存率、生存期、生存風(fēng)險(xiǎn)比/危險(xiǎn)比（Hazard ratio，HR）等[2]。其中，生存率指標(biāo)主要有總生存率（overall survival，OS）和無(wú)進(jìn)展生存率（progression-free survival，PFS）等，有時(shí)也會(huì)報(bào)告在某個(gè)時(shí)間點(diǎn)患者存活比例（時(shí)點(diǎn)生存率）等；生存期指標(biāo)主要有平均生存期、中位生存期、無(wú)病存活期、無(wú)“事”存活期（event free survival，EFS）等；但在實(shí)踐中，HR是用來(lái)比較不同組別至事件時(shí)間干預(yù)效應(yīng)最常用、最合適的指標(biāo)[9]，但遺憾的是有不少研究并未直接報(bào)告HR，但可以通過(guò)其他統(tǒng)計(jì)量轉(zhuǎn)換或從Kaplan-Meier（K-M）曲線中提取來(lái)獲得，計(jì)算HR及其方差的主要方法如表1所示[9]。

表1 HR及其方差的計(jì)算方法

2 效應(yīng)量計(jì)算

在系統(tǒng)評(píng)價(jià)和Meta分析中，生存數(shù)據(jù)最合適的效應(yīng)量是HR，如果不能獲得IPD，則需要小心采用相應(yīng)方法從研究中水平來(lái)獲得HR及相應(yīng)統(tǒng)計(jì)量[10]，本文重點(diǎn)討論從概括統(tǒng)計(jì)量轉(zhuǎn)換或從K-M曲線中估算HR。

2.1 時(shí)序檢驗(yàn)Tierney等[10]提供了基于研究的報(bào)告信息如何從單個(gè)研究中計(jì)算HR及相應(yīng)統(tǒng)計(jì)量，并給出大量公式、詳細(xì)解釋、實(shí)用例子和基于Excel開(kāi)發(fā)的計(jì)算工具，非常有參考意義。例如，假設(shè)能從時(shí)序檢驗(yàn)（log-rank test）獲得實(shí)驗(yàn)組相對(duì)于治療組的lnHR可通過(guò)(O-E)/V來(lái)估計(jì)，標(biāo)準(zhǔn)誤為數(shù)，E為實(shí)驗(yàn)組期望發(fā)生事件數(shù)，(O-E)為logrank統(tǒng)計(jì)量，V為log-rank統(tǒng)計(jì)量的方差。一項(xiàng)新輔助化療治療口咽癌研究[11]的時(shí)序檢驗(yàn)結(jié)果如表2所示，則根據(jù)公式很容易的獲得HR的對(duì)數(shù)及相應(yīng)方差或標(biāo)準(zhǔn)誤，分別為:

表2 新輔助化療干預(yù)口咽癌的時(shí)序檢驗(yàn)結(jié)果

2.2 Cox比例風(fēng)險(xiǎn)回歸模型在納入Meta分析的隨機(jī)對(duì)照或隊(duì)列研究中，常用Cox比例風(fēng)險(xiǎn)回歸模型（Cox proportional hazards regression models，簡(jiǎn)稱Cox回歸模型）來(lái)計(jì)算HR[5]，一般會(huì)直接給出校正后的HR值及其95%可信區(qū)間。有時(shí)，多臂研究中存在兩個(gè)及以上的干預(yù)措施，但只報(bào)告相對(duì)于同一參照（如安慰劑）的相對(duì)干預(yù)效應(yīng)，如一項(xiàng)新輔助化療治療壺腹周圍癌術(shù)后患者的臨床研究[12]，將性別、吸煙狀態(tài)、干預(yù)措施等作為調(diào)節(jié)因素納入Cox回歸模型進(jìn)行分析，選取干預(yù)措施相應(yīng)結(jié)果如表3所示，共分為三個(gè)治療組，分別為觀察組、氟脲嘧啶（fluorouracil）聯(lián)合亞葉酸（folinic acid）組、吉西他濱（gemcitabine）組。

表3 新輔助化治療壺腹周圍癌術(shù)后的Cox回歸模型分析結(jié)果

針對(duì)表3中的數(shù)據(jù)，可以按下列方法獲得兩兩比較的效應(yīng)量（如HR對(duì)數(shù)尺度）及其標(biāo)準(zhǔn)誤。首先，根據(jù)Tierney法[10]分別獲得Bvs. A和Cvs. A風(fēng)險(xiǎn)比HR對(duì)數(shù)尺度的相應(yīng)標(biāo)準(zhǔn)誤：

其次，根據(jù)Woods法[13]獲得Cvs. B風(fēng)險(xiǎn)比的對(duì)數(shù)尺度及相應(yīng)近似標(biāo)準(zhǔn)誤:

2.3 Kaplan-Meier曲線K-M曲線法是一種用于評(píng)估風(fēng)險(xiǎn)函數(shù)的非參數(shù)方法，不需要對(duì)生存時(shí)間的分布進(jìn)行特殊指定[6]，并圖示化結(jié)果顯示數(shù)據(jù)的風(fēng)險(xiǎn)函數(shù)估計(jì)。如果原始研究采用K-M曲線報(bào)告了OS、PFS等生存結(jié)局，則可從中獲得隨時(shí)間變化的生存比例等數(shù)據(jù)。如Guyot等[14]建立了一種算法，可從已發(fā)表的K-M曲線中重建“偽（pseudo）”IPD，從而獲得相應(yīng)的HR等效應(yīng)量，可用于進(jìn)行經(jīng)典的Meta分析或NMA。國(guó)內(nèi)已有學(xué)者[15-17]介紹了采用圖形數(shù)據(jù)提取軟件（如Engauge Digitizer軟件）處理K-M曲線、提取數(shù)據(jù)，并計(jì)算HR及其方差的具體操作過(guò)程，均可供參考，本文不再贅述。

2.4 二分類數(shù)據(jù)或計(jì)數(shù)數(shù)據(jù)轉(zhuǎn)換有的研究可能報(bào)告了時(shí)點(diǎn)生存率，但未報(bào)告獲HR等其他數(shù)據(jù)，如果沒(méi)有截尾值，可通過(guò)相對(duì)危險(xiǎn)度（RR）或者比值比（OR）替代計(jì)算HR，但RR或者OR并未考慮時(shí)間因素，與HR相比丟失了一些重要的信息，故這種處理方法不作為常規(guī)選擇。新近Watkins[18]提出了一種合并二分類數(shù)據(jù)與HR的簡(jiǎn)單方法，與Woods[13]法異曲同工，其中也提供了二分類數(shù)據(jù)（計(jì)數(shù)或比例）轉(zhuǎn)換為HR的方法：為簡(jiǎn)單起見(jiàn)，假設(shè)每個(gè)兩臂研究中的臂為k(k=1,2)，每個(gè)臂的總?cè)藬?shù)為nk，時(shí)間為t*上觀察到事件發(fā)生人數(shù)為rk(t*)，或比例pk(t*)=rk(t*)/nk(t*)，則HR的對(duì)數(shù)尺度及相應(yīng)標(biāo)準(zhǔn)誤分別為：。式中O為每個(gè)研究中實(shí)驗(yàn)組事件發(fā)生

3 生存數(shù)據(jù)NMA模型

3.1 圖論法

3.1.1 基本模型由Rücker等[19-21]基于電網(wǎng)理論提出。假設(shè)在NMA中有n個(gè)不同干預(yù)措施（點(diǎn)），存在m個(gè)比較（線），顯然，如果納入的均為兩臂研究，則m為納入分析研究的個(gè)數(shù)，假設(shè)干預(yù)措施（點(diǎn)1，…，n）和研究（邊1，…，m）按固定順序任意編號(hào)。令向量y=(y1,L,ym)T和w=(w1,L,wm)T分別表示觀測(cè)效應(yīng)量及其抽樣方差倒數(shù)，并定義對(duì)角線m×m矩陣W=diag(w)，含有研究（邊）在其對(duì)角線線上的權(quán)重wi，令Wy為經(jīng)方差倒加權(quán)觀測(cè)效應(yīng)(wiyi)i的點(diǎn)乘向量。為了給定比較，首先要定義網(wǎng)絡(luò)結(jié)構(gòu)，Rücker引入幾個(gè)非常重要的概念，一是邊點(diǎn)發(fā)生矩陣（edgevertex incidence），是m×n設(shè)計(jì)矩陣，行表示研究（邊），在同一行中，某列用“1”表示干預(yù)措施（點(diǎn)）所在研究（邊）的起點(diǎn)，“-1”表示終點(diǎn)；二是拉普拉斯矩陣（Laplacian matrix），是n×n矩陣，定義為L(zhǎng)=BTWB；三是拉普拉斯矩陣的摩爾-彭若斯廣義逆（Moore-Penrose pseudoinverse）矩陣L+=(L-J)-1+J/n，其中，J為含1的n×n矩陣。通過(guò)這些矩陣，可以計(jì)算直接比較和間接比較的方差，進(jìn)而估計(jì)Q統(tǒng)計(jì)量、成對(duì)比較的效應(yīng)量等。

3.1.2 軟件實(shí)現(xiàn)如果能夠獲得每個(gè)研究中配對(duì)比較干預(yù)措施的HR及其方差或標(biāo)準(zhǔn)誤，該模型可通過(guò)R軟件的netmeta包擬合來(lái)實(shí)現(xiàn)。netmeta包是基于頻率框下NMA的R軟件擴(kuò)展包，由Rücker等編寫(xiě)，目前版本為1.3-0，于2021-1-18發(fā)布，需要在3.5.0及以上版本的R軟件運(yùn)行[21]。該包的數(shù)據(jù)輸入格式、計(jì)算編程極其簡(jiǎn)單易行，感興趣的讀者可閱讀其自帶幫助文件[21]或相關(guān)文獻(xiàn)介紹[19,22,23]。

3.2 參數(shù)化生存曲線法

3.2.1 基本模型參數(shù)化生存曲線（parametric survival curves）法由Ouwens等[24]與Jansen[25]提出，是用參數(shù)化生存曲線的參數(shù)替代HR建模。假設(shè)ln(hjkt)表示第j個(gè)研究第k個(gè)干預(yù)措施在時(shí)間點(diǎn)t的風(fēng)險(xiǎn)率（Hazard rate）；向量為研究特定效應(yīng)，用來(lái)表示基線干預(yù)（baseline treatment）b的尺度參數(shù)v和形狀參數(shù)θ；向量表示第j個(gè)研究中干預(yù)措施 威布爾生存曲線在尺度參數(shù) 和形狀參數(shù)θ上的差異，可通過(guò)合并由參照干預(yù)（reference treatment）表達(dá)的效應(yīng)估計(jì)來(lái)刻畫(huà)； ∑表示協(xié)方差矩陣，用來(lái)表示異質(zhì)性；σ1和σ2分別表示在尺度參數(shù)v和形狀參數(shù)θ上的異質(zhì)性， ρ表示兩參數(shù)的相關(guān)性，則兩臂研究的NMA隨機(jī)效應(yīng)模型為[24]：

該模型雖然以威布爾分布為例，但從原理上而言，也適用于其它參數(shù)化生存分布[24]，如龔帕茲分布（Gompertz distribution）、對(duì)數(shù)邏輯斯諦分布（Log-logistic distribution）、對(duì)數(shù)正態(tài)分布（log-normal distribution）等。

3.2.2 軟件實(shí)現(xiàn)Ouwens等[24]提出的模型在貝葉斯分析框架下由WinBUGS/OpenBUGS、JAGS、R等軟件來(lái)擬合。實(shí)際上，文獻(xiàn)[24]中給出了模型的WinBUGS程序代碼、數(shù)據(jù)輸入格式、實(shí)例數(shù)據(jù)分析等，可供參考；研究者們也可以根據(jù)分析目的、數(shù)據(jù)類型、軟件功能和自己對(duì)軟件的掌握程度等方面來(lái)綜合考慮選擇其他合適的軟件。

3.3 分段多項(xiàng)式模型

3.3.1 基本模型分段多項(xiàng)式（Fractional polynomials，F(xiàn)P）是一族靈活參數(shù)模型，用于描述變量間的相關(guān)性，模型假定隨時(shí)間變化干預(yù)效應(yīng)不呈線性變化[26]。Jansen等[25,27]提出將FP模型應(yīng)于NMA中，如將第s個(gè)研究中第j個(gè)干預(yù)措施M階次FP結(jié)構(gòu)的測(cè)量平均值定義為：

其中μs表示研究平均值，δj表示干預(yù)措施j相對(duì)于參照干預(yù)的研究和時(shí)間特定效應(yīng)，協(xié)方差矩陣∑表示干預(yù)效應(yīng)參數(shù)δmsj的研究間異質(zhì)性。令t0: =ln(t)；模型中感興趣的協(xié)變量（時(shí)間）冪次p從S={-2,-1,-0.5,0,0.5,1,2,3}中選??；階次為m=1,L,M，在醫(yī)學(xué)應(yīng)用中，一階和二階轉(zhuǎn)換最常見(jiàn)，故M選1和2即可。

3.3.2 軟件實(shí)現(xiàn)Jansen等[25]提出的模型可以在貝葉斯分析框架下，由WinBUGS/OpenBUGS、JAGS、R等軟件來(lái)擬合。文獻(xiàn)[25]中給出了模型的WinBUGS程序代碼、數(shù)據(jù)輸入格式、實(shí)例數(shù)據(jù)分析等，可供參考。

3.4 其他模型針對(duì)生存數(shù)據(jù)，有研究者探索一些適用于結(jié)構(gòu)更為復(fù)雜的數(shù)據(jù)的統(tǒng)計(jì)模型，如Woods等[13]基于研究中每個(gè)臂的事件風(fēng)險(xiǎn)對(duì)數(shù)尺度及其方差或標(biāo)準(zhǔn)誤，建立了合并計(jì)數(shù)和HR統(tǒng)計(jì)量的NMA模型；Saramago等[28]建立了合并IPD和AD不同數(shù)據(jù)類型的NMA模型；Cope等[29]認(rèn)為Ouwens[24]和Jansen[25,27]提出的模型存在一定的不足，并對(duì)其進(jìn)行改良，提出基于生存函數(shù)參數(shù)的兩步法即多元NMA策略。在實(shí)踐中，研究者可以參考這些模型，并可根據(jù)實(shí)際情況靈活選用文獻(xiàn)中所提供的分析代碼。如只獲得了HR統(tǒng)計(jì)量，則可從Woods[13]模型代碼中選擇合HR統(tǒng)計(jì)量的部分用于分析自己的數(shù)據(jù)。

4 結(jié)語(yǔ)

當(dāng)前，尚無(wú)明確的指南或推薦意見(jiàn)用于指導(dǎo)生存數(shù)據(jù)進(jìn)行NMA。理想情況下，直接利用IPD建模進(jìn)行統(tǒng)計(jì)分析。如果不能獲得IPD，則針對(duì)AD有兩種主要建模方法，一是擬合HR數(shù)據(jù)，但需要注意的是，合并分別由時(shí)序檢驗(yàn)和Cox回歸得到的HR統(tǒng)計(jì)量時(shí)要小心，因?yàn)楹笳咭话闶堑玫搅诵Ｕ?；一是擬合基于K-M曲線重建的偽IPD。建議在實(shí)踐中要報(bào)告這兩種分析方法的結(jié)果：如NMA網(wǎng)絡(luò)中含有足夠的K-M曲線數(shù)據(jù)，推薦基于分段多項(xiàng)式模型展示結(jié)果；如NMA網(wǎng)絡(luò)中沒(méi)有足夠的K-M曲線數(shù)據(jù)，則推薦展示基于HR的NMA分析結(jié)果。