當前已發表的介紹網狀Meta分析的方法學文章主要集中于二分類變量，尚無針對連續性變量的相關報道。本文介紹如何使用R、GeMTC和Stata軟件聯合實現連續性變量網狀Meta分析的計算及繪制收斂診斷圖、軌道密度圖、森林圖、等級概率圖、內部關系總結圖、網狀關系圖、貢獻圖和發表偏倚檢測圖。
引用本文： 趙坤, 陳凌霄, 田金徽, 曾憲濤, 張超, 拜爭剛. 使用R、GeMTC和STATA軟件實現連續變量的網狀Meta分析. 中國循證醫學雜志, 2015, 15(7): 861-868. doi: 10.7507/1672-2531.20150145 復制
網狀Meta分析（network meta-analysis，NMA）因可將間接比較（indirect comparison）和直接比較（head to head comparison）的證據進行合并，更好地滿足臨床實踐的需求，而日益受到青睞[1]。當前能夠實現NMA軟件有很多，但各有利弊[2]：ADDIS軟件操作簡單、數據處理能力強、繪圖較為全面，但數據輸入十分繁瑣[3]；Microsoft Excel軟件的NetMetaXL宏命令將數據提取與輸入、結果輸出和圖形繪制集于一體，但僅適用于二分類變量[4]；BUGS軟件處理能力強，但要求使用者有較好的編程功底且采用其他軟件調用時也容易出錯，初學者難以掌握[5-9]。
R軟件是一套完整的數據處理、計算和繪圖軟件系統，是屬于GNU系統的自由、免費、源代碼開放的軟件，可通過自身程序包或調用其他軟件完成NMA[10]。GeMTC（Generate Mixed Treatment Comparisons）軟件是由van Valkenhoef等研發的MTC軟件改編而來的專用于制作NMA的非編程軟件，但其已停止更新，且不能制作網狀關系圖，數據僅顯示小數點后兩位，有一定的局限性[11]。目前國內已發表的方法學文章主要集中于介紹二分類變量的NMA[3-11]，尚無介紹連續性變量的NMA的論文。本文將介紹如何實現連續型變量的NMA，所選擇的軟件為GeMTC、R和Stata。
1 聯合GeMTC和R軟件實現
1.1 軟件安裝及程序包加載
R軟件當前最新版本為R-3.1.2，獲取網址為http://www.r.project.org；GeMTC軟件最新版本為0.14.3，獲取網址為http://drugis.org/software/addis1/gemtc。GeMTC軟件在運行時需要Java環境，Java可從http://www.Java.com/zh_CN/下載。本文采用OpenBUGS軟件執行數據分析，因此還需從http://www.openbugs.net/w/Downloads下載該軟件。所有軟件下載完成后，雙擊安裝包進行安裝。
軟件安裝完畢后，雙擊桌面的R軟件圖標啟動R的交互式窗口（R—GUI）；接著在菜單欄處選擇程序包下的安裝程序包→選擇CRAN鏡像→gemtc，安裝完成后可由library(gemtc)命令完成gemtc程序包加載。此外還需加載coda程序包、lattice程序包，具體方法如上。
1.2 數據處理及錄入
本文以Multiple-treatments Meta-analysis of a network of interventions[12]的數據為例進行演示，處理后的數據如表 1所示。數據處理完成后，采用GeMTC軟件錄入（操作簡便且不易出錯）。
表1
使用數據
1.2.1 錄入干預措施組
選擇Data下的Treatment進行干預措施的數據輸入，可點擊紅色框區的Add treatment進行添加、編輯與刪除，ID編號與表格一致，輸入1，具體描述為文法拉辛（圖 1）。西酞普蘭與帕羅西汀的添加方法與此相同。
圖1
干預措施的錄入
1.2.2 錄入研究組
選擇Data下的Study進行研究數據輸入，可點擊Study進行添加，ID編號與表格一致，可知1號研究對應的干預措施為1和2，因此ID為1，選擇干預措施1、2并打勾。在紅色框區選擇數據類型并對該數據集命名，本例選擇CONTINOUS（連續性變量），如圖 2所示。其余研究同理進行添加。
圖2
研究組數據錄入
1.2.3 錄入結局指標
連續型變量錄入的是mean、sd和n。圖 3為數據錄入完成后的界面，需注意的是ID為自動排序，因此編號為10的排列第二，且數據列表中干預措施的排列順序與在研究組錄入時的點擊順序有關。錄入數據完成后，點擊Save進行保存及文件命名，本例命名為depression。
圖3
數據完成界面
1.3 數據分析
采用R軟件gemtc程序包調用OpenBUGS軟件進行NMA以及相關圖形的繪制。
1.3.1 讀取錄入數據
讀取數據的命令為：
network＜-read.mtc.network("C:/Users/Administrator/Desktop/depression.gemtc")
注意括號內容應與1.2.3的文件命名及保存位置一致，命令前的network＜-相當于賦值，不能省略，命令符號應用英文輸入狀態且不應有空格，以下同理。如需查看數據可輸入“print(network)”命令。
1.3.2 設置model
設置model的命令為：model＜- mtc.model(network,type = "consistency",n.chain = 3) 。其中network為network數據，type為是否選取一致性模型（注：在模型的選擇上，通常認為由一致性模型得出的結果更加接近真實情況），n.chain為迭代運算中鏈的條數（一般來說，保證n≥3的鏈條數目使結果更加可信）。
1.3.3 分析數據
在加載程序包之前，應確保openbugs安裝在C盤，否則不能調用。執行數據分析的命令為：
results＜-mtc.run(model,sampler ="BRugs",n.adapt = 5000,n.iter = 20000,thin = 1)
命令中，“ n.adapt ”為預迭代次數，“ n.iter ”為迭代運算次數，“ thin ”為步長。
1.4 匯總結果
輸入命令summary(results)，查看均數差與百分數，結果如圖 4所示。
圖4
Meta分析圖
1.5 圖形繪制
1.5.1 收斂診斷圖
收斂診斷圖（convergence diagnostics）的繪制方法為通過輸入命令“gelman.plot(results)”實現，圖形如圖 5所示。保存格式及位置可以自行選擇。
圖5
收斂診斷圖
1.5.2 森林圖
森林圖（forest）的繪制命令為forest(results)，命令執行后結果見圖 6。
圖6
森林圖
1.5.3 軌跡密度圖
軌跡密度圖（trace and density）的實現命令為“plot(results)”，結果如圖 7。
圖7
軌道密度圖
1.5.4 等級概率圖
等級概率圖（rankk probability）分別有以下三種表示方式（圖 8、圖 9），分別對應的命令為1.ranks＜- rank.probability(results)，輸入命令print(ranks)查看。
2.plot(ranks)
3.plot(ranks,beside=TRUE)
圖8
風險概率圖
圖9
軌道密度圖
1.5.5 內部關系總結圖
輸入命令summary (network)，可查看每種干預措施的研究數量以及試驗臂數對應的研究數量。如圖 10所示，本研究中研究文法拉辛的研究有8篇，2臂試驗9篇，3臂試驗1篇。
圖10
內部關系總結圖
1.5.6 網狀關系圖
因為利用plot(network)制作出的網狀關系圖的命令僅反映該研究之間的直接比較關系，無法通過頂點反映研究數量的變化，因此可選擇R軟件或Stata繪制網狀關系圖（network plot）[9]，因Stata還可制作貢獻圖（contribution plot）及發表偏倚檢測圖（publication bias test），因此本文使用Stata進行連續性變量的NMA的詳細介紹。
2 應用Stata軟件實現
2.1 Stata下載與操作
用戶可從http://www.stata.com/stata13/購買最新版Stata 13.0，本文使用的為該版本[
2.2 模塊加載
加載metan模塊的命令為：ssc install metan。加載圖形分析network graphs模塊的命令為：net from http://www.mtm.uoi.gr。點擊network_graphs后按照提示進行安裝。
2.3 數據錄入
首先需要將表 1中的數據整理成2臂試驗的格式[7]，整理后的數據如表 2所示。整理完成后，輸入到Stata中。
表2
Stata應用數據格式
2.4 圖形繪制
2.4.1 網狀關系圖
在命令輸入框輸入networkplot t1 t2[13]，如圖 11，頂點與線條的粗細均表示研究數量的多少，與樣本量無關，這應與R軟件相區別[10]。用戶可以使用Graph Editor進行網狀控制圖的編輯，調整字體的大小、線條的顏色等。
圖11
網狀關系圖
2.4.2 貢獻圖
輸入命令metan n1 x1 s1 n2 x2 s2，可得出各個研究的SMD以及seSMD，在Stata中表示為_ES,_seES。輸入命令netweight _ES _seES t1 t2[13]，結果如圖 12。
圖12
貢獻圖
2.4.3 發表偏倚檢測圖
輸入命令netfunnel _ES _seES t1 t2,bycomparison[13]，結果如圖 13。
圖13
發表偏倚檢測圖
3 結語
盡管目前關于NMA的方法學文獻不勝枚舉，但僅限于以二分類變量作為示例進行講解，對于初學者而言，連續性變量的NMA方法無疑成為一大難題。在市面上主流的軟件中，WinBUGS/OpenBUGS盡管功能強大且變化靈活，可惜其用戶界面不夠友好，錄入數據繁瑣易錯，對初學者不易上手，且作圖功能薄弱；Stata可進行頻率學派的運算，不失為一種較好的方法，但非統計專業出身的醫學生對于協方差方差矩陣的構建理解起來具有一定困難且其運算結果的穩健程度尚未經過大規模驗證[14-16]，不建議初學者先行學習，但是Sanlanti團隊開發的畫圖命令堪稱現階段最易上手的方法[13]；R軟件免費開源，不涉及版權問題，統計作圖功能全面強大，通過GeMTC包調用相關的貝葉斯軟件進行運算，上手雖有一定難度，但比之原軟件要容易不少。本文將R軟件與Stata軟件聯合使用，首先利用GeMTC輸入連續性變量數據并將其錄入R軟件中，然后通過輸入命令制作Meta分析圖、收斂診斷圖、軌道密度圖、森林圖等，由于R軟件的局限性，本文又使用Stata軟件制作網狀控制圖、貢獻圖以及發表偏倚檢測圖，這種方法不僅節省了時間，而且使NMA的圖形更加全面，以供研究者選擇合適的圖形用于發表論文使用。但這種方法的局限性在于GeMTC軟件錄入數據僅適用于小樣本的數據，對于大樣本的數據還需使用R軟件[11]。
本文結合GeMTC簡單的錄入數據操作與R軟件強勁的統計分析與繪圖能力，同時利用Stata詳細介紹了連續性變量的NMA的圖形繪制，希望能夠對完善NMA的方法學有所幫助。
網狀Meta分析（network meta-analysis，NMA）因可將間接比較（indirect comparison）和直接比較（head to head comparison）的證據進行合并，更好地滿足臨床實踐的需求，而日益受到青睞[1]。當前能夠實現NMA軟件有很多，但各有利弊[2]：ADDIS軟件操作簡單、數據處理能力強、繪圖較為全面，但數據輸入十分繁瑣[3]；Microsoft Excel軟件的NetMetaXL宏命令將數據提取與輸入、結果輸出和圖形繪制集于一體，但僅適用于二分類變量[4]；BUGS軟件處理能力強，但要求使用者有較好的編程功底且采用其他軟件調用時也容易出錯，初學者難以掌握[5-9]。
R軟件是一套完整的數據處理、計算和繪圖軟件系統，是屬于GNU系統的自由、免費、源代碼開放的軟件，可通過自身程序包或調用其他軟件完成NMA[10]。GeMTC（Generate Mixed Treatment Comparisons）軟件是由van Valkenhoef等研發的MTC軟件改編而來的專用于制作NMA的非編程軟件，但其已停止更新，且不能制作網狀關系圖，數據僅顯示小數點后兩位，有一定的局限性[11]。目前國內已發表的方法學文章主要集中于介紹二分類變量的NMA[3-11]，尚無介紹連續性變量的NMA的論文。本文將介紹如何實現連續型變量的NMA，所選擇的軟件為GeMTC、R和Stata。
1 聯合GeMTC和R軟件實現
1.1 軟件安裝及程序包加載
R軟件當前最新版本為R-3.1.2，獲取網址為http://www.r.project.org；GeMTC軟件最新版本為0.14.3，獲取網址為http://drugis.org/software/addis1/gemtc。GeMTC軟件在運行時需要Java環境，Java可從http://www.Java.com/zh_CN/下載。本文采用OpenBUGS軟件執行數據分析，因此還需從http://www.openbugs.net/w/Downloads下載該軟件。所有軟件下載完成后，雙擊安裝包進行安裝。
軟件安裝完畢后，雙擊桌面的R軟件圖標啟動R的交互式窗口（R—GUI）；接著在菜單欄處選擇程序包下的安裝程序包→選擇CRAN鏡像→gemtc，安裝完成后可由library(gemtc)命令完成gemtc程序包加載。此外還需加載coda程序包、lattice程序包，具體方法如上。
1.2 數據處理及錄入
本文以Multiple-treatments Meta-analysis of a network of interventions[12]的數據為例進行演示，處理后的數據如表 1所示。數據處理完成后，采用GeMTC軟件錄入（操作簡便且不易出錯）。
表1
使用數據
1.2.1 錄入干預措施組
選擇Data下的Treatment進行干預措施的數據輸入，可點擊紅色框區的Add treatment進行添加、編輯與刪除，ID編號與表格一致，輸入1，具體描述為文法拉辛（圖 1）。西酞普蘭與帕羅西汀的添加方法與此相同。
圖1
干預措施的錄入
1.2.2 錄入研究組
選擇Data下的Study進行研究數據輸入，可點擊Study進行添加，ID編號與表格一致，可知1號研究對應的干預措施為1和2，因此ID為1，選擇干預措施1、2并打勾。在紅色框區選擇數據類型并對該數據集命名，本例選擇CONTINOUS（連續性變量），如圖 2所示。其余研究同理進行添加。
圖2
研究組數據錄入
1.2.3 錄入結局指標
連續型變量錄入的是mean、sd和n。圖 3為數據錄入完成后的界面，需注意的是ID為自動排序，因此編號為10的排列第二，且數據列表中干預措施的排列順序與在研究組錄入時的點擊順序有關。錄入數據完成后，點擊Save進行保存及文件命名，本例命名為depression。
圖3
數據完成界面
1.3 數據分析
采用R軟件gemtc程序包調用OpenBUGS軟件進行NMA以及相關圖形的繪制。
1.3.1 讀取錄入數據
讀取數據的命令為：
network＜-read.mtc.network("C:/Users/Administrator/Desktop/depression.gemtc")
注意括號內容應與1.2.3的文件命名及保存位置一致，命令前的network＜-相當于賦值，不能省略，命令符號應用英文輸入狀態且不應有空格，以下同理。如需查看數據可輸入“print(network)”命令。
1.3.2 設置model
設置model的命令為：model＜- mtc.model(network,type = "consistency",n.chain = 3) 。其中network為network數據，type為是否選取一致性模型（注：在模型的選擇上，通常認為由一致性模型得出的結果更加接近真實情況），n.chain為迭代運算中鏈的條數（一般來說，保證n≥3的鏈條數目使結果更加可信）。
1.3.3 分析數據
在加載程序包之前，應確保openbugs安裝在C盤，否則不能調用。執行數據分析的命令為：
results＜-mtc.run(model,sampler ="BRugs",n.adapt = 5000,n.iter = 20000,thin = 1)
命令中，“ n.adapt ”為預迭代次數，“ n.iter ”為迭代運算次數，“ thin ”為步長。
1.4 匯總結果
輸入命令summary(results)，查看均數差與百分數，結果如圖 4所示。
圖4
Meta分析圖
1.5 圖形繪制
1.5.1 收斂診斷圖
收斂診斷圖（convergence diagnostics）的繪制方法為通過輸入命令“gelman.plot(results)”實現，圖形如圖 5所示。保存格式及位置可以自行選擇。
圖5
收斂診斷圖
1.5.2 森林圖
森林圖（forest）的繪制命令為forest(results)，命令執行后結果見圖 6。
圖6
森林圖
1.5.3 軌跡密度圖
軌跡密度圖（trace and density）的實現命令為“plot(results)”，結果如圖 7。
圖7
軌道密度圖
1.5.4 等級概率圖
等級概率圖（rankk probability）分別有以下三種表示方式（圖 8、圖 9），分別對應的命令為1.ranks＜- rank.probability(results)，輸入命令print(ranks)查看。
2.plot(ranks)
3.plot(ranks,beside=TRUE)
圖8
風險概率圖
圖9
軌道密度圖
1.5.5 內部關系總結圖
輸入命令summary (network)，可查看每種干預措施的研究數量以及試驗臂數對應的研究數量。如圖 10所示，本研究中研究文法拉辛的研究有8篇，2臂試驗9篇，3臂試驗1篇。
圖10
內部關系總結圖
1.5.6 網狀關系圖
因為利用plot(network)制作出的網狀關系圖的命令僅反映該研究之間的直接比較關系，無法通過頂點反映研究數量的變化，因此可選擇R軟件或Stata繪制網狀關系圖（network plot）[9]，因Stata還可制作貢獻圖（contribution plot）及發表偏倚檢測圖（publication bias test），因此本文使用Stata進行連續性變量的NMA的詳細介紹。
2 應用Stata軟件實現
2.1 Stata下載與操作
用戶可從http://www.stata.com/stata13/購買最新版Stata 13.0，本文使用的為該版本[
2.2 模塊加載
加載metan模塊的命令為：ssc install metan。加載圖形分析network graphs模塊的命令為：net from http://www.mtm.uoi.gr。點擊network_graphs后按照提示進行安裝。
2.3 數據錄入
首先需要將表 1中的數據整理成2臂試驗的格式[7]，整理后的數據如表 2所示。整理完成后，輸入到Stata中。
表2
Stata應用數據格式
2.4 圖形繪制
2.4.1 網狀關系圖
在命令輸入框輸入networkplot t1 t2[13]，如圖 11，頂點與線條的粗細均表示研究數量的多少，與樣本量無關，這應與R軟件相區別[10]。用戶可以使用Graph Editor進行網狀控制圖的編輯，調整字體的大小、線條的顏色等。
圖11
網狀關系圖
2.4.2 貢獻圖
輸入命令metan n1 x1 s1 n2 x2 s2，可得出各個研究的SMD以及seSMD，在Stata中表示為_ES,_seES。輸入命令netweight _ES _seES t1 t2[13]，結果如圖 12。
圖12
貢獻圖
2.4.3 發表偏倚檢測圖
輸入命令netfunnel _ES _seES t1 t2,bycomparison[13]，結果如圖 13。
圖13
發表偏倚檢測圖
3 結語
盡管目前關于NMA的方法學文獻不勝枚舉，但僅限于以二分類變量作為示例進行講解，對于初學者而言，連續性變量的NMA方法無疑成為一大難題。在市面上主流的軟件中，WinBUGS/OpenBUGS盡管功能強大且變化靈活，可惜其用戶界面不夠友好，錄入數據繁瑣易錯，對初學者不易上手，且作圖功能薄弱；Stata可進行頻率學派的運算，不失為一種較好的方法，但非統計專業出身的醫學生對于協方差方差矩陣的構建理解起來具有一定困難且其運算結果的穩健程度尚未經過大規模驗證[14-16]，不建議初學者先行學習，但是Sanlanti團隊開發的畫圖命令堪稱現階段最易上手的方法[13]；R軟件免費開源，不涉及版權問題，統計作圖功能全面強大，通過GeMTC包調用相關的貝葉斯軟件進行運算，上手雖有一定難度，但比之原軟件要容易不少。本文將R軟件與Stata軟件聯合使用，首先利用GeMTC輸入連續性變量數據并將其錄入R軟件中，然后通過輸入命令制作Meta分析圖、收斂診斷圖、軌道密度圖、森林圖等，由于R軟件的局限性，本文又使用Stata軟件制作網狀控制圖、貢獻圖以及發表偏倚檢測圖，這種方法不僅節省了時間，而且使NMA的圖形更加全面，以供研究者選擇合適的圖形用于發表論文使用。但這種方法的局限性在于GeMTC軟件錄入數據僅適用于小樣本的數據，對于大樣本的數據還需使用R軟件[11]。
本文結合GeMTC簡單的錄入數據操作與R軟件強勁的統計分析與繪圖能力，同時利用Stata詳細介紹了連續性變量的NMA的圖形繪制，希望能夠對完善NMA的方法學有所幫助。
