網狀Meta回歸模型可用于分析重要的效應修飾因素對干預效應的影響，它能在頻率學或貝葉斯框架下實現。本文以實例介紹如何使用Stata軟件的mvmeta命令在頻率學框架下實現網狀Meta回歸，并簡單討論網狀Meta回歸的應用。
引用本文： 張天嵩. 基于頻率學框架的網狀Meta回歸實現與應用. 中國循證醫學雜志, 2015, 15(8): 988-992. doi: 10.7507/1672-2531.20150165 復制
網狀Meta分析（network meta-analyses，NMA）可同時進行直接與間接比較，將一系列不同治療方法的隨機臨床試驗數據匯總，并就給定的治療終點進行點和可信區間估計，同時對不相關性進行評估 [1-3]，但治療終點結果可能會隨不同的人群或干預特征如劑量或療程而發生變化，要探討它們之間的關系，可將Meta回歸的思想擴展到NMA中，稱為網狀Meta回歸（Network meta-regression），它可以在頻率學或貝葉斯框架下實現。本文以實例介紹如何使用Stata軟件的mvmeta命令在頻率學框架下實現網狀Meta回歸，并簡單討論網狀Meta回歸的應用。
1 資料與方法
1.1 數據來源
本文所用數據來源于治療類風濕關節炎生物制劑的Meta分析 [4]，包含Placebo、Abatacept、Adalimumab、Anakinra、Etanercept、Golimumab、Infliximab、Rituximab和Tocilizumab共9種藥物，其納入患者的指標為美國風濕病學會ACR50標準（患者關節腫脹及觸痛的個數有50%改善），原文共納入31個研究。本研究針對報告了平均病程的28個研究（含1個三臂研究，27個兩臂研究）進行分析。數據按表 1所示格式輸入Stata數據管理器中，命名為acr50.dta存儲于C盤根目錄下，其中“studyid”表示各研究，“duration”表示每個研究的平均病程，A、B、C、D、E、F、G、H、I分別表示上述9種藥物，以n*和r*分別表示每個研究中每個臂總人數和達到ACR50的人數，缺失數據（空格）以圓點表示。
表1
治療類風濕關節炎生物制劑Meta分析的具體數據
1.2 網狀Meta回歸模型
White等 [5, 6]建立了多元隨機效應Meta回歸模型（multivariate random effects meta-regression model）可以用于網狀Meta回歸分析。假設第i個研究中有相對于某一共同參照有p個比較，相應p個效應估計值yi為（1×p）向量，其相應方差-協方差矩陣為Si，則模型為yi～N(μi,Si)，μi～N(βXi,∑)。式中，μi表示研究平均效應量，Xi表示研究水平協變量，∑為研究間方差-協方差矩陣。
1.3 模型擬合
該模型擬合可由Stata軟件mvmeta命令實現 [5, 7]，該命令下載與安裝可以按參考文獻 [3, 8]的方法進行，其命令行基本操作格式為 [5]：mvmeta b V xvars [if] [in] [,options]，具體解釋和使用方法可以參閱相關文獻 [3, 5]，本文不再贅述。
為方便計算，筆者編寫了.do文件，在擬合網狀Meta回歸模型的同時，根據累積排序曲線下面積（surface under the cumulative ranking curve，SUCRA）值對各種干預藥物的療效進行排序 [9]，其具體操作步驟和主要解釋見框1。有興趣的讀者，可以在Stata軟件的菜單欄中按“Windows → Do-file Editor → New Do-file Editor”步驟操作，打開“Do-file”編輯器，將框1中代碼復制到已打開的“Do-file”編輯器中，并在編輯器的“Tools”下拉菜單中選擇“Excute(do)”執行即可。
框 1 進行網狀Meta回歸分析、干預措施療效排序的Stata代碼
**第一步：讀入數據 use C:\acr50.dta,clear **第二步：數據預處理 *建立效應量、方差－協方差矩陣 foreach trt in A B C D E F G H I{ if “`trt’”==”A” continue gen y`trt’ = log(r`trt’/(n`trt’-r`trt’)) - log(rA/(nA-rA)) gen S`trt’`trt’ = 1/r`trt’ + 1/(n`trt’-r`trt’) + 1/rA + 1/(nA-rA) foreach trt2 in A B C D E F G H I{ if “`trt2’”==”A” continue if “`trt2’”＞”`trt’” gen S`trt’`trt2’ = 1/rA + 1/(nA-rA) if !mi(r`trt’) & !mi(r`trt2’) } } format y* S* %6.2g mat P8=0.5*(I(8) +J(8,8,1) ) **第三步：擬合Meta回歸模型，對干預藥物療效進行排序 *擬合Meta回歸模型 sum duration gen x=duration-r(mean) mvmeta y S x,bscov(prop P8) pbest(max,all zero gen(probm) reps(10000) ) eform *療效排序 sucra probm*,mvmeta lab(A B C D E F G H I) noplot
為考察平均病程對干預效應的影響，本次研究采用經典的NMA（擬合一致性模型），并進行療效排序。只需要將框1中第三步中的代碼做如下修改即可。
mvmeta y S,bscov(prop P8) pbest(max,all zero gen(prob) reps(10000) ) eform
sucra prob*,mvmeta lab(A B C D E F G H I) noplot
請注意，擬合一致性模型和網狀Meta回歸模型，均基于跨研究間干預措施配對比較異質性相同的假設。
2 結果
2.1 模型擬合結果
模型擬合的迭代過程從略，只報告主要結果。擬合一致性模型、網狀Meta回歸模型結果見表 2。由表 2可見，相對于藥物A（安慰劑），其他治療類風濕關節炎的生物制劑可以提高ACR50達標率，其差異均有統計學意義（P≤ 0.001），但經過平均病程協變量（即具有8.52年平均病程時的效應量）校正后，藥物E（Etanercept）與藥物A相比，ACR50達標率差異未再出現統計學意義（P=0.806），提示病程可能與藥物E存在交互作用。
表2
網狀Meta分析及網狀Meta回歸分析結果
2.2 療效排序結果
干預藥物療效排序結果見表 3，其中SUCRA是累積排序曲線下面積的數字化表示，“PrBest”表示每個干預藥物療效最佳的概率。從表 3可知，在經典NMA結果中，藥物E的SUCRA值最大，療效最佳的可能性為89.0%，但經過平均病程協變量校正后，其療效最佳的可能性降為36.3%，也提示平均病程對藥物作用存在影響。
表3
干預藥物療效排序
3 討論
醫學研究中，干預措施呈現多樣化發展，面對多種干預措施，兩兩比較的經典Meta分析已不能為選擇最佳干預措施提供方法學支持，但NMA使得多種干預措施之間的比較成為可能 [3, 10]。NMA涉及同質性、相似性和一致性三個基本假設 [3, 11-13]，雖然這些假設各有自己的不同術語和定義 [14]，但實質上均是經典Meta分析同質性思想在NMA的擴展應用，因此本文所指的“異質性”可以解釋為異質性和/或不一致性 [15]。
NMA包括多個不同干預比較的不同試驗，因此較之于經典的Meta分析在臨床、方法學等方面更可能存在多樣性：如患者、干預措施、研究地點、隨訪時間等臨床特征，以及隨機、盲法、分配隱藏等研究設計方法等各種效應修飾因素，不僅存在于某對特定干預比較之中，也存在不同干預比較之間 [3, 13]，必然會對NMA結果的真實性產生影響。
Meta回歸常被用于探討Meta分析中潛在的效應修飾因素（稱為協變量）對干預效應量的影響，有眾多的Meta回歸模型可用于NMA中 [15]，實現方法可以基于頻率學或貝葉斯框架。貝葉斯框架，如分層NMA模型納入協變量可以由WinBUGS軟件擬合，但可能因選擇不同先驗需要進行敏感性分析，或易出現蒙特卡洛誤差；而頻率學方法，如White提出的多元Meta回歸方法則可避免上述不足，且具有計算速度快的優點 [16]，可使用Stata軟件的mvmeta命令輕松實現，本文即采用此方法。
協變量可以根據試驗的臨床特征，或根據臨床設計的偏倚風險等來選擇。協變量數據類型可分二分類或連續型，但連續型協變量需要采用中心化協變量值 [15]，也就是對其進行測量中心化（measurement centering），測量中心化方法分為總均數中心化和組均數中心化兩種。針對所引用的數據實例，選擇納入每個研究的平均病程這一試驗臨床特征為協變量，采用總均數中心化方法，即每個研究的平均病程減去所有納入研究的平均病程從而產生一個新變量×作為擬合模型的協變量。
從表 1的數據可以假設每種藥物對不同病程的類風濕關節炎患者干預效應不同，較長的病程可能獲得更高的治療應答。為檢驗病程對干預效應的影響，本文選擇兩個指標：一是經過協變量“校正”后的每種藥物的干預效應量，一是整個網狀中所有藥物療效“最佳”排序，前者可通過mvmeta命令實現，后者可由sucra命令實現 [9]。該研究發現，在控制了病程這一效應修飾因素后，藥物E的ACR50達標率降低，與參照藥物A（安慰劑）差異無統計學意義；而療效排序結果顯示，當模型中納入平均病程，藥物E療效“最佳”的概率由89.0%降為36.3%，提示其具有較高的報告應答率，其“更有效”的可能原因是納入的類風濕關節炎患者具有較長的病程，從而證實較長的病程可能獲得更高的治療應答這一假設。
網狀Meta回歸分析與簡單Meta回歸分析也會遇到相同的問題 [14]：以聚合數據為協變量時會產生生態偏倚（ecological bias）；納入研究少時統計效能較低；存在異質性時采用固定效應模型則假陽性率高等。
總之，如果NMA的同質性、相似性和一致性三個基本假設難以成立，或眾多的效應修飾因素（稱為協變量）可能會影響到干預效應量，則可以使用網狀Meta回歸解釋異質性與不一致性來源，探討效應修飾因素如何影響干預效應量等，但應注意其不足之處，以達到合理使用網狀Meta回歸分析的目的。
網狀Meta分析（network meta-analyses，NMA）可同時進行直接與間接比較，將一系列不同治療方法的隨機臨床試驗數據匯總，并就給定的治療終點進行點和可信區間估計，同時對不相關性進行評估 [1-3]，但治療終點結果可能會隨不同的人群或干預特征如劑量或療程而發生變化，要探討它們之間的關系，可將Meta回歸的思想擴展到NMA中，稱為網狀Meta回歸（Network meta-regression），它可以在頻率學或貝葉斯框架下實現。本文以實例介紹如何使用Stata軟件的mvmeta命令在頻率學框架下實現網狀Meta回歸，并簡單討論網狀Meta回歸的應用。
1 資料與方法
1.1 數據來源
本文所用數據來源于治療類風濕關節炎生物制劑的Meta分析 [4]，包含Placebo、Abatacept、Adalimumab、Anakinra、Etanercept、Golimumab、Infliximab、Rituximab和Tocilizumab共9種藥物，其納入患者的指標為美國風濕病學會ACR50標準（患者關節腫脹及觸痛的個數有50%改善），原文共納入31個研究。本研究針對報告了平均病程的28個研究（含1個三臂研究，27個兩臂研究）進行分析。數據按表 1所示格式輸入Stata數據管理器中，命名為acr50.dta存儲于C盤根目錄下，其中“studyid”表示各研究，“duration”表示每個研究的平均病程，A、B、C、D、E、F、G、H、I分別表示上述9種藥物，以n*和r*分別表示每個研究中每個臂總人數和達到ACR50的人數，缺失數據（空格）以圓點表示。
表1
治療類風濕關節炎生物制劑Meta分析的具體數據
1.2 網狀Meta回歸模型
White等 [5, 6]建立了多元隨機效應Meta回歸模型（multivariate random effects meta-regression model）可以用于網狀Meta回歸分析。假設第i個研究中有相對于某一共同參照有p個比較，相應p個效應估計值yi為（1×p）向量，其相應方差-協方差矩陣為Si，則模型為yi～N(μi,Si)，μi～N(βXi,∑)。式中，μi表示研究平均效應量，Xi表示研究水平協變量，∑為研究間方差-協方差矩陣。
1.3 模型擬合
該模型擬合可由Stata軟件mvmeta命令實現 [5, 7]，該命令下載與安裝可以按參考文獻 [3, 8]的方法進行，其命令行基本操作格式為 [5]：mvmeta b V xvars [if] [in] [,options]，具體解釋和使用方法可以參閱相關文獻 [3, 5]，本文不再贅述。
為方便計算，筆者編寫了.do文件，在擬合網狀Meta回歸模型的同時，根據累積排序曲線下面積（surface under the cumulative ranking curve，SUCRA）值對各種干預藥物的療效進行排序 [9]，其具體操作步驟和主要解釋見框1。有興趣的讀者，可以在Stata軟件的菜單欄中按“Windows → Do-file Editor → New Do-file Editor”步驟操作，打開“Do-file”編輯器，將框1中代碼復制到已打開的“Do-file”編輯器中，并在編輯器的“Tools”下拉菜單中選擇“Excute(do)”執行即可。
框 1 進行網狀Meta回歸分析、干預措施療效排序的Stata代碼
**第一步：讀入數據 use C:\acr50.dta,clear **第二步：數據預處理 *建立效應量、方差－協方差矩陣 foreach trt in A B C D E F G H I{ if “`trt’”==”A” continue gen y`trt’ = log(r`trt’/(n`trt’-r`trt’)) - log(rA/(nA-rA)) gen S`trt’`trt’ = 1/r`trt’ + 1/(n`trt’-r`trt’) + 1/rA + 1/(nA-rA) foreach trt2 in A B C D E F G H I{ if “`trt2’”==”A” continue if “`trt2’”＞”`trt’” gen S`trt’`trt2’ = 1/rA + 1/(nA-rA) if !mi(r`trt’) & !mi(r`trt2’) } } format y* S* %6.2g mat P8=0.5*(I(8) +J(8,8,1) ) **第三步：擬合Meta回歸模型，對干預藥物療效進行排序 *擬合Meta回歸模型 sum duration gen x=duration-r(mean) mvmeta y S x,bscov(prop P8) pbest(max,all zero gen(probm) reps(10000) ) eform *療效排序 sucra probm*,mvmeta lab(A B C D E F G H I) noplot
為考察平均病程對干預效應的影響，本次研究采用經典的NMA（擬合一致性模型），并進行療效排序。只需要將框1中第三步中的代碼做如下修改即可。
mvmeta y S,bscov(prop P8) pbest(max,all zero gen(prob) reps(10000) ) eform
sucra prob*,mvmeta lab(A B C D E F G H I) noplot
請注意，擬合一致性模型和網狀Meta回歸模型，均基于跨研究間干預措施配對比較異質性相同的假設。
2 結果
2.1 模型擬合結果
模型擬合的迭代過程從略，只報告主要結果。擬合一致性模型、網狀Meta回歸模型結果見表 2。由表 2可見，相對于藥物A（安慰劑），其他治療類風濕關節炎的生物制劑可以提高ACR50達標率，其差異均有統計學意義（P≤ 0.001），但經過平均病程協變量（即具有8.52年平均病程時的效應量）校正后，藥物E（Etanercept）與藥物A相比，ACR50達標率差異未再出現統計學意義（P=0.806），提示病程可能與藥物E存在交互作用。
表2
網狀Meta分析及網狀Meta回歸分析結果
2.2 療效排序結果
干預藥物療效排序結果見表 3，其中SUCRA是累積排序曲線下面積的數字化表示，“PrBest”表示每個干預藥物療效最佳的概率。從表 3可知，在經典NMA結果中，藥物E的SUCRA值最大，療效最佳的可能性為89.0%，但經過平均病程協變量校正后，其療效最佳的可能性降為36.3%，也提示平均病程對藥物作用存在影響。
表3
干預藥物療效排序
3 討論
醫學研究中，干預措施呈現多樣化發展，面對多種干預措施，兩兩比較的經典Meta分析已不能為選擇最佳干預措施提供方法學支持，但NMA使得多種干預措施之間的比較成為可能 [3, 10]。NMA涉及同質性、相似性和一致性三個基本假設 [3, 11-13]，雖然這些假設各有自己的不同術語和定義 [14]，但實質上均是經典Meta分析同質性思想在NMA的擴展應用，因此本文所指的“異質性”可以解釋為異質性和/或不一致性 [15]。
NMA包括多個不同干預比較的不同試驗，因此較之于經典的Meta分析在臨床、方法學等方面更可能存在多樣性：如患者、干預措施、研究地點、隨訪時間等臨床特征，以及隨機、盲法、分配隱藏等研究設計方法等各種效應修飾因素，不僅存在于某對特定干預比較之中，也存在不同干預比較之間 [3, 13]，必然會對NMA結果的真實性產生影響。
Meta回歸常被用于探討Meta分析中潛在的效應修飾因素（稱為協變量）對干預效應量的影響，有眾多的Meta回歸模型可用于NMA中 [15]，實現方法可以基于頻率學或貝葉斯框架。貝葉斯框架，如分層NMA模型納入協變量可以由WinBUGS軟件擬合，但可能因選擇不同先驗需要進行敏感性分析，或易出現蒙特卡洛誤差；而頻率學方法，如White提出的多元Meta回歸方法則可避免上述不足，且具有計算速度快的優點 [16]，可使用Stata軟件的mvmeta命令輕松實現，本文即采用此方法。
協變量可以根據試驗的臨床特征，或根據臨床設計的偏倚風險等來選擇。協變量數據類型可分二分類或連續型，但連續型協變量需要采用中心化協變量值 [15]，也就是對其進行測量中心化（measurement centering），測量中心化方法分為總均數中心化和組均數中心化兩種。針對所引用的數據實例，選擇納入每個研究的平均病程這一試驗臨床特征為協變量，采用總均數中心化方法，即每個研究的平均病程減去所有納入研究的平均病程從而產生一個新變量×作為擬合模型的協變量。
從表 1的數據可以假設每種藥物對不同病程的類風濕關節炎患者干預效應不同，較長的病程可能獲得更高的治療應答。為檢驗病程對干預效應的影響，本文選擇兩個指標：一是經過協變量“校正”后的每種藥物的干預效應量，一是整個網狀中所有藥物療效“最佳”排序，前者可通過mvmeta命令實現，后者可由sucra命令實現 [9]。該研究發現，在控制了病程這一效應修飾因素后，藥物E的ACR50達標率降低，與參照藥物A（安慰劑）差異無統計學意義；而療效排序結果顯示，當模型中納入平均病程，藥物E療效“最佳”的概率由89.0%降為36.3%，提示其具有較高的報告應答率，其“更有效”的可能原因是納入的類風濕關節炎患者具有較長的病程，從而證實較長的病程可能獲得更高的治療應答這一假設。
網狀Meta回歸分析與簡單Meta回歸分析也會遇到相同的問題 [14]：以聚合數據為協變量時會產生生態偏倚（ecological bias）；納入研究少時統計效能較低；存在異質性時采用固定效應模型則假陽性率高等。
總之，如果NMA的同質性、相似性和一致性三個基本假設難以成立，或眾多的效應修飾因素（稱為協變量）可能會影響到干預效應量，則可以使用網狀Meta回歸解釋異質性與不一致性來源，探討效應修飾因素如何影響干預效應量等，但應注意其不足之處，以達到合理使用網狀Meta回歸分析的目的。
