BUGSnet程序包是一款功能強大、用于實現貝葉斯網狀Meta分析的R語言程序包。該程序包基于JAGS軟件并可根據公認的報告規范（PRISMA、ISPOR-AMPC-NCA和NICE-DSU）進行高質量的貝葉斯網絡Meta分析。本文借助類固醇佐劑治療天皰瘡的網狀Meta分析實例展示BUGSnet程序包處理連續性數據或二分類數據實現貝葉斯網狀Meta分析的操作流程。
引用本文： 向虹, 岳磊, 趙紅, 李華, 謝天宏, 楊婷, 王杰, 張宇昊, 謝忠平. 應用R語言BUGSnet程序包實現貝葉斯網狀Meta分析. 中國循證醫學雜志, 2022, 22(5): 600-608. doi: 10.7507/1672-2531.202111140 復制
近年來，網狀Meta分析（network meta-analysis，NMA）發展迅速，研究者針對不同類型的研究數據開發出對應的解決方案，NMA分析變得越來越簡單，相關發表文獻也呈逐年遞增趨勢[1-2]。Tian等[3]統計了Web of Science數據庫中納入研究NMA方法學文章的引用次數，發現引用次數最高的是Hutton等[4]發表的PRISMA-NMA報告指南，該指南可幫助Meta作者規范撰寫報告并提高報告的完整性。引用次數第2高的是Lu等[5]使用馬爾可夫鏈蒙特卡羅（Markov chain Monte Carlo，MCMC）[6]方法提出的一種分層貝葉斯模型，該模型解決了間接證據之間的比較問題，為NMA貝葉斯分析奠定了基礎。由此可見，NMA作者關心的是NMA分析方法能否按照指南標準對納入研究的數據進行規范和完整性的研究，進而基于證據提供可靠的決策。盡管已經開發出一些NMA分析軟件，但由于存在收費、使用困難、功能不完全、無具體指南標準和可視化效果不良等問題而可能受到使用限制[7-9]。為更好地實施NMA分析報告，Béliveau等[10]創建了一個R程序包BUGSnet（Bayesian inference using gibbs sampling to conduct a network meta-analysis）。該程序包是一款基于R軟件平臺的開源、免費、功能齊備的軟件包，能夠實現貝葉斯NMA分析，并滿足當前Meta分析報告規范的要求。BUGSnet程序包能夠處理不同類型的資料包括計數數據、連續性數據或二分類數據，可用于描述證據網狀、估計模型擬合及評估異質性和不一致性等。此外，該程序包定期更新，并在GitHub中提供使用幫助。本文以實例展示該程序包的具體操作。
1 資料與方法
1.1 數據來源
數據來源于Lee等[11]關于類固醇佐劑治療天皰瘡的NMA研究。該研究納入10個試驗，包含592例患者，評估了7種保留類固醇佐劑的方法治療天皰瘡的效果（表1）。結局指標包含二分類數據（疾病緩解病例數、疾病復發病例數和不良事件相關退出病例數）和連續性數據（類固醇的累積使用劑量）。
表1
納入研究的基本特征[11]
1.2 軟件及程序包的安裝
本文使用的R軟件版本為R x64 4.1.2（下載地址：https://www.r-project.org）。使用R軟件控制臺R Studio操作（下載地址：https://www.rstudio.com/products/rstudio/）。BUGSnet程序包依賴JAGS（Just Another Gibbs Sampler）的廣義線性模型實現貝葉斯NMA，因此需調用JAGS-4.3.0軟件（下載地址：https://sourceforge.net/projects/mcmc-jags/）。上述軟件下載完成后，雙擊軟件安裝包可完成軟件安裝。
完成上述軟件安裝后，打開R Studio，在腳本區輸入以下代碼安裝Rtool。
install.packages("pkgbuild")
pkgbuild::has_build_tools()
接著，在R Studio腳本區輸入以下代碼完成BUGSnet程序包的安裝。
install.packages(c("remotes", "knitr"))
remotes::install_github("audrey-b/BUGSnet@v1.1.0", upgrade = TRUE, build_vignettes = TRUE, dependencies = TRUE)
安裝“xlsx”、“dplyr”和“tidyr”程序包以便于數據的錄入、整理和變換。
install.packages(c("xlsx", "dplyr", "tidyr"))
1.3 數據的錄入和準備
數據按表1所示，將中文漢字改為對應的英文后整理為Excel表格格式，命名為“data_total.xlsx”存儲于演示目錄下。
加載以下程序包準備數據的變換。
library(xlsx)
library(dplyr)
library(tidyr)
輸入以下代碼完成數據錄入。
data_total<-read.xlsx("data_total.xlsx", 1)
輸入以下代碼完成不同結局指標的數據提取。
data_remission<- data_total[, c(1:4)]
data_relapse<-data_total[, c(1:3, 5)]
data_withdrawal<- data_total[, c(1:3, 6)]
data_dose<-data_total[, c(1:3, 7:8)]
輸入以下代碼規范數據錄入格式，即去除缺失值、并將數據存儲到對應的結局指標數據集。
data_remission<-data_remission %>% rename(sampleSize=Number.of.patients, events= Number.of.remissions)
data_relapse<-data_relapse %>% drop_na() %>% rename(sampleSize=Number.of.patients, events=Number.of.relapses)
data_withdrawal<-data_withdrawal %>% drop_na() %>% rename(sampleSize= Number.of.patients, events=Number.of.adverse.event.related.withdrawals)
data_dose<-data_dose %>% drop_na() %>% rename(sampleSize=Number.of.patients, mean=Cumulative.dose.of.steroid_Mean, SD=Cumulative.dose.of.steroid_SD)
在上述代碼中，“data_remission”代表緩解病例數據，為主要療效指標；“data_relapse”代表復發病例數據，為次要療效指標；“data_withdrawal”代表不良事件相關退出病例數據，為次要安全性指標；“data_dose”代表固醇的累積使用劑量，為主要安全性指標，與前三個指標的二分類數據類型不同的是，該指標類型為連續數據。
1.4 BUGSnet程序包數據分析
在完成數據錄入后，正式進入數據分析。如表2所示，BUGSnet程序包的函數可執行包括數據預處理、各種數據分析及圖形或表格結果的輸出等命令。
表2
BUGSnet程序包可實現的功能
在完成數據載入后，可通過BUGSnet程序包將整理好的數據進行NMA分析。在RStudio腳本區輸入以下代碼完成BUGSnet程序包的加載和變量名稱的賦值，以便于進一步分析處理。由于“data_remission”、“data_relapse”和“data_withdrawal”均為二分類數據類型，數據分析方法相似，為避免重復，在后續代碼演示時以“data_remission”為例。
library(BUGSnet)
data_remission <- data.prep(arm.data = data_remission, varname.t = "Treatment", varname.s = "Study")
data_dose <- data.prep(arm.data = data_dose, varname.t = "Treatment", varname.s = "Study")
使用net.tab()函數或net.plot()函數獲得證據網狀的描述性統計表格或圖形。“type.outcome”參數指示數據的類型，可根據結局指標的數據類型進行設置：二分類數據（binomial），連續性數據（continuous）和計數數據（rate/rate2）。
輸入以下代碼獲得證據網狀關系圖（圖1），參數“node.colour”可指定節點的顏色。
圖1
net.plot()函數生成網絡關系圖
a：疾病緩解，b：累積類固醇劑量。
par(mfrow = c(1,2))
net.plot(data_remission, node.colour = "dark grey")
net.plot(data_dose, node.colour = "dark grey")
nma.model()函數用于設置擬合模型，用戶可選擇使用固定效應或隨機效應模型。固定效應模型代碼如下：
fixed_effects_model_remission <- nma.model(data_remission, outcome="events", N="sampleSize", reference="Steroid alone", family="binomial",link="logit", effects="fixed")
fixed_effects_model_dose <- nma.model(data=data_dose, outcome="mean", N="sampleSize", sd="SD", reference="Steroid alone", family="normal", link="identity", effects="fixed")
隨機效應模型代碼如下：
random_effects_model_remission <- nma.model(data_remission, outcome="events", N="sampleSize", reference="Steroid alone", family="binomial", link="logit", effects="random")
random_effects_model_dose <- nma.model(data=data_dose, outcome="mean", N="sampleSize", sd="SD", reference="Steroid alone", family="normal", link="identity", effects="random")
上述命令中，參數“outcome”用于指定結局變量；參數“reference”用于設置對照組，通常是某種安慰劑或對照藥物等；參數“family”用于指定結果的分布類型，通常有“binomial”（二分類數據），“normal”（連續性數據），“poisson”（計數數據）可供選擇；參數“link”用于指定NMA模型使用的函數，其中“logit”用于二分類數據的比值比（Odds Ratio），“log”用于二分類數據的風險比（Risk Ratio）或計數數據的比率（Rate Ratio）；“cloglog”用于二分類數據的危險比（Hazard Ratio）；“identity”用于連續性數據；“effects”用于設置效應模型的類型。
生成NMA模型后，運行nma.run()函數命令執行貝葉斯NMA，該函數運行需要用戶指定MCMC算法的老化（n.burnin）、迭代（n.iter）和適應次數（n.adapt）及監控的變量。
對于固定效應模型的運行代碼如下。
set.seed(20222022)
fixed_effects_results_remission <- nma.run(fixed_effects_model_remission, n.adapt=1000, n.burnin=1000, n.iter=10000)
fixed_effects_results_dose <- nma.run(fixed_effects_model_dose, n.adapt=1000, n.burnin=1000, n.iter=10000)
對于隨機效應模型的運行代碼如下。
random_effects_results_remission <- nma.run(random_effects_model_remission, n.adapt=1000, n.burnin=1000, n.iter=10000)
random_effects_results_dose <- nma.run(random_effects_model_dose, n.adapt=1000, n.burnin=1000, n.iter=10000)
以下代碼可進行NMA分析的效應模型評估，即評估模型的擬合和識別潛在異常值，并輸出杠桿圖（圖2）。
圖2
nma.fit()函數生成的圖和擬合統計數據
leverage
par(mfrow = c(2,2))
對于固定效應模型的擬合代碼如下。
nma.fit(fixed_effects_results_remission)
nma.fit(fixed_effects_results_dose)
對于隨機效應模型的擬合代碼如下。
nma.fit(random_effects_results_remission)
nma.fit(random_effects_results_dose)
1.5 圖形繪制
1.5.1 網狀關系圖
網狀關系圖的繪制如前述，結果見圖1。
1.5.2 治療排名結果圖
使用nma.rank()函數可獲得治療排名的表格和圖形結果。將擬合最優模型結果輸入到該函數中即可獲得對應的排名結果，可參考以下代碼輸出結果，參數“largerbetter”若為“TRUE”則代表更大的結果與更好的治療效果相關，反之亦然。
sucra_out_remission<- nma.rank(random_effects_results_remission, largerbetter=TRUE)
sucra_out_dose<- nma.rank(fixed_effects_results_dose, largerbetter=FALSE)
nma.rank()函數自帶變量顏色，若各個組結局指標的變量名及其數量相同，則可直接在nma.rank()函數中將參數“sucra.palette”設定為所需的顏色集。如果“sucra_out_remission”和“sucra_out_dose”的“Treatment”變量總數及變量名均不完全相同，為統一各個圖片的圖注顏色，可使用以下代碼對“Treatment”中的不同變量設置特定的顏色。指定顏色需用到ggplot2程序包的“scale_colour_manual()”函數，如未安裝ggplot2程序包可使用“install.packages(“ggplot2”)”代碼進行安裝。
ggplot2程序包的加載代碼如下：
library(ggplot2)
各個變量顏色的賦值代碼如下：
z<-scale_colour_manual(values=c(AZA="#FF6600", MMF="#FF0088", "Steroid alone"="#00FF00", CP_P="#007700", Cyclosporine="#6633CC", "DCP_C(6M)"="#00CCFF", "DCP_C(12M)"="#FFFF00", Rituximab="#AA0000"))
繪制累積排序概率圖下面積（surface under the cumulative ranking，SUCRA）的代碼如下：
s1<- sucra_out_remission$sucraplot+z
s4<- sucra_out_dose$sucraplot+z
繪制排序概率圖（rankogram）的代碼如下：
r1<- sucra_out_remission$rankogram+theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
r4<- sucra_out_dose$rankogram+theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
輸入以下代碼，完成上述圖片的輸出（圖3）。
圖3
nma.rank()函數生成的SUCRA圖和rankogram圖
library(gridExtra)
grid.arrange(s1, s4, r1, r4, ncol = 2)
上述治療排名結果圖也可表格形式輸出，例如對于疾病緩解數據“sucra_out_remission”，不僅可使用“sucra_out_remission$sucratable”輸出SUCRA圖對應的表格，也可使用“sucra_out_remission$ranktable”輸出排序概率圖對應的表格。
1.5.3 排名表熱圖
nma.league()函數可獲得排名表熱圖，可比較任何一對干預之間的相對效果。將擬合最優模型結果輸入到nma.league()函中，如下面代碼所示：參數“central.tdcy”是用來衡量相對有效性的統計數據，可選擇“mean”或“median”；參數“order”表示圖片顯示的治療順序，本例以1.5.2治療排名結果圖中的輸出結果來排序；參數“log.scale”代表數據呈現方式，如果為“TRUE”則表示以對數刻度報告結果；參數“low.colour”、“mid.colour”及“high.colour”可設置熱圖的低、中、高值的顏色；參數“digits”設定小數點后的顯示位數，相關代碼設置如下。
league.out_remission<- nma.league(random_effects_results_remission, central.tdcy="median", order = sucra_out_remission$order, log.scale =TRUE, low.colour = "springgreen4", mid.colour = "white", high.colour = "red", digits = 2)
league.out_dose<- nma.league(fixed_effects_results_dose, central.tdcy="median", order = sucra_out_dose$order, log.scale =TRUE, low.colour = "springgreen4", mid.colour = "white", high.colour = "red", digits = 2)
使用以下代碼將排名表熱圖輸出結果放置在同一頁，輸出結果見圖4。
圖4
nma.league()函數生成排名表熱圖
l1<-league.out_remission$heatplot
l4<-league.out_dose$heatplot
grid.arrange(l1, l4, ncol = 1)
排名表熱圖也可表格的形式輸出，以疾病緩解數據為例，使用以下代碼可輸出兩種形式的排名表。
league.out_remission$table
league.out_remission$longtable
1.5.4 森林圖
nma.forest()函數可輸出森林圖，代碼如下。
f1<- nma.forest(random_effects_results_remission, central.tdcy="median", log.scale =TRUE,comparator = "Steroid alone")
f4<- nma.forest(random_effects_results_remission, central.tdcy="median", log.scale =TRUE,comparator = "Steroid alone")
其中，參數“central.tdcy”和“log.scale”與1.5.3排名表熱圖設置含義相同；參數“comparator”用于設置對照組，通常是某種安慰劑或對照藥物等。
接著，利用如下代碼將圖形輸出，命令執行后結果見圖5。
圖5
nma.forest()函數生成森林圖
grid.arrange(f1, f4, ncol = 2)
2 結果
2.1 網狀關系特征信息及網狀關系圖
BUGSnet程序包的net.tab()函數可輸出網狀、干預和比較特征的統計數據。網狀特征（charactertics）包含了網狀的連通性、網狀中的處理數量及事件數量等一些網狀關系的基本特征。干預特征（intervention）為按治療細分的結果和樣本量數據。比較特征（comparison）為按不同治療比較的結果和樣本量數據。
net.plot()函數可按作者需求輸出美觀的網狀關系圖，各個結局指標的網狀關系圖見圖1，節點大小與干預組納入研究的數量成正比，節點連線的粗細與節點兩端直接比較的研究數量成正比。如果節點之間能形成一個封閉環，代表這些研究能同時參與比較。該結果與示例數據來源的網狀分析圖一致[11]。
2.2 效應模型的選擇與評估
nma.fit()函數可進行模型擬合并識別潛在的異常值，該函數輸出的杠桿圖和偏差信息準則（deviance information criterion，DIC）值可幫助確定最優效應模型。利用該函數生成了各個結局指標的固定效應和隨機效應模型的杠桿圖（leverage plots）。圖中顯示了有效參數數量（pD）、總殘差偏差（Dres）和DIC值等信息（圖2）。我們可利用杠桿圖比較模型的優劣，從而幫助分析者選擇更為合理的效應模型，并且nma.fit()函數能識別出對DIC值變化有影響的數據點，篩選出潛在的異常值。
如圖2所示，杠桿圖表示在所有i次試驗和k個臂中的每一個的杠桿率（leverageik）與貝葉斯偏差殘差（Wik）之間的比較，可在擬合模型時突出潛在的異常值，即如果數據點位于紫色弧線（x2+y=3）之外，則說明該數據點可能導致模型擬合不佳。對于疾病緩解這一組，固定效應模型組有兩個數據點在紫色弧線外。通過查詢R輸出結果發現，該數據點為Beissert 2006的MMF和AZA干預組（表1）。其余各結局指標的固定和隨機效應組均未出現異常值。因此，對于疾病緩解這一組，我們選擇隨機效應模型，或剔除異常值；對于其他各個組固定效應模型和隨機效應模型的選擇采用“模型擬合”越低越好，模型越簡約越好的原則。因此，在累積類固醇劑量組，選擇固定效應模型（DIC更低，模型擬合越好）。
2.3 療效排序結果
nma.rank()可比較各個干預措施處理的后驗概率，實現NMA分析結果的排序比較。該函數中參數“largerbetter”代表所輸出的結果越大是否與更好的治療相關，在這里我們根據研究的意義將疾病緩解設置為“真（TRUE）”，而“累積類固醇劑量”這些值越高風險越大的結局指標設置為“假（FALSE）”，以便直接生成與各個結局指標一致的干預處理分析結果。
如圖3所示，SUCRA和rankogram圖分別以曲線圖和柱狀圖形式直觀地顯示了各個干預措施的排序概率。利用這兩個圖，我們可迅速判斷不同結局指標的最優和最差的干預措施。例如，緩解疾病最有效的措施是“Rituximab”，最差的措施是“DCP_C（6M）”；累積類固醇劑量使用最少的是“Rituximab”，最多的是“Steroid alone”。與示例數據相比，兩個組推薦的最優效果干預措施均與原分析一致[11]，進而確定最優治療方案為“Rituximab”。
2.4 相對效應的估計
nma.league()函數可實現相對效應估計，并能夠比較任何一對干預措施之間的相對效果。如圖4所示的排行榜方陣熱圖，即排名表熱圖，該圖包含所有可能干預對的相對有效性及其95%的可信區間[12]。圖中“Treatment”代表干預措施，“Comparator”代表對照措施。符號（**）表明在95%的可信區間水平下，“Treatment”與其“Comparator”之間存在統計學上的差異。在以下干預措施中，與對照組（Steroid alone）相比，對于緩解疾病最有效的顯著性措施是“Rituximab”，累積類固醇劑量使用最少的措施是“Rituximab”、“AZA”、“CP_P”，這與示例數據結果相似[11]。
2.5 森林圖
nma.forest()函數可輸出森林圖，對藥物在不同研究和不同效應量水平的合并結果進行比較，并顯示了有效的干預措施。如圖5所示，以Steroid alone為參照，將各個結局指標的各個干預措施中的研究進行合并，繪制出森林圖。與之前的排名和排序結果相同的是，通過各組森林圖判斷，“Rituximab”是最有效的疾病緩解干預措施，而“Rituximab”、“AZA”和“CP_P”是類固醇累積使用量最少的干預措施，也是相對安全的治療方式。
3 討論
近年來，隨著越來越多NMA報告的發表，對NMA分析規范及分析手段的需求日益增長。良好的NMA分析流程和規范化的NMA分析結果不僅有助于幫助讀者更好地理解相關療效與風險，也有助于臨床醫生及指南制定者進行臨床證據決策和標準治療指南的制訂[13-15]。
評價具有NMA分析功能軟件的標準取決于其數據處理和圖形展現的能力[16]。相比付費軟件，開源的免費軟件由于其使用靈活及更新速度快，可能更容易滿足普通NMA分析研究人員的要求。貝葉斯算法理論在排序功能上更加優秀，以R語言平臺開發的軟件包的算法多數基于貝葉斯算法理論，另一部分采用頻率算法理論[17]。
當前NMA分析仍然面臨一些挑戰，由于文獻來源、語言、二次分析研究人員等的偏差可能會導致證據決策的誤判[18]。因此，所有NMA分析研究人員都需要按照標準報告規范指南進行資料的收集和數據的分析。BUGSnet程序包是按照ISPOR-AMCP-NPC[19]、PRISMA[20]和NICE-DSU（http://nicedsu.org.uk/）指南標準設計的一款R語言程序包。該程序包在前人的基礎上完善了功能，僅用一款程序包便可實現NMA分析、模型擬合、圖片繪制和表格繪制等。此外，由于是基于R語言平臺，使用友好、無版權困擾。本文以Lee等[11]關于類固醇佐劑治療天皰瘡的NMA數據為例，利用BUGSnet程序包進行NMA分析，得出主要治療結果的證據判斷是與原文完全一致。但由于Lee等[11]采用的是頻率法，本文所展現的BUGSnet程序調用貝葉斯分析程序，理論上當研究數量足夠大時，兩種算法的結果一樣，然而由于示例數據的樣本可能偏少，造成部分次要結果存在一定的差異。
綜上，我們認為BUGSnet程序包不乏是一款準確、實用、方便的NMA分析軟件包。本文提供了清晰的代碼和代碼解釋，以實例的形式演示，方便NMA分析研究人員理解和使用。
近年來，網狀Meta分析（network meta-analysis，NMA）發展迅速，研究者針對不同類型的研究數據開發出對應的解決方案，NMA分析變得越來越簡單，相關發表文獻也呈逐年遞增趨勢[1-2]。Tian等[3]統計了Web of Science數據庫中納入研究NMA方法學文章的引用次數，發現引用次數最高的是Hutton等[4]發表的PRISMA-NMA報告指南，該指南可幫助Meta作者規范撰寫報告并提高報告的完整性。引用次數第2高的是Lu等[5]使用馬爾可夫鏈蒙特卡羅（Markov chain Monte Carlo，MCMC）[6]方法提出的一種分層貝葉斯模型，該模型解決了間接證據之間的比較問題，為NMA貝葉斯分析奠定了基礎。由此可見，NMA作者關心的是NMA分析方法能否按照指南標準對納入研究的數據進行規范和完整性的研究，進而基于證據提供可靠的決策。盡管已經開發出一些NMA分析軟件，但由于存在收費、使用困難、功能不完全、無具體指南標準和可視化效果不良等問題而可能受到使用限制[7-9]。為更好地實施NMA分析報告，Béliveau等[10]創建了一個R程序包BUGSnet（Bayesian inference using gibbs sampling to conduct a network meta-analysis）。該程序包是一款基于R軟件平臺的開源、免費、功能齊備的軟件包，能夠實現貝葉斯NMA分析，并滿足當前Meta分析報告規范的要求。BUGSnet程序包能夠處理不同類型的資料包括計數數據、連續性數據或二分類數據，可用于描述證據網狀、估計模型擬合及評估異質性和不一致性等。此外，該程序包定期更新，并在GitHub中提供使用幫助。本文以實例展示該程序包的具體操作。
1 資料與方法
1.1 數據來源
數據來源于Lee等[11]關于類固醇佐劑治療天皰瘡的NMA研究。該研究納入10個試驗，包含592例患者，評估了7種保留類固醇佐劑的方法治療天皰瘡的效果（表1）。結局指標包含二分類數據（疾病緩解病例數、疾病復發病例數和不良事件相關退出病例數）和連續性數據（類固醇的累積使用劑量）。
表1
納入研究的基本特征[11]
1.2 軟件及程序包的安裝
本文使用的R軟件版本為R x64 4.1.2（下載地址：https://www.r-project.org）。使用R軟件控制臺R Studio操作（下載地址：https://www.rstudio.com/products/rstudio/）。BUGSnet程序包依賴JAGS（Just Another Gibbs Sampler）的廣義線性模型實現貝葉斯NMA，因此需調用JAGS-4.3.0軟件（下載地址：https://sourceforge.net/projects/mcmc-jags/）。上述軟件下載完成后，雙擊軟件安裝包可完成軟件安裝。
完成上述軟件安裝后，打開R Studio，在腳本區輸入以下代碼安裝Rtool。
install.packages("pkgbuild")
pkgbuild::has_build_tools()
接著，在R Studio腳本區輸入以下代碼完成BUGSnet程序包的安裝。
install.packages(c("remotes", "knitr"))
remotes::install_github("audrey-b/BUGSnet@v1.1.0", upgrade = TRUE, build_vignettes = TRUE, dependencies = TRUE)
安裝“xlsx”、“dplyr”和“tidyr”程序包以便于數據的錄入、整理和變換。
install.packages(c("xlsx", "dplyr", "tidyr"))
1.3 數據的錄入和準備
數據按表1所示，將中文漢字改為對應的英文后整理為Excel表格格式，命名為“data_total.xlsx”存儲于演示目錄下。
加載以下程序包準備數據的變換。
library(xlsx)
library(dplyr)
library(tidyr)
輸入以下代碼完成數據錄入。
data_total<-read.xlsx("data_total.xlsx", 1)
輸入以下代碼完成不同結局指標的數據提取。
data_remission<- data_total[, c(1:4)]
data_relapse<-data_total[, c(1:3, 5)]
data_withdrawal<- data_total[, c(1:3, 6)]
data_dose<-data_total[, c(1:3, 7:8)]
輸入以下代碼規范數據錄入格式，即去除缺失值、并將數據存儲到對應的結局指標數據集。
data_remission<-data_remission %>% rename(sampleSize=Number.of.patients, events= Number.of.remissions)
data_relapse<-data_relapse %>% drop_na() %>% rename(sampleSize=Number.of.patients, events=Number.of.relapses)
data_withdrawal<-data_withdrawal %>% drop_na() %>% rename(sampleSize= Number.of.patients, events=Number.of.adverse.event.related.withdrawals)
data_dose<-data_dose %>% drop_na() %>% rename(sampleSize=Number.of.patients, mean=Cumulative.dose.of.steroid_Mean, SD=Cumulative.dose.of.steroid_SD)
在上述代碼中，“data_remission”代表緩解病例數據，為主要療效指標；“data_relapse”代表復發病例數據，為次要療效指標；“data_withdrawal”代表不良事件相關退出病例數據，為次要安全性指標；“data_dose”代表固醇的累積使用劑量，為主要安全性指標，與前三個指標的二分類數據類型不同的是，該指標類型為連續數據。
1.4 BUGSnet程序包數據分析
在完成數據錄入后，正式進入數據分析。如表2所示，BUGSnet程序包的函數可執行包括數據預處理、各種數據分析及圖形或表格結果的輸出等命令。
表2
BUGSnet程序包可實現的功能
在完成數據載入后，可通過BUGSnet程序包將整理好的數據進行NMA分析。在RStudio腳本區輸入以下代碼完成BUGSnet程序包的加載和變量名稱的賦值，以便于進一步分析處理。由于“data_remission”、“data_relapse”和“data_withdrawal”均為二分類數據類型，數據分析方法相似，為避免重復，在后續代碼演示時以“data_remission”為例。
library(BUGSnet)
data_remission <- data.prep(arm.data = data_remission, varname.t = "Treatment", varname.s = "Study")
data_dose <- data.prep(arm.data = data_dose, varname.t = "Treatment", varname.s = "Study")
使用net.tab()函數或net.plot()函數獲得證據網狀的描述性統計表格或圖形。“type.outcome”參數指示數據的類型，可根據結局指標的數據類型進行設置：二分類數據（binomial），連續性數據（continuous）和計數數據（rate/rate2）。
輸入以下代碼獲得證據網狀關系圖（圖1），參數“node.colour”可指定節點的顏色。
圖1
net.plot()函數生成網絡關系圖
a：疾病緩解，b：累積類固醇劑量。
par(mfrow = c(1,2))
net.plot(data_remission, node.colour = "dark grey")
net.plot(data_dose, node.colour = "dark grey")
nma.model()函數用于設置擬合模型，用戶可選擇使用固定效應或隨機效應模型。固定效應模型代碼如下：
fixed_effects_model_remission <- nma.model(data_remission, outcome="events", N="sampleSize", reference="Steroid alone", family="binomial",link="logit", effects="fixed")
fixed_effects_model_dose <- nma.model(data=data_dose, outcome="mean", N="sampleSize", sd="SD", reference="Steroid alone", family="normal", link="identity", effects="fixed")
隨機效應模型代碼如下：
random_effects_model_remission <- nma.model(data_remission, outcome="events", N="sampleSize", reference="Steroid alone", family="binomial", link="logit", effects="random")
random_effects_model_dose <- nma.model(data=data_dose, outcome="mean", N="sampleSize", sd="SD", reference="Steroid alone", family="normal", link="identity", effects="random")
上述命令中，參數“outcome”用于指定結局變量；參數“reference”用于設置對照組，通常是某種安慰劑或對照藥物等；參數“family”用于指定結果的分布類型，通常有“binomial”（二分類數據），“normal”（連續性數據），“poisson”（計數數據）可供選擇；參數“link”用于指定NMA模型使用的函數，其中“logit”用于二分類數據的比值比（Odds Ratio），“log”用于二分類數據的風險比（Risk Ratio）或計數數據的比率（Rate Ratio）；“cloglog”用于二分類數據的危險比（Hazard Ratio）；“identity”用于連續性數據；“effects”用于設置效應模型的類型。
生成NMA模型后，運行nma.run()函數命令執行貝葉斯NMA，該函數運行需要用戶指定MCMC算法的老化（n.burnin）、迭代（n.iter）和適應次數（n.adapt）及監控的變量。
對于固定效應模型的運行代碼如下。
set.seed(20222022)
fixed_effects_results_remission <- nma.run(fixed_effects_model_remission, n.adapt=1000, n.burnin=1000, n.iter=10000)
fixed_effects_results_dose <- nma.run(fixed_effects_model_dose, n.adapt=1000, n.burnin=1000, n.iter=10000)
對于隨機效應模型的運行代碼如下。
random_effects_results_remission <- nma.run(random_effects_model_remission, n.adapt=1000, n.burnin=1000, n.iter=10000)
random_effects_results_dose <- nma.run(random_effects_model_dose, n.adapt=1000, n.burnin=1000, n.iter=10000)
以下代碼可進行NMA分析的效應模型評估，即評估模型的擬合和識別潛在異常值，并輸出杠桿圖（圖2）。
圖2
nma.fit()函數生成的圖和擬合統計數據
leverage
par(mfrow = c(2,2))
對于固定效應模型的擬合代碼如下。
nma.fit(fixed_effects_results_remission)
nma.fit(fixed_effects_results_dose)
對于隨機效應模型的擬合代碼如下。
nma.fit(random_effects_results_remission)
nma.fit(random_effects_results_dose)
1.5 圖形繪制
1.5.1 網狀關系圖
網狀關系圖的繪制如前述，結果見圖1。
1.5.2 治療排名結果圖
使用nma.rank()函數可獲得治療排名的表格和圖形結果。將擬合最優模型結果輸入到該函數中即可獲得對應的排名結果，可參考以下代碼輸出結果，參數“largerbetter”若為“TRUE”則代表更大的結果與更好的治療效果相關，反之亦然。
sucra_out_remission<- nma.rank(random_effects_results_remission, largerbetter=TRUE)
sucra_out_dose<- nma.rank(fixed_effects_results_dose, largerbetter=FALSE)
nma.rank()函數自帶變量顏色，若各個組結局指標的變量名及其數量相同，則可直接在nma.rank()函數中將參數“sucra.palette”設定為所需的顏色集。如果“sucra_out_remission”和“sucra_out_dose”的“Treatment”變量總數及變量名均不完全相同，為統一各個圖片的圖注顏色，可使用以下代碼對“Treatment”中的不同變量設置特定的顏色。指定顏色需用到ggplot2程序包的“scale_colour_manual()”函數，如未安裝ggplot2程序包可使用“install.packages(“ggplot2”)”代碼進行安裝。
ggplot2程序包的加載代碼如下：
library(ggplot2)
各個變量顏色的賦值代碼如下：
z<-scale_colour_manual(values=c(AZA="#FF6600", MMF="#FF0088", "Steroid alone"="#00FF00", CP_P="#007700", Cyclosporine="#6633CC", "DCP_C(6M)"="#00CCFF", "DCP_C(12M)"="#FFFF00", Rituximab="#AA0000"))
繪制累積排序概率圖下面積（surface under the cumulative ranking，SUCRA）的代碼如下：
s1<- sucra_out_remission$sucraplot+z
s4<- sucra_out_dose$sucraplot+z
繪制排序概率圖（rankogram）的代碼如下：
r1<- sucra_out_remission$rankogram+theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
r4<- sucra_out_dose$rankogram+theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1))
輸入以下代碼，完成上述圖片的輸出（圖3）。
圖3
nma.rank()函數生成的SUCRA圖和rankogram圖
library(gridExtra)
grid.arrange(s1, s4, r1, r4, ncol = 2)
上述治療排名結果圖也可表格形式輸出，例如對于疾病緩解數據“sucra_out_remission”，不僅可使用“sucra_out_remission$sucratable”輸出SUCRA圖對應的表格，也可使用“sucra_out_remission$ranktable”輸出排序概率圖對應的表格。
1.5.3 排名表熱圖
nma.league()函數可獲得排名表熱圖，可比較任何一對干預之間的相對效果。將擬合最優模型結果輸入到nma.league()函中，如下面代碼所示：參數“central.tdcy”是用來衡量相對有效性的統計數據，可選擇“mean”或“median”；參數“order”表示圖片顯示的治療順序，本例以1.5.2治療排名結果圖中的輸出結果來排序；參數“log.scale”代表數據呈現方式，如果為“TRUE”則表示以對數刻度報告結果；參數“low.colour”、“mid.colour”及“high.colour”可設置熱圖的低、中、高值的顏色；參數“digits”設定小數點后的顯示位數，相關代碼設置如下。
league.out_remission<- nma.league(random_effects_results_remission, central.tdcy="median", order = sucra_out_remission$order, log.scale =TRUE, low.colour = "springgreen4", mid.colour = "white", high.colour = "red", digits = 2)
league.out_dose<- nma.league(fixed_effects_results_dose, central.tdcy="median", order = sucra_out_dose$order, log.scale =TRUE, low.colour = "springgreen4", mid.colour = "white", high.colour = "red", digits = 2)
使用以下代碼將排名表熱圖輸出結果放置在同一頁，輸出結果見圖4。
圖4
nma.league()函數生成排名表熱圖
l1<-league.out_remission$heatplot
l4<-league.out_dose$heatplot
grid.arrange(l1, l4, ncol = 1)
排名表熱圖也可表格的形式輸出，以疾病緩解數據為例，使用以下代碼可輸出兩種形式的排名表。
league.out_remission$table
league.out_remission$longtable
1.5.4 森林圖
nma.forest()函數可輸出森林圖，代碼如下。
f1<- nma.forest(random_effects_results_remission, central.tdcy="median", log.scale =TRUE,comparator = "Steroid alone")
f4<- nma.forest(random_effects_results_remission, central.tdcy="median", log.scale =TRUE,comparator = "Steroid alone")
其中，參數“central.tdcy”和“log.scale”與1.5.3排名表熱圖設置含義相同；參數“comparator”用于設置對照組，通常是某種安慰劑或對照藥物等。
接著，利用如下代碼將圖形輸出，命令執行后結果見圖5。
圖5
nma.forest()函數生成森林圖
grid.arrange(f1, f4, ncol = 2)
2 結果
2.1 網狀關系特征信息及網狀關系圖
BUGSnet程序包的net.tab()函數可輸出網狀、干預和比較特征的統計數據。網狀特征（charactertics）包含了網狀的連通性、網狀中的處理數量及事件數量等一些網狀關系的基本特征。干預特征（intervention）為按治療細分的結果和樣本量數據。比較特征（comparison）為按不同治療比較的結果和樣本量數據。
net.plot()函數可按作者需求輸出美觀的網狀關系圖，各個結局指標的網狀關系圖見圖1，節點大小與干預組納入研究的數量成正比，節點連線的粗細與節點兩端直接比較的研究數量成正比。如果節點之間能形成一個封閉環，代表這些研究能同時參與比較。該結果與示例數據來源的網狀分析圖一致[11]。
2.2 效應模型的選擇與評估
nma.fit()函數可進行模型擬合并識別潛在的異常值，該函數輸出的杠桿圖和偏差信息準則（deviance information criterion，DIC）值可幫助確定最優效應模型。利用該函數生成了各個結局指標的固定效應和隨機效應模型的杠桿圖（leverage plots）。圖中顯示了有效參數數量（pD）、總殘差偏差（Dres）和DIC值等信息（圖2）。我們可利用杠桿圖比較模型的優劣，從而幫助分析者選擇更為合理的效應模型，并且nma.fit()函數能識別出對DIC值變化有影響的數據點，篩選出潛在的異常值。
如圖2所示，杠桿圖表示在所有i次試驗和k個臂中的每一個的杠桿率（leverageik）與貝葉斯偏差殘差（Wik）之間的比較，可在擬合模型時突出潛在的異常值，即如果數據點位于紫色弧線（x2+y=3）之外，則說明該數據點可能導致模型擬合不佳。對于疾病緩解這一組，固定效應模型組有兩個數據點在紫色弧線外。通過查詢R輸出結果發現，該數據點為Beissert 2006的MMF和AZA干預組（表1）。其余各結局指標的固定和隨機效應組均未出現異常值。因此，對于疾病緩解這一組，我們選擇隨機效應模型，或剔除異常值；對于其他各個組固定效應模型和隨機效應模型的選擇采用“模型擬合”越低越好，模型越簡約越好的原則。因此，在累積類固醇劑量組，選擇固定效應模型（DIC更低，模型擬合越好）。
2.3 療效排序結果
nma.rank()可比較各個干預措施處理的后驗概率，實現NMA分析結果的排序比較。該函數中參數“largerbetter”代表所輸出的結果越大是否與更好的治療相關，在這里我們根據研究的意義將疾病緩解設置為“真（TRUE）”，而“累積類固醇劑量”這些值越高風險越大的結局指標設置為“假（FALSE）”，以便直接生成與各個結局指標一致的干預處理分析結果。
如圖3所示，SUCRA和rankogram圖分別以曲線圖和柱狀圖形式直觀地顯示了各個干預措施的排序概率。利用這兩個圖，我們可迅速判斷不同結局指標的最優和最差的干預措施。例如，緩解疾病最有效的措施是“Rituximab”，最差的措施是“DCP_C（6M）”；累積類固醇劑量使用最少的是“Rituximab”，最多的是“Steroid alone”。與示例數據相比，兩個組推薦的最優效果干預措施均與原分析一致[11]，進而確定最優治療方案為“Rituximab”。
2.4 相對效應的估計
nma.league()函數可實現相對效應估計，并能夠比較任何一對干預措施之間的相對效果。如圖4所示的排行榜方陣熱圖，即排名表熱圖，該圖包含所有可能干預對的相對有效性及其95%的可信區間[12]。圖中“Treatment”代表干預措施，“Comparator”代表對照措施。符號（**）表明在95%的可信區間水平下，“Treatment”與其“Comparator”之間存在統計學上的差異。在以下干預措施中，與對照組（Steroid alone）相比，對于緩解疾病最有效的顯著性措施是“Rituximab”，累積類固醇劑量使用最少的措施是“Rituximab”、“AZA”、“CP_P”，這與示例數據結果相似[11]。
2.5 森林圖
nma.forest()函數可輸出森林圖，對藥物在不同研究和不同效應量水平的合并結果進行比較，并顯示了有效的干預措施。如圖5所示，以Steroid alone為參照，將各個結局指標的各個干預措施中的研究進行合并，繪制出森林圖。與之前的排名和排序結果相同的是，通過各組森林圖判斷，“Rituximab”是最有效的疾病緩解干預措施，而“Rituximab”、“AZA”和“CP_P”是類固醇累積使用量最少的干預措施，也是相對安全的治療方式。
3 討論
近年來，隨著越來越多NMA報告的發表，對NMA分析規范及分析手段的需求日益增長。良好的NMA分析流程和規范化的NMA分析結果不僅有助于幫助讀者更好地理解相關療效與風險，也有助于臨床醫生及指南制定者進行臨床證據決策和標準治療指南的制訂[13-15]。
評價具有NMA分析功能軟件的標準取決于其數據處理和圖形展現的能力[16]。相比付費軟件，開源的免費軟件由于其使用靈活及更新速度快，可能更容易滿足普通NMA分析研究人員的要求。貝葉斯算法理論在排序功能上更加優秀，以R語言平臺開發的軟件包的算法多數基于貝葉斯算法理論，另一部分采用頻率算法理論[17]。
當前NMA分析仍然面臨一些挑戰，由于文獻來源、語言、二次分析研究人員等的偏差可能會導致證據決策的誤判[18]。因此，所有NMA分析研究人員都需要按照標準報告規范指南進行資料的收集和數據的分析。BUGSnet程序包是按照ISPOR-AMCP-NPC[19]、PRISMA[20]和NICE-DSU（http://nicedsu.org.uk/）指南標準設計的一款R語言程序包。該程序包在前人的基礎上完善了功能，僅用一款程序包便可實現NMA分析、模型擬合、圖片繪制和表格繪制等。此外，由于是基于R語言平臺，使用友好、無版權困擾。本文以Lee等[11]關于類固醇佐劑治療天皰瘡的NMA數據為例，利用BUGSnet程序包進行NMA分析，得出主要治療結果的證據判斷是與原文完全一致。但由于Lee等[11]采用的是頻率法，本文所展現的BUGSnet程序調用貝葉斯分析程序，理論上當研究數量足夠大時，兩種算法的結果一樣，然而由于示例數據的樣本可能偏少，造成部分次要結果存在一定的差異。
綜上，我們認為BUGSnet程序包不乏是一款準確、實用、方便的NMA分析軟件包。本文提供了清晰的代碼和代碼解釋，以實例的形式演示，方便NMA分析研究人員理解和使用。
