本文主要是使用BUGSnet程序包对2007年发表在Lancet上的《Incident diabetes in clinical trials of antihypertensive drugs: a network meta-analysis》一文结果进行复现。
公众号回:diabetes 获取效应数据。
1. 程序包加载及数据导入
install.packages(c("remotes", "knitr"))
remotes::install_github("audrey-b/BUGSnet@v1.1.0", upgrade = TRUE,
build_vignettes = TRUE, dependencies = TRUE)
library(BUGSnet)
library(readr)
data <- read_csv("diabetes.csv")
2. 干预措施信息比较
在进行分析时首先使用data.prep
函数对数据进行规整,主要用于干预措施(treatment)和研究(study)。
data <- data.prep(arm.data = data,
varname.t = "treatment",
varname.s = "study")
BUGSnet
程序包的 net.tab()
函数可输出网状、干预和比较特征的统计数据,其中:
-
network 输出的结果包含了网状的连通性、网状中的处理数量及事件数量等一些网状关系的基本特征
-
intervention 输出的结果包含了不同干预措施的总发生事件数,总样本量等信息
-
comparison 输出的结果包含了按不同治疗比较的事件数和总样本量等信息
network.char <- net.tab(data = data,
outcome = "responders",
N = "sampleSize",
type.outcome = "binomial",
time = NULL)
> network.char$network
# A tibble: 12 × 2
Characteristic Value
<chr> <chr>
1 Number of Interventions 6
2 Number of Studies 22
3 Total Number of Patients in Network 154176
4 Total Possible Pairwise Comparisons 15
5 Total Number of Pairwise Comparisons With Direct Data 14
6 Is the network connected? TRUE
7 Number of Two-arm Studies 18
8 Number of Multi-Arms Studies 4
9 Total Number of Events in Network 10962
10 Number of Studies With No Zero Events 22
11 Number of Studies With At Least One Zero Event 0
12 Number of Studies with All Zero Events 0
> network.char$intervention
# A tibble: 6 × 7
treatment n.studies n.events n.patients min.outcome max.outcome av.outcome
<chr> <int> <int> <int> <dbl> <dbl> <dbl>
1 ACE 8 1618 23351 0.0291 0.171 0.0693
2 ARB 5 1189 14185 0.00510 0.136 0.0838
3 BBlocker 9 2705 36150 0.0261 0.173 0.0748
4 CCB 9 2791 38809 0.0366 0.167 0.0719
5 Diuretic 8 973 18699 0.0229 0.0858 0.0520
6 Placebo 9 1686 22982 0.0154 0.185 0.0734
> network.char$comparison
# A tibble: 14 × 5
comparison n.studies n.patients n.outcomes proportion
<chr> <int> <int> <int> <dbl>
1 ACE vs. BBlocker 3 15158 1022 0.0674
2 ACE vs. CCB 3 12597 538 0.0427
3 ACE vs. Diuretic 2 16488 759 0.0460
4 ACE vs. Placebo 3 17893 1929 0.108
5 ARB vs. BBlocker 1 7999 562 0.0703
6 ARB vs. CCB 1 10161 1535 0.151
7 ARB vs. Diuretic 1 392 9 0.0230
8 ARB vs. Placebo 2 9778 573 0.0586
9 BBlocker vs. CCB 5 44974 3361 0.0747
10 BBlocker vs. Diuretic 2 8752 241 0.0275
11 BBlocker vs. Placebo 1 3315 71 0.0214
12 CCB vs. Diuretic 2 15739 768 0.0488
13 CCB vs. Placebo 1 9711 331 0.0341
14 Diuretic vs. Placebo 3 7343 384 0.0523
3. 网络证据图( Network plot)
网络证据图使用net.plot()
函数进行绘制,不需进行过多设置,便可以绘制比较好的图片。
net.plot(data)
也可以进行复杂的图形绘制:
net.plot(data,
node.scale=4,
edge.scale=1.5,
flag="Placebo",
label.offset1=0.5,
label.offset2=0.5,
study.counts=T,
graph.scale = TRUE,
node.lab.cex = 1,
edge.lab.cex = 1,
node.colour = "steelblue",
edge.colour = "black",
edge.lab.colour = "blue",
flag.edge.colour = "black",
layout = "layout_in_circle",
layout.params = NULL)
4. 模型选择(固定vs随机,一致vs不一致)
nma.model()
函数是该程序包的核心函数,用来构建NMA模型可以运行固定效应模型、随机效应模型及一致性模型、不一致性模型;并提供了模型比较图,来进行模型选择。
4.1 固定效应模型和随机效应模型选择
fixmodel <- nma.model(data=data,#固定效应模型
outcome = "responders",
N = "sampleSize",
reference="Placebo",
family="binomial",
link="log",
effects="fixed")
ranmodel <- nma.model(data = data,#随机效应模型
outcome = "responders",
N = "sampleSize",
reference="Placebo",
family="binomial",
link="log",
effects="random")
fixresults <- nma.run(fixmodel,
n.adapt=1000,
n.burnin=10000,
n.iter=50000)
ranresults <- nma.run(ranmodel,
n.adapt=1000,
n.burnin=10000,
n.iter=50000)
上述代码中,
-
参数“outcome”用于指定结局变量;
-
参数reference用于设置对照组,通常是某种安慰剂或对照药物等;
-
参数family用于指定结果的分布类型,通常有“binomial”(二分类数据),“normal”(连续性数据),“poisson”(计数数据)可供选择;
-
参数 link用于指 定NMA 模型使用的函数,其中“logit”用于二分类数据的比值比(Odds Ratio),“log”用于二分类数据的风险比(Risk Ratio)或计数数据的比率(RateRatio);“cloglog”用于二分类数据的危险比(Hazard Ratio);identity用于连续性数据;
-
effects用于设置效应模型的类型。
使用图形展示模型结果:
par(mfrow=c(1,2))
nma.fit(fixresults, main ="fixed model" )
nma.fit(ranresults, main= "random model")
可使用nma.diag()
函数绘制模型诊断图:
nma.diag(ranresults)
4.2 一致性模型和不一性模型选择
根据DIC最小原则,上述模型结果显示随机效应模型拟合更佳,可进一步选择一致和不一致模型:
inconmodel <- nma.model(data = data,
outcome = "responders",
N = "sampleSize",
reference="Placebo",
family="binomial",
link="log",
type = "inconsistency",
effects="random")
inconresults<- nma.run(inconmodel,
n.adapt=5000,
n.burnin=3,
n.iter=20000)
conresults<- nma.model(data = data,
outcome = "responders",
N = "sampleSize",
reference="Placebo",
family = "binomial",
link = "log",
effects = "random",
type="consistency",)
conresults <- nma.run(conresults,
n.adapt=5000,
n.burnin=3,
n.iter=20000)
对结果进行可视化:
par(mfrow=c(1,2))
nma.fit(conresults,main = "consistency" )
nma.fit(inconresults,main = "inconsistency")
根据DIC最小原则,最终选择随机效应模型和一致性模型。
5. 森林图 (Forest plot)
5.1 网状meta分析森林图
使用nma.forest()
函数输出森林图,对药物在不同研究和不同效应量水平的合并结果进行比较,并显示了有效的干预措施。
nma.forest(nma = ranresults,
log.scale=FALSE,
central.tdcy="mean",
comparator = "Placebo")
5.2 成对比较森林图
除了绘制网状meta结果的森林图外,该程序包还可以调用metafor程序包使用pma
对两两干预措施进行直接比较,并绘制森林图,具体代码如下:
pma(data = data,
type.outcome="binomial",
sm="OR",
name.trt1 = "CCB",
name.trt2 = "BBlocker",
outcome = "responders",
N = "sampleSize")
$summary
comparison n.studies i.squared re.estimate re.lci re.uci fe.estimate fe.lci fe.uci
1 CCB vs. BBlocker 5 0.6249776 1.230353 1.075073 1.408062 1.277449 1.189803 1.371551
$forest
$forest$xlim
[1] -0.5615787 0.5615787
$forest$addrows.below.overall
[1] 0
$forest$colgap
[1] 2mm
$forest$colgap.left
[1] 2mm
$forest$colgap.right
[1] 2mm
$forest$colgap.studlab
[1] 2mm
$forest$colgap.forest
[1] 2mm
$forest$colgap.forest.left
[1] 2mm
$forest$colgap.forest.right
[1] 2mm
$forest$studlab
[1] "1" "5" "14" "17" "21"
$forest$TE.format
[1] "" "" "" "0.10" "0.38" "0.17" "0.16" "0.02"
$forest$seTE.format
[1] "" "" "" "0.2332" "0.0577" "0.0594" "0.0950" "0.1480"
$forest$cluster.format
[1] "" "" ""
$forest$effect.format
[1] "1.28" "1.23" "" "1.11" "1.47" "1.19" "1.18" "1.02"
$forest$ci.format
[1] "[1.19; 1.37]" "[1.08; 1.41]" "" "[0.70; 1.75]" "[1.31; 1.64]" "[1.06; 1.33]"
[7] "[0.98; 1.42]" "[0.77; 1.37]"
$raw
Number of studies combined: k = 5
Number of observations: o = 44974
Number of events: e = 3361
OR 95%-CI z p-value
Common effect model 1.2774 [1.1898; 1.3716] 6.75 < 0.0001
Random effects model 1.2304 [1.0751; 1.4081] 3.01 0.0026
Quantifying heterogeneity:
tau^2 = 0.0130; tau = 0.1138; I^2 = 62.5% [0.6%; 85.8%]; H = 1.63 [1.00; 2.66]
Test of heterogeneity:
Q d.f. p-value
10.67 4 0.0306
Details on meta-analytical method:
- Mantel-Haenszel method
- DerSimonian-Laird estimator for tau^2
- Mantel-Haenszel estimator used in calculation of Q and tau^2 (like RevMan 5)
6. 累计疗效排序图(SUCRA plot)
使用nma.rank()
函数进行干预措施的排序,结果可产生rankogram
、longtable
和sucraplot
三种形式。
sucra.out <- nma.rank(ranresults,
largerbetter=T,
sucra.palette= "Set1")
sucra.out$rankogram
sucra.out$sucraplot
league.out$table
write.csv(league.out$table,"effectR.csv"
7. 赛联表及赛联热图( league heat plot)
两两比较使用nma.league()
函数可实现,结果可产出赛联表和赛联热图。
league.out <- nma.league(ranresults,
central.tdcy="mean",
high.colour = "cornflowerblue")
league.out$heatplot
也可以产出赛联表,并导出EXCEL。
league.out$table
write.csv(league.out$table,"effectR.csv")
参考文献
1. Béliveau A, Boyne DJ, Slater J, Brenner D, Arora P. BUGSnet: an R package to facilitate the conduct and reporting of Bayesian network Meta-analyses. BMC Med Res Methodol. 2019;19(1):196. Published 2019 Oct 22. doi:10.1186/s12874-019-0829-2.
2. .https://bugsnetsoftware.github.io/instructions.html