Bayesian FE & RE NMA for HR data (via gemtc package) - Result Generation

Introduction

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

Prepare the environment

knitr::opts_chunk$set(echo = TRUE)

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

Load in the data

# load example data
data("hr_data", package = "gemtcPlus")

Plan the model

#Plan model
model_plan <- plan_hr(bth.model = "FE", 
                          n.chain = 3, 
                          n.iter = 6000, 
                          thin = 1,
                          n.adapt = 1000, 
                          link = "identity",
                          linearModel = "fixed")

Ready the data

Fit the model

model  <- nma_fit(model_input = model_input)
## Warning in (function (network, type = "consistency", factor = 2.5, n.chain = 4,
## : Likelihood can not be inferred. Defaulting to normal.
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 5
##    Total graph size: 132
## 
## Initializing model

Post processing

ggs_traceplot(ggs(model$samples))

ggs_density(ggs(model$samples))

summary(model)
## 
## Results on the Mean Difference scale
## 
## Iterations = 1:6000
## 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.12835 0.05539 0.0004129      0.0004129
## d.C.B  0.19726 0.08320 0.0006201      0.0006313
## d.C.D  0.03578 0.12182 0.0009080      0.0008890
## d.C.E -0.09242 0.06593 0.0004914      0.0004910
## d.D.F  0.20852 0.19508 0.0014540      0.0014541
## 
## 2. Quantiles for each variable:
## 
##           2.5%      25%      50%      75%    97.5%
## d.B.A -0.23636 -0.16594 -0.12837 -0.09106 -0.01975
## d.C.B  0.03411  0.14194  0.19794  0.25263  0.36011
## d.C.D -0.20209 -0.04794  0.03582  0.11856  0.27699
## d.C.E -0.22187 -0.13704 -0.09198 -0.04728  0.03591
## d.D.F -0.17147  0.07815  0.20894  0.33972  0.59366
## 
## -- Model fit (residual deviance):
## 
##      Dbar        pD       DIC 
##  5.481882  4.980620 10.462501 
## 
## 7 data points, ratio 0.7831, I^2 = 0%
plot(model_input$fitting_data)

get_mtc_sum(model)
##     DIC   pD resDev dataPoints
## 1 10.46 4.98   5.48          7

Random Effects example

Update model plan and re-run fit.

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

Fit the model

model  <- nma_fit(model_input = model_input)
## Warning in (function (network, type = "consistency", factor = 2.5, n.chain = 4,
## : Likelihood can not be inferred. Defaulting to normal.
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 7
##    Unobserved stochastic nodes: 5
##    Total graph size: 132
## 
## 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))

Figure Brooks-Gelman-Rubin convergence diagnostic (Rhat)

ggs_Rhat(ggs(model$samples))

Figure Auto-correlation plot

ggs_autocorrelation(ggs(model$samples))

Figure Running means

ggs_running(ggs(model$samples))

Produce outputs of interest

Posterior summaries (log-scale)

The contrasts in this model are log-hazard ratios (which correspond to differences in log-hazard rates).

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 Mean Difference scale
## 
## Iterations = 1:6000
## 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.12835 0.05539 0.0004129      0.0004129
## d.C.B  0.19726 0.08320 0.0006201      0.0006313
## d.C.D  0.03578 0.12182 0.0009080      0.0008890
## d.C.E -0.09242 0.06593 0.0004914      0.0004910
## d.D.F  0.20852 0.19508 0.0014540      0.0014541
## 
## 2. Quantiles for each variable:
## 
##           2.5%      25%      50%      75%    97.5%
## d.B.A -0.23636 -0.16594 -0.12837 -0.09106 -0.01975
## d.C.B  0.03411  0.14194  0.19794  0.25263  0.36011
## d.C.D -0.20209 -0.04794  0.03582  0.11856  0.27699
## d.C.E -0.22187 -0.13704 -0.09198 -0.04728  0.03591
## d.D.F -0.17147  0.07815  0.20894  0.33972  0.59366
## 
## -- Model fit (residual deviance):
## 
##      Dbar        pD       DIC 
##  5.481882  4.980620 10.462501 
## 
## 7 data points, ratio 0.7831, I^2 = 0%

In the example here, the chain length seems borderline (sufficient for posterior means and medians, but rather a bit too small for stable 95% credible intervals).

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 10.46 4.98   5.48          7

Hazard ratio estimates

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

Table Hazard ratios A vs other treatments

HR <- get_mtc_newVsAll(model, new.lab = "A", transform = "exp", digits = 2)
HR
##   Comparator  Med CIlo CIup
## 1          B 0.88 0.79 0.98
## 2          C 1.07 0.88 1.30
## 3          D 1.03 0.76 1.41
## 4          E 1.17 0.93 1.48
## 5          F 0.84 0.51 1.37

Table Probability A better than other treatments (better meaning smaller HR)

get_mtc_probBetter(model, new.lab = "A", smaller.is.better = TRUE, sort.by = "effect")
##   New Comparator probNewBetter
## 1   A          B         0.989
## 5   A          F         0.756
## 3   A          D         0.415
## 2   A          C         0.244
## 4   A          E         0.090

Figure Forest plot A vs other treatments

plot_mtc_forest(x = HR, lab = "Hazard ratio A vs others", sort.by = "effect")

Table Cross-tabulation of HRs

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 1.14 (1.02, 1.27) 0.93 (0.77, 1.14) 0.97 (0.71, 1.32) 0.85 (0.67, 1.08) 1.19 (0.73, 1.94)
B 0.88 (0.79, 0.98) B 0.82 (0.7, 0.97) 0.85 (0.64, 1.14) 0.75 (0.61, 0.92) 1.05 (0.65, 1.68)
C 1.07 (0.88, 1.3) 1.22 (1.03, 1.43) C 1.04 (0.82, 1.32) 0.91 (0.8, 1.04) 1.28 (0.82, 2.01)
D 1.03 (0.76, 1.41) 1.18 (0.88, 1.57) 0.96 (0.76, 1.22) D 0.88 (0.69, 1.13) 1.23 (0.84, 1.81)
E 1.17 (0.93, 1.48) 1.34 (1.08, 1.64) 1.1 (0.96, 1.25) 1.14 (0.89, 1.45) E 1.4 (0.89, 2.19)
F 0.84 (0.51, 1.37) 0.95 (0.59, 1.54) 0.78 (0.5, 1.22) 0.81 (0.55, 1.19) 0.71 (0.46, 1.12) 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()

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 7.061111e-02 0.120333333 0.253888889 0.34872222 0.201944444 0.004500000
## B 5.555556e-05 0.001333333 0.007222222 0.12822222 0.455000000 0.408166667
## C 4.072222e-02 0.423777778 0.355500000 0.15322222 0.025722222 0.001055556
## D 1.165000e-01 0.211666667 0.255555556 0.24111111 0.154833333 0.020333333
## E 7.158333e-01 0.196111111 0.069555556 0.01583333 0.002666667 0.000000000
## F 5.627778e-02 0.046777778 0.058277778 0.11288889 0.159833333 0.565944444

Extract model code (e.g. for Appendix)

cat(model$model$code)
## model {
##  # Likelihood for arm-based data
##  ## OMITTED
##  # Likelihood for contrast-based data (univariate for 2-arm trials)
##  for(i in studies.r2) {
##      for (k in 2:na[i]) {
##          mest[i, k] <- delta[i, k]
##      }
##      m[i, 2] ~ dnorm(mest[i, 2], prec[i, 2])
##      prec[i, 2] <- 1 / (e[i, 2] * e[i, 2])
##  
##      dev[i, 1] <- pow(m[i, 2] - mest[i, 2], 2) * prec[i, 2]
##  }
##  # Likelihood for contrast-based data (multivariate for multi-arm trials)
##  ## OMITTED
## 
##  # Fixed effect model
##  for (i in studies) {
##      delta[i, 1] <- 0
##      for (k in 2:na[i]) {
##          delta[i, k] <- d[t[i, 1], t[i, k]]
##      }
##  }
## 
##  # Relative effect matrix
##  d[1, 1] <- 0
##  d[1, 2] <- -d.B.A
##  d[1, 3] <- -d.B.A + -d.C.B
##  d[1, 4] <- -d.B.A + -d.C.B + d.C.D
##  d[1, 5] <- -d.B.A + -d.C.B + d.C.E
##  d[1, 6] <- -d.B.A + -d.C.B + d.C.D + d.D.F
##  for (i in 2:nt) {
##      for (j in 1:nt) {
##          d[i, j] <- d[1, j] - d[1, i]
##      }
##  }
## 
##  prior.prec <- pow(re.prior.sd, -2)
## 
##  # Study baseline priors
##  ## OMITTED
## 
##  # Effect parameter priors
##  d.B.A ~ dnorm(0, prior.prec)
##  d.C.B ~ dnorm(0, prior.prec)
##  d.C.D ~ dnorm(0, prior.prec)
##  d.C.E ~ dnorm(0, prior.prec)
##  d.D.F ~ dnorm(0, prior.prec)
##  
## }

Session info

BEE repository: /tmp/Rtmp9XTquv/Rbuild14e01251655f/gemtcPlus/vignettes

date()
## [1] "Fri Oct 11 04:20:22 2024"
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            reshape_0.8.9        
## [76] minqa_1.2.8           jsonlite_1.8.9        R6_2.5.1             
## [79] plyr_1.8.9            statnet.common_4.10.0 R2WinBUGS_2.1-22.1