NMA for Binary data (2-arm trials)

Introduction

This vignette provides a short example of a Bayesian NMA for Binary data. The model fit relies on the gemtc package, pre- and post-processing is done with gemtcPlus.

Prepare the environment

library(dplyr)
library(gemtc)        
library(gemtcPlus)     
library(ggmcmc)

Load in the data

data("binary_data", package = "gemtcPlus") # This should be binary

Plan the model

model_plan <- plan_binary(bth.model = "RE", 
                          n.chain = 3, 
                          n.iter = 6000, 
                          thin = 1,
                          n.adapt = 1000, 
                          link = "logit",
                          bth.prior =  mtc.hy.prior(type = "var",
                                                   distr = "dlnorm",-4.18, 1 / 1.41 ^ 2)
                          )

Ready the data

model_input <- nma_pre_proc(binary_data, plan = model_plan)

Figure Network plot

plot(model_input$fitting_data)

Fit the model

model  <- nma_fit(model_input = model_input)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 22
##    Total graph size: 400
## 
## Initializing model

Post processing

Inspect convergence

The ggmcmc package provides ggplot2 versions of all major convergence plots and diagnostics.

Figure Traceplot

ggs_traceplot(ggs(model$samples))

Figure Densityplot

ggs_density(ggs(model$samples))

Many more diagnostic plots are available, such as Brooks-Gelman-Rubin convergence diagnostic (Rhat), auto-correlation plot, or running means.

ggs_Rhat(ggs(model$samples))

ggs_autocorrelation(ggs(model$samples))

ggs_running(ggs(model$samples))

Produce outputs of interest

Posterior summaries (log-scale)

The contrasts in this model are log-odds ratios.

Unfortunately, gemtc does not provide an estimate of the effective sample size (n.eff). Instead, a time-series SE is given. As a rule of thumb, the length of the MCMC is sufficient if the time-series SE is smaller than 2%(-5%) of the posterior SD.

summary(model)
## 
## Results on the Log Odds Ratio scale
## 
## Iterations = 1001:7000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 6000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean     SD Naive SE Time-series SE
## d.B.A 0.9773 0.2995 0.002232       0.003498
## d.B.D 1.4322 0.3948 0.002943       0.012389
## d.B.E 0.9654 0.4820 0.003592       0.019097
## d.E.C 0.7142 0.3439 0.002563       0.005321
## d.E.F 0.4904 0.4369 0.003257       0.006454
## sd.d  0.3134 0.1982 0.001477       0.009943
## 
## 2. Quantiles for each variable:
## 
##           2.5%    25%    50%    75%  97.5%
## d.B.A  0.37256 0.8096 0.9747 1.1418 1.6031
## d.B.D  0.50365 1.2225 1.4796 1.6877 2.0876
## d.B.E -0.14996 0.6977 1.0288 1.2942 1.7559
## d.E.C  0.05764 0.5108 0.6990 0.8975 1.4686
## d.E.F -0.39221 0.2435 0.4897 0.7429 1.3653
## sd.d   0.04960 0.1629 0.2747 0.4232 0.7962
## 
## -- Model fit (residual deviance):
## 
##     Dbar       pD      DIC 
## 21.72859 15.76110 37.48968 
## 
## 16 data points, ratio 1.358, I^2 = 31%

To judge overall model fit, the residual deviance should be compared to the number of independent data points (which can be done via a small utility function in gemtcPlus).

get_mtc_sum(model)
##     DIC    pD resDev dataPoints
## 1 37.49 15.76  21.73         16

Odds ratio (OR) estimates

Assume new treatment is “A” and is to be compared vs all other treatments.

Table Odds ratios treatment A vs other treatments

OR <- get_mtc_newVsAll(model, new.lab = "A", transform = "exp", digits = 2)
OR
##   Comparator  Med CIlo CIup
## 1          B 2.65 1.45 4.97
## 2          C 0.46 0.19 1.92
## 3          D 0.61 0.27 2.00
## 4          E 0.94 0.40 3.76
## 5          F 0.58 0.19 3.17

Table Probability A better than other treatments (better meaning larger OR)

get_mtc_probBetter(model, new.lab = "A", smaller.is.better = FALSE, sort.by = "effect")
##   New Comparator probNewBetter
## 1   A          B         0.996
## 4   A          E         0.452
## 5   A          F         0.208
## 3   A          D         0.150
## 2   A          C         0.103

Figure Forest plot A vs other treatments

plot_mtc_forest(x = OR, 
                lab = "Odds ratio A vs others",
                breaks = c(0.125, 0.25, 0.5, 0.8, 1, 1.25, 2, 4, 8, 12), 
                sort.by = "effect")

Table Cross-tabulation of ORs

ctab <- round(exp(relative.effect.table(model)), 2)
pander::pandoc.table(as.data.frame(ctab), split.tables = Inf)
  A B C D E F
A A 0.38 (0.2, 0.69) 2.15 (0.52, 5.37) 1.65 (0.5, 3.65) 1.06 (0.27, 2.52) 1.72 (0.32, 5.37)
B 2.65 (1.45, 4.97) B 5.7 (1.71, 12.39) 4.39 (1.65, 8.07) 2.8 (0.86, 5.79) 4.55 (0.95, 12.93)
C 0.46 (0.19, 1.92) 0.18 (0.08, 0.59) C 0.77 (0.42, 1.57) 0.5 (0.23, 0.94) 0.82 (0.24, 2.35)
D 0.61 (0.27, 2) 0.23 (0.12, 0.6) 1.3 (0.64, 2.41) D 0.64 (0.3, 1.14) 1.05 (0.31, 2.87)
E 0.94 (0.4, 3.76) 0.36 (0.17, 1.16) 2.01 (1.06, 4.34) 1.56 (0.88, 3.39) E 1.63 (0.68, 3.92)
F 0.58 (0.19, 3.17) 0.22 (0.08, 1.05) 1.23 (0.42, 4.12) 0.95 (0.35, 3.24) 0.61 (0.26, 1.48) F

Treatment rankings

rk  <- rank.probability(model, preferredDirection = 1)
mrk <- reshape2::melt(rk[,], varnames = c("Treatment", "Rank"), value.name = "Probability")

fig <- ggplot(data = mrk) +
  geom_line(aes(Rank, Probability, color = Treatment, linetype = Treatment), size = 1.5) +
  theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Figure Rankogram

plot(fig)

Table Rank probabilities

colnames(rk) <- paste("Rank", 1:ncol(rk))
rk
## Rank probability; preferred direction = 1
##         Rank 1      Rank 2      Rank 3     Rank 4      Rank 5      Rank 6
## A 6.188889e-02 0.065222222 0.124555556 0.21911111 0.526333333 0.002888889
## B 5.555556e-05 0.001166667 0.002666667 0.01533333 0.037222222 0.943555556
## C 5.616111e-01 0.308055556 0.089555556 0.03038889 0.008166667 0.002222222
## D 8.822222e-02 0.361611111 0.446111111 0.08711111 0.015777778 0.001166667
## E 9.444444e-04 0.009833333 0.081500000 0.52188889 0.356222222 0.029611111
## F 2.872778e-01 0.254111111 0.255611111 0.12616667 0.056277778 0.020555556

Node-splitting (inconsistency assessment)

nsplit <- mtc.nodesplit(model$model$network)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 23
##    Total graph size: 509
## 
## Initializing model
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 23
##    Total graph size: 510
## 
## Initializing model
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 23
##    Total graph size: 509
## 
## Initializing model
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 23
##    Total graph size: 509
## 
## Initializing model
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 23
##    Total graph size: 509
## 
## Initializing model
## 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 16
##    Unobserved stochastic nodes: 22
##    Total graph size: 397
## 
## Initializing model
summary(nsplit)
## Node-splitting analysis of inconsistency
## ========================================
## 
##    comparison  p.value  CrI               
## 1  d.B.D       0.075775                   
## 2  -> direct            1.9 (0.32, 3.4)   
## 3  -> indirect          -0.41 (-2.7, 1.8) 
## 4  -> network           1.2 (-0.70, 2.8)  
## 5  d.B.E       0.073600                   
## 6  -> direct            -0.60 (-2.5, 1.2) 
## 7  -> indirect          1.7 (-0.28, 3.8)  
## 8  -> network           0.58 (-1.5, 2.2)  
## 9  d.C.D       0.964350                   
## 10 -> direct            -0.21 (-2.6, 2.2) 
## 11 -> indirect          -0.14 (-3.3, 3.1) 
## 12 -> network           -0.20 (-1.9, 1.6) 
## 13 d.C.E       0.960175                   
## 14 -> direct            -0.77 (-3.2, 1.7) 
## 15 -> indirect          -0.85 (-4.1, 2.2) 
## 16 -> network           -0.78 (-2.6, 0.88)
## 17 d.D.E       0.230750                   
## 18 -> direct            0.098 (-2.0, 2.2) 
## 19 -> indirect          -1.3 (-3.6, 0.71) 
## 20 -> network           -0.59 (-2.2, 0.89)
plot(summary(nsplit))

Forest plot

HR_i    <- get_mtc_newVsAll(model, new.lab = "A", transform = "exp", digits = 2)
plot_mtc_forest(HR_i)
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_segment()`).

Extract model code (e.g. for Appendix)

cat(model$code)

Session info

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggmcmc_1.5.1.1  ggplot2_3.5.1   tidyr_1.3.1     gemtcPlus_1.0.0
##  [5] R2jags_0.8-5    rjags_4-16      gemtc_1.0-2     coda_0.19-4.1  
##  [9] dplyr_1.1.4     rmarkdown_2.28 
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.5          xfun_0.48             bslib_0.8.0          
##  [4] GGally_2.2.1          CompQuadForm_1.4.3    lattice_0.22-6       
##  [7] tzdb_0.4.0            numDeriv_2016.8-1.1   mathjaxr_1.6-0       
## [10] vctrs_0.6.5           tools_4.4.1           generics_0.1.3       
## [13] parallel_4.4.1        tibble_3.2.1          fansi_1.0.6          
## [16] highr_0.11            pkgconfig_2.0.3       Matrix_1.7-0         
## [19] RColorBrewer_1.1-3    truncnorm_1.0-9       lifecycle_1.0.4      
## [22] farver_2.1.2          compiler_4.4.1        stringr_1.5.1        
## [25] munsell_0.5.1         htmltools_0.5.8.1     sys_3.4.3            
## [28] buildtools_1.0.0      sass_0.4.9            yaml_2.3.10          
## [31] crayon_1.5.3          pillar_1.9.0          nloptr_2.1.1         
## [34] jquerylib_0.1.4       MASS_7.3-61           cachem_1.1.0         
## [37] meta_7.0-0            metadat_1.2-0         boot_1.3-31          
## [40] abind_1.4-8           nlme_3.1-166          ggstats_0.7.0        
## [43] network_1.18.2        tidyselect_1.2.1      digest_0.6.37        
## [46] slam_0.1-53           stringi_1.8.4         reshape2_1.4.4       
## [49] pander_0.6.5          purrr_1.0.2           labeling_0.4.3       
## [52] maketools_1.3.1       forcats_1.0.0         splines_4.4.1        
## [55] fastmap_1.2.0         grid_4.4.1            colorspace_2.1-1     
## [58] cli_3.6.3             metafor_4.6-0         magrittr_2.0.3       
## [61] utf8_1.2.4            readr_2.1.5           withr_3.0.1          
## [64] scales_1.3.0          igraph_2.0.3          lme4_1.1-35.5        
## [67] hms_1.1.3             evaluate_1.0.1        knitr_1.48           
## [70] Rglpk_0.6-5.1         rlang_1.1.4           Rcpp_1.0.13          
## [73] glue_1.8.0            xml2_1.3.6            minqa_1.2.8          
## [76] jsonlite_1.8.9        R6_2.5.1              plyr_1.8.9           
## [79] statnet.common_4.10.0 R2WinBUGS_2.1-22.1

```