本文轉載自：張天嵩，董聖傑，楊智榮，武珊珊，田金徽，孫鳳*. 網路Meta分析研究進展系列(十三):生存資料的網路Meta分析[J]. 中國循證心血管醫學雜誌. 2021;13(6):649-652

感謝《中國循證心血管醫學雜誌》和孫鳳教授授權醫咖會轉載。

在生物醫學領域，有很多研究感興趣的結局是到某個事件發生的時間，發生的事件可以是不利結果（如死亡、腫瘤復發等）、有益結果（如懷孕、出院等）、中性結果（如停止母乳時間）“，這類資料稱為

至事件時間資料（time-to-event data）

，傳統上稱為生存資料或生存資料 （survival data）、隨訪資料（follow-up data）。

對該類資料行的統計分析稱為生存分析，主要考察每個研究物件出現某一結局所經歷的時間，廣泛應用於惡性腫瘤、慢性疾病等隨訪研究中事件分析，比如疾病的發生、復發、轉移、傷口的癒合、某種症狀的消失等。

生存資料具有獨特的特點

：至事件發生的時間趨向於偏態分佈，則採用正態似然不適合；在隨訪期內並非每例患者都出現待觀察的事件，即有截尾觀測。

在生存資料中，

結局事件和發生時間同等重要

，主要的分析方法為風險或生存函式。針對生存資料進行經典的Meta分析和網路Meta分析（network meta-analysis， NMA），與連續型資料和二分類資料相比更加棘手，既不能簡單地將時間作為連續型資料處理，也不能把事件簡單地作為二分類資料處理。

生存資料NMA，一般可分為

一步法和兩步法

兩種分析策略，一步法適用個體參與者資料（individual participant data，IPD）；兩步法適用於IPD和聚合資料（aggregate data，AD），

目前多使用兩步法

。

本文主要梳理

基於兩步法生存資料NMA的相關方法和模型

，以期為研究人員瞭解和開展相關工作提供參考。

1、生存資料型別

1.1 個體參與者資料

理想情況下，系統評價員能夠獲得納人NMA每個原始研究的IPD，收集到原始研究的起點事件和終點事件及相關人數、生存時間、影響因素（協變數）等。IPD既直接用於一步法NMA建模；也可採用適合的模型計算出單個研究的尺度引數和位置引數，然後進行NMA效應量合併。

1.2 聚合資料

如果不能獲得原始研究的IPD，則需要依賴於研究水平的聚合資料如相應的概括統計量（summary statistics）。

一般情況下，報告的生存資料效應指標多為生存率、生存期、生存風險比/危險比（Hazard ratio，HR）等

。

其中，生存率指標主要有總生存率（overall survival，OS）和無進展生存率（progression-free survival，PFS）等，有時也會報告在某個時間點患者存活比例（時點生存率）等；生存期指標主要有平均生存期、中位生存期、無病存活期、無“事”存活期（event free survival，EFS）等。

但

在實踐中，

HR是用來比較不同組別至事件時間干預效應最常用、最合適的指標，但遺憾的是有不少研究並未直接報告HR

，但可以透過其他統計量轉換或從Kaplan-Meier （K-M）曲線中提取來獲得，計算HR及其方差的主要方法如表1所示。

表1。 HR及其方差的計算方法

2、效應量計算

在系統評價和Meta分析中，生存資料最合適的效應量是HR，如果不能獲得IPD，則需要小心採用相應方法從研究中水平來獲得HR及相應統計量，本文重點討論

從概括統計量轉換或從K-M曲線中估算HR。

2.1 時序檢驗

Tierney等［1］提供了基於研究的報告資訊如何從單個研究中計算HR及相應統計量，並給出大量公式、詳細解釋、實用例子和基於 Excel開發的計算工具，非常有參考意義。

例如，假設能從時序檢驗（log-rank test）獲得實驗組相對於治療組的lnHR可透過（O-E）/V來估計，標準誤為1/√V。式中O為每個研究中實驗組事件發生數，E為實驗組期望發生事件數，（O-E）為log-rank統計量，V為log-rank統計量的方差。

一項新輔助化療治療口咽癌研究的時序檢驗結果如表2所示，則根據公式很容易的獲得HR的對數及相應方差或標準誤，分別為：

2.2 Cox比例風險迴歸模型

在納人Meta分析的隨機對照或佇列研究中，常用Cox比例風險迴歸模型（Cox proportional hazards regression models，簡稱Cox迴歸模型）來計算HR，一般會直接給出校正後的HR值及其95%可信區間。

有時，多臂研究中存在兩個及以上的干預措施，但只報告相對於同一參照（如安慰劑）的相對干預效應，如一項新輔助化療治療壺腹周圍癌術後患者的臨床研究，將性別、吸菸狀態、干預措施等作為調節因素納入Cox迴歸模型進行分析，選取干預措施相應結果如表3所示，共分為三個治療組，分別為觀察組、氟脲嘧啶（fluorouracil）聯合亞葉酸（folinic acid）組、吉西他濱（gemcitabine）組。

針對表3中的資料，可以按下列方法獲得兩兩比較的效應量（如HR對數尺度）及其標準誤。首先，根據Tierney法［1］分別獲得B vs。 A和C vs。 A風險比HR對數尺度的相應標準誤：

其次，根據Woods法獲得C vs。 B風險比的對數尺度及相應近似標準誤：

2.3 Kaplan-Meier曲線

K-M曲線法是一種用於評估風險函式的非引數方法，不需要對生存時間的分佈進行特殊指定，並圖示化結果顯示資料的風險函式估計。如果原始研究採用K-M曲線報告了OS、PFS等生存結局，則可從中獲得隨時間變化的生存比例等資料。

如Guyot等建立了一種演算法，可從已發表的K-M曲線中重建“偽（pseudo）”IPD，從而獲得相應的HR等效應量，可用於進行經典的Meta分析或NMA。國內已有學者［2］介紹了採用圖形資料提取軟體（如Engauge Digitizer軟體）處理K-M曲線、提取資料，並計算HR及其方差的具體操作過程，均可供參考，本文不再贅述。

2.4 二分類資料或計數資料轉換

有的研究可能報告了時點生存率，但未報告獲HR等其他資料，如果沒有截尾值，可透過相對危險度（RR）或者比值比（OR）替代計算HR，但RR或者OR並未考慮時間因素，與HR相比丟失了一些重要的資訊，故這種處理方法不作為常規選擇。

新近Watkins提出了一種合併二分類資料與HR的簡單方法，與Woods法異曲同工，其中也提供了二分類資料（計數或比例）轉換為HR的方法：為簡單起見，假設每個兩臂研究中的臂為k（k=1，2），每個臂的總人數為nk，時間為t*上觀察到事件發生人數為rk（t*），或比例pk（t*）=rk（t*）/nk（t*），則HR的對數尺度及相應標準誤分別為：

3、生存資料NMA模型

3.1 圖論法

3.1.1 基本模型

由Rücker等基於電網理論提出。假設在NMA中有n個不同干預措施（點），存在m個比較（線），顯然，如果納入的均為兩臂研究，則m為納入分析研究的個數，假設干預措施（點1，。。。，n）和研究（邊1，。。。，m）按固定順序任意編號。

令向量y=（y1，L，ym）T和w=（w1，。。。L，wm）T分別表示觀測效應量及其抽樣方差倒數，並定義對角線m m矩陣W=diag（w），含有研究（邊）在其對角線線上的權重wi，令Wy為經方差倒加權觀測效應（wiyi）i的點乘向量。

為了給定比較，首先要定義網路結構，Rücker引入幾個非常重要的概念，一是邊點發生矩陣（edge-vertex incidence），是m n設計矩陣，行表示研究（邊），在同一行中，某列用“1”表示干預措施（點）所在研究（邊）的起點，“-1”表示終點；

二是拉普拉斯矩陣（Laplacian matrix），是n n矩陣，定義為L=BTWB；三是拉普拉斯矩陣的摩爾-彭若斯廣義逆（Moore-Penrosepseudoinverse）矩陣L+=（L-）-1+J/n，其中，J為含1的n n矩陣。透過這些矩陣，可以計算直接比較和間接比較的方差，進而估計Q統計量、成對比較的效應量等。

3.1.2 軟體實現

如果能夠獲得每個研究中配對比較干預措施的HR及其方差或標準誤，該模型可透過

R軟體的netmeta包

擬合來實現。netmeta包是基於頻率框下NMA的R軟體擴充套件包，由Rücker等編寫，目前版本為1。3-0，於2021-1-18釋出，需要在3。5。0及以上版本的R軟體執行［3］。

該包的資料輸人格式、計算程式設計極其簡單易行，感興趣的讀者可閱讀其自帶幫助檔案或相關文獻介紹。

3.2 引數化生存曲線法

3.2.1 基本模型

引數化生存曲線（parametricsurvival curves）法由Ouwens等與Jansen提出，是用引數化生存曲線的引數替代HR建模。假設ln（hjkt）表示第j個研究第k個干預措施在時間點的風險率（Hazard rate）；向量

為研究特定效應，用來表示基線干預（baseline treatment）b的尺度引數v。和形狀引數θ；向量

表示第j個研究中干預措施威布林生存曲線在尺度引數和形狀引數θ上的差異，可透過合併由參照干預 （reference treatment）表達的效應估計來刻畫；∑表示協方差矩陣，用來表示異質性；σ1和σ2分別表示在尺度引數v和形狀引數θ上的異質性，ρ表示兩引數的相關性，則兩臂研究的NMA隨機效應模型為：

該模型雖然以威布林分佈為例，但從原理上而言，也適用於其它引數化生存分佈，如龔帕茲分佈（Gompertz distribution）、對數邏輯斯諦分佈（Log-logistic distribution）、對數正態分佈（log-normal distribution）等。

3.2.2 軟體實現

Ouwens等提出的模型在貝葉斯分析框架下由WinBUGS/OpenBUGS、JAGS、R等軟體來擬合。實際上，

文獻中[4]給出了模型的WinBUCS程式程式碼、資料輸入格式、例項資料分析等

，可供參考；研究者們也可以根據分析目的、資料型別、軟體功能和自己對軟體的掌握程度等方面來綜合考慮選擇其他合適的軟體。

3.3 分段多項式模型

3.3.1 基本模型

分段多項式（Fractional polynomials，FP）是一族靈活引數模型，用於描述變數間的相關性，模型假定隨時間變化干預效應不呈線性變化。Jansen等提出將FP模型應於NMA中，如將第s個研究中第j個干預措施M階次FP結構的測量平均值定義為：

其中us表示研究平均值，δj表示干預措施j相對於參照干預的研究和時間特定效應，協方差矩陣∑表示干預效應引數δmsj的研究間異質性。

令t=ln（t）；模型中感興趣的協變數（時間）冪次p從S={-2，-1，-0。5，0，0。5，1，2，3}中選取；階次為m=1，L，在醫學應用中，一階和二階轉換最常見，故M選1和2即可。

3.3.2 軟體實現

Jansen等提出的模型可以在貝葉斯分析框架下，由WinBUGS/OpenBUGS、JAGS、R等軟體來擬合。

文獻中[5]給出了模型的WinBUCS程式程式碼、資料輸入格式、例項資料分析等

，可供參考。

3.4 其他模型

針對生存資料，有研究者探索一些適用於結構更為複雜的資料的統計模型，如Woods等基於研究中每個臂的事件風險對數尺度及其方差或標準誤，建立了合併計數和HR統計量的NMA模型；Saramago等建立了合併IPD和AD不同資料型別的NMA模型；Cope等認為 Ouwens和Jansen提出的模型存在一定的不足，並對其進行改良，提出基於生存函式引數的兩步法即多元NMA策略。

在實踐中，研究者可以參考這些模型，並可根據實際情況靈活選用文獻中所提供的分析程式碼。如只獲得了HR統計量，則可從Woods模型程式碼中選擇合HR統計量的部分用於分析自己的資料。

4、結語

當前，尚無明確的指南或推薦意見用於指導生存資料進行NMA。理想情況下，直接利用IPD建模進行統計分析。

如果不能獲得IPD，則針對 AD有兩種主要建模方法

，一是擬合HR資料，但需要注意的是，合併分別由時序檢驗和Cox迴歸得到的HR統計量時要小心，因為後者一般是得到了校正；一是擬合基於K-M曲線重建的偽IPD。

建議在實踐中要報告這兩種分析方法的結果

：如NMA網路中含有足夠的K-M曲線資料，推薦基於分段多項式模型展示結果；如NMA網路中沒有足夠的K-M曲線資料，則推薦展示基於HR的NMA分析結果。

參考文獻：

1。 Trials。 2007；8（1）：16。

2。 中國循證心血管醫學雜誌。 2016；18（1）：2-6。

3。

https：//cran。r-project。org/web/packages/netmeta/

4。 Res Synth Methods。 2010；1（3-4）：258-71。

5。 BMC Med Res Methodol。 2011；11：61。

本文轉載自：

中國循證心血管醫學雜誌。 2021；13（6）：649-652。