tests/testthat/test-example_thrombolytics.R

In multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data

# Generated by vignette example_thrombolytics.Rmd: do not edit by hand
# Instead edit example_thrombolytics.Rmd and then run precompile.R

skip_on_cran()



params <-
list(run_tests = FALSE)

## ---- code=readLines("children/knitr_setup.R"), include=FALSE-------------------------------------


## ---- eval = FALSE--------------------------------------------------------------------------------
## library(multinma)
## options(mc.cores = parallel::detectCores())

## ----setup, echo = FALSE--------------------------------------------------------------------------
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), 
             "true" =, "warn" = 2, 
             parallel::detectCores())
options(mc.cores = nc)


## -------------------------------------------------------------------------------------------------
head(thrombolytics)


## -------------------------------------------------------------------------------------------------
thrombo_net <- set_agd_arm(thrombolytics, 
                           study = studyn,
                           trt = trtc,
                           r = r, 
                           n = n)
thrombo_net

## ---- include=FALSE, eval=params$run_tests--------------------------------------------------------
# Make trtf factor to order treatments in same way as Dias analysis - needed to
# recreate inconsistency analyses
trts <- dplyr::distinct(thrombolytics, trtn, trtc)
trts <- dplyr::arrange(trts, trtn)
thrombolytics$trtf <- factor(thrombolytics$trtn, levels = trts$trtn, labels = trts$trtc)
thrombo_net2 <- set_agd_arm(thrombolytics, 
                            study = studyn,
                            trt = trtf,
                            r = r, 
                            n = n)


## ---- eval=FALSE----------------------------------------------------------------------------------
## plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE)

## ----thrombo_net_plot, echo=FALSE-----------------------------------------------------------------
plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(legend.margin = ggplot2::margin(l = 4, unit = "lines"))


## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))


## -------------------------------------------------------------------------------------------------
thrombo_fit <- nma(thrombo_net, 
                   trt_effects = "fixed",
                   prior_intercept = normal(scale = 100),
                   prior_trt = normal(scale = 100))


## -------------------------------------------------------------------------------------------------
thrombo_fit


## ---- eval=FALSE----------------------------------------------------------------------------------
## # Not run
## print(thrombo_fit, pars = c("d", "mu"))


## ----thrombo_pp_plot------------------------------------------------------------------------------
plot_prior_posterior(thrombo_fit, prior = "trt")


## -------------------------------------------------------------------------------------------------
(dic_consistency <- dic(thrombo_fit))


## ----thrombo_resdev_plot, fig.width=12------------------------------------------------------------
plot(dic_consistency)


## -------------------------------------------------------------------------------------------------
thrombo_fit_ume <- nma(thrombo_net, 
                       consistency = "ume",
                       trt_effects = "fixed",
                       prior_intercept = normal(scale = 100),
                       prior_trt = normal(scale = 100))
thrombo_fit_ume


## -------------------------------------------------------------------------------------------------
dic_consistency
(dic_ume <- dic(thrombo_fit_ume))


## ----thrombo_devdev_plot--------------------------------------------------------------------------
plot(dic_consistency, dic_ume, show_uncertainty = FALSE)


## ---- eval=!params$run_tests----------------------------------------------------------------------
## thrombo_nodesplit <- nma(thrombo_net,
##                          consistency = "nodesplit",
##                          trt_effects = "fixed",
##                          prior_intercept = normal(scale = 100),
##                          prior_trt = normal(scale = 100))

## ---- include=FALSE, eval=params$run_tests--------------------------------------------------------
# Run node-splits with treatments ordered as per Dias
thrombo_nodesplit <- nma(thrombo_net2, 
                         consistency = "nodesplit",
                         trt_effects = "fixed",
                         prior_intercept = normal(scale = 100),
                         prior_trt = normal(scale = 100))


## -------------------------------------------------------------------------------------------------
summary(thrombo_nodesplit)


## ----thrombo_nodesplit, fig.width = 7-------------------------------------------------------------
plot(thrombo_nodesplit)


## ----thrombo_nodesplit_omega, fig.height = 6------------------------------------------------------
plot(thrombo_nodesplit, pars = "omega", stat = "halfeye", ref_line = 0) +
  ggplot2::aes(y = comparison) +
  ggplot2::facet_null()


## ----thrombo_releff-------------------------------------------------------------------------------
(thrombo_releff <- relative_effects(thrombo_fit, all_contrasts = TRUE))
plot(thrombo_releff, ref_line = 0)


## ----thrombo_ranks--------------------------------------------------------------------------------
(thrombo_ranks <- posterior_ranks(thrombo_fit))
plot(thrombo_ranks)

## ----thrombo_rankprobs----------------------------------------------------------------------------
(thrombo_rankprobs <- posterior_rank_probs(thrombo_fit))
plot(thrombo_rankprobs)

## ----thrombo_cumrankprobs-------------------------------------------------------------------------
(thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE))
plot(thrombo_cumrankprobs)


## ----thrombo_tests, include=FALSE, eval=params$run_tests------------------------------------------
#--- Test against TSD 4 results ---
library(testthat)
library(dplyr)
library(tidyr)

test_that("Reference trt is SK", {
  expect_equivalent(levels(thrombo_net$treatments)[1], "SK")
})

tol <- 0.05
tol_dic <- 0.1

# Relative effects
tsd_releff <- tribble(
~trt_b    , ~trt       , ~mean , ~sd  , ~lower, ~upper,
"SK"      , "t-PA"     , 0.002 , 0.030, -0.06 , 0.06  ,
"SK"      , "Acc t-PA" , -0.177, 0.043, -0.26 , -0.09 ,
"SK"      , "SK + t-PA", -0.049, 0.046, -0.14 , 0.04  ,
"SK"      , "r-PA"     , -0.124, 0.060, -0.24 , -0.01 ,
"SK"      , "PTCA"     , -0.476, 0.101, -0.67 , -0.28 ,
"SK"      , "UK"       , -0.203, 0.221, -0.64 , 0.23  ,
"SK"      , "ASPAC"    , 0.016 , 0.037, -0.06 , 0.09  ,
"t-PA"    , "PTCA"     , -0.478, 0.104, -0.68 , -0.27 ,
"t-PA"    , "UK"       , -0.206, 0.221, -0.64 , 0.23  ,
"t-PA"    , "ASPAC"    , 0.013 , 0.037, -0.06 , 0.09  ,
"Acc t-PA", "r-PA"     , 0.054 , 0.055, -0.05 , 0.16  ,
"Acc t-PA", "TNK"      , 0.005 , 0.064, -0.12 , 0.13  ,
"Acc t-PA", "PTCA"     , -0.298, 0.098, -0.49 , -0.11 ,
"Acc t-PA", "UK"       , -0.026, 0.221, -0.46 , 0.41  ,
"Acc t-PA", "ASPAC"    , 0.193 , 0.056, 0.08  , 0.30  ) %>% 
  mutate(.trt_b = ordered(trt_b, levels = levels(thrombo_net$treatments)),
         .trt = ordered(trt, levels = levels(thrombo_net$treatments)),
         rev = if_else(.trt_b > .trt, -1, 1),
         .l = lower, .u = upper,
         lower = if_else(.trt_b > .trt, .u, .l),
         upper = if_else(.trt_b > .trt, .l, .u),
         trt_b = if_else(.trt_b > .trt, .trt, .trt_b),
         trt = if_else(.trt_b > .trt, .trt_b, .trt),
         lab = paste0("d[", trt, " vs. ", trt_b, "]")) %>% 
  arrange(trt_b, trt) %>% 
  mutate_at(vars(mean, lower, upper), ~.*rev)

thrombo_releff_summary <- as.data.frame(thrombo_releff) %>% 
  filter(parameter %in% tsd_releff$lab)

test_that("FE relative effects", {
  expect_equivalent(thrombo_releff_summary$mean, tsd_releff$mean, tolerance = tol)
  expect_equivalent(thrombo_releff_summary$sd, tsd_releff$sd, tolerance = tol)
  expect_equivalent(thrombo_releff_summary$`2.5%`, tsd_releff$lower, tolerance = tol)
  expect_equivalent(thrombo_releff_summary$`97.5%`, tsd_releff$upper, tolerance = tol)
})

test_that("SUCRAs", {
  thrombo_ranks <- posterior_ranks(thrombo_fit, sucra = TRUE)
  thrombo_rankprobs <- posterior_rank_probs(thrombo_fit, sucra = TRUE)
  thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE, sucra = TRUE)
  
  expect_equal(thrombo_ranks$summary$sucra, thrombo_rankprobs$summary$sucra)
  expect_equal(thrombo_ranks$summary$sucra, thrombo_cumrankprobs$summary$sucra)
})

# DIC
test_that("DIC", {
  expect_equivalent(dic_consistency$resdev, 105.9, tolerance = tol_dic)
  expect_equivalent(dic_consistency$pd, 58, tolerance = tol_dic)
  expect_equivalent(dic_consistency$dic, 163.9, tolerance = tol_dic)
})

# Relative effects (UME)

# FE UME model, so no differences by reference treatment, no multi-arm correction
tsd_ume <- tribble(
~trt_b    , ~trt       , ~mean , ~sd  , ~lower, ~upper,
"SK"      , "t-PA"     , -0.004, 0.030, -0.06 , 0.06  ,
"SK"      , "Acc t-PA" , -0.158, 0.049, -0.25 , -0.06 ,
"SK"      , "SK + t-PA", -0.044, 0.047, -0.14 , 0.05  ,
"SK"      , "r-PA"     , -0.060, 0.089, -0.23 , 0.11  ,
"SK"      , "PTCA"     , -0.665, 0.185, -1.03 , -0.31 ,
"SK"      , "UK"       , -0.369, 0.518, -1.41 , 0.63  ,
"SK"      , "ASPAC"    , 0.005 , 0.037, -0.07 , 0.08  ,
"t-PA"    , "PTCA"     , -0.544, 0.417, -1.38 , 0.25  ,
"t-PA"    , "UK"       , -0.294, 0.347, -0.99 , 0.37  ,
"t-PA"    , "ASPAC"    , -0.290, 0.361, -1.01 , 0.41  ,
"Acc t-PA", "r-PA"     , 0.019 , 0.066, -0.11 , 0.15  ,
"Acc t-PA", "TNK"      , 0.006 , 0.064, -0.12 , 0.13  ,
"Acc t-PA", "PTCA"     , -0.216, 0.119, -0.45 , 0.02  ,
"Acc t-PA", "UK"       , 0.146 , 0.358, -0.54 , 0.86  ,
"Acc t-PA", "ASPAC"    , 1.405 , 0.417, 0.63  , 2.27  ) %>% 
  mutate(.trt_b = ordered(trt_b, levels = levels(thrombo_net$treatments)),
         .trt = ordered(trt, levels = levels(thrombo_net$treatments)),
         rev = if_else(.trt_b > .trt, -1, 1),
         .l = lower, .u = upper,
         lower = if_else(.trt_b > .trt, .u, .l),
         upper = if_else(.trt_b > .trt, .l, .u),
         trt_b = if_else(.trt_b > .trt, .trt, .trt_b),
         trt = if_else(.trt_b > .trt, .trt_b, .trt),
         lab = paste0("d[", trt, " vs. ", trt_b, "]")) %>%
  arrange(trt_b, trt) %>%
  mutate_at(vars(mean, lower, upper), ~.*rev)

thrombo_ume_releff <- summary(thrombo_fit_ume, pars = "d")

test_that("UME relative effects", {
  expect_equivalent(thrombo_ume_releff$summary$mean, tsd_ume$mean, tolerance = tol)
  expect_equivalent(thrombo_ume_releff$summary$sd, tsd_ume$sd, tolerance = tol)
  expect_equivalent(thrombo_ume_releff$summary$`2.5%`, tsd_ume$lower, tolerance = tol)
  expect_equivalent(thrombo_ume_releff$summary$`97.5%`, tsd_ume$upper, tolerance = tol)
})

# DIC (UME)
test_that("UME DIC", {
  expect_equivalent(dic_ume$resdev, 99.7, tolerance = tol_dic)
  expect_equivalent(dic_ume$pd, 65, tolerance = tol_dic)
  expect_equivalent(dic_ume$dic, 164.7, tolerance = tol_dic)
})


# Node-splitting
trt_code <- c("SK", "t-PA", "Acc t-PA", "SK + t-PA", "r-PA",
              "TNK", "PTCA", "UK", "ASPAC")

dias2010_nodesplit_est <- tribble(
  ~trt1, ~trt2, ~net_mean, ~net_sd, ~dir_mean, ~dir_sd, ~ind_mean, ~ind_sd, ~omega_mean, ~omega_sd, ~p_value,
 1, 2,  0.002,  0.030,  0.000,  0.030,  0.189,  0.235, -0.190,  0.236,  0.42,
 1, 3, -0.177,  0.043, -0.158,  0.048, -0.247,  0.092,  0.088,  0.104,  0.39,
 1, 5, -0.124,  0.060, -0.060,  0.089, -0.175,  0.081,  0.115,  0.121,  0.34,
 1, 7, -0.475,  0.101, -0.666,  0.185, -0.393,  0.120, -0.272,  0.222,  0.22,
 1, 8, -0.203,  0.219, -0.369,  0.518, -0.168,  0.244, -0.207,  0.575,  0.73,
 1, 9,  0.016,  0.037,  0.009,  0.037,  0.424,  0.252, -0.413,  0.253,  0.10,
 2, 7, -0.477,  0.104, -0.545,  0.417, -0.475,  0.108, -0.073,  0.432,  0.88,
 2, 8, -0.205,  0.220, -0.295,  0.347, -0.144,  0.290, -0.155,  0.452,  0.74,
 2, 9,  0.014,  0.037,  0.006,  0.037,  0.471,  0.241, -0.468,  0.241,  0.05,
 3, 4,  0.128,  0.054,  0.126,  0.054,  0.630,  0.697, -0.506,  0.696,  0.47,
 3, 5,  0.053,  0.056,  0.019,  0.066,  0.135,  0.101, -0.116,  0.120,  0.34,
 3, 7, -0.298,  0.097, -0.216,  0.118, -0.477,  0.174,  0.260,  0.211,  0.21,
 3, 8, -0.026,  0.220,  0.143,  0.356, -0.136,  0.288,  0.277,  0.461,  0.55,
 3, 9,  0.193,  0.056,  1.409,  0.415,  0.165,  0.057,  1.239,  0.420,  0.001
) %>% 
  mutate(trt1 = forcats::fct_relevel(factor(trt1, levels = 1:9, labels = trt_code), 
                                     !! levels(thrombo_net2$treatments)),
         trt2 = forcats::fct_relevel(factor(trt2, levels = 1:9, labels = trt_code), 
                                     !! levels(thrombo_net2$treatments)))

for (i in 1:nrow(dias2010_nodesplit_est)) {
  if (as.numeric(dias2010_nodesplit_est$trt1[i]) > as.numeric(dias2010_nodesplit_est$trt2[i])) {
    dias2010_nodesplit_est[i, c("trt1", "trt2")] <- dias2010_nodesplit_est[i, c("trt2", "trt1")]
    dias2010_nodesplit_est[i, c("net_mean", "dir_mean", "ind_mean")] <- -dias2010_nodesplit_est[i, c("net_mean", "dir_mean", "ind_mean")]
  }
}

dias2010_nodesplit_est <- arrange(dias2010_nodesplit_est, trt1, trt2)

thrombo_nodesplit_est <- as_tibble(summary(thrombo_nodesplit), nest = FALSE) %>% 
  mutate(parameter = stringr::str_replace(parameter, "(^d_)?(.+)\\[.+\\]$", "\\2")) %>% 
  pivot_wider(names_from = "parameter", 
              values_from = c("mean", "sd"),
              names_glue = "{parameter}_{.value}",
              id_cols = c("trt1", "trt2", "p_value")) %>% 
  mutate(across(where(is.numeric), unname)) %>% 
  select(!!! colnames(dias2010_nodesplit_est))

test_that("Node-splitting estimates", {
  # expect_equal(thrombo_nodesplit_est, dias2010_nodesplit_est, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$net_mean, dias2010_nodesplit_est$net_mean, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$net_sd, dias2010_nodesplit_est$net_sd, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$ind_mean, dias2010_nodesplit_est$ind_mean, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$ind_sd, dias2010_nodesplit_est$ind_sd, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$dir_mean, dias2010_nodesplit_est$dir_mean, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$dir_sd, dias2010_nodesplit_est$dir_sd, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$omega_mean, dias2010_nodesplit_est$omega_mean, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$omega_sd, dias2010_nodesplit_est$omega_sd, tolerance = tol)
  expect_equal(thrombo_nodesplit_est$p_value, dias2010_nodesplit_est$p_value, tolerance = tol)
})

dias2010_nodesplit_dic <- tribble(
  ~trt1, ~trt2, ~resdev,   ~pd,   ~dic,
      1,     2,   106.4,  58.7,  165.1,
      1,     3,   106.2,  58.7,  165.0,
      1,     5,   106.1,  58.7,  164.8,
      1,     7,   105.5,  58.7,  164.2,
      1,     8,   106.9,  58.7,  165.6,
      1,     9,   104.3,  58.7,  163.0,
      2,     7,   106.9,  58.7,  165.6,
      2,     8,   106.8,  58.7,  165.5,
      2,     9,   103.3,  58.7,  162.0,
      3,     4,   106.5,  58.7,  165.2,
      3,     5,   106.1,  58.7,  164.8,
      3,     7,   105.5,  58.7,  164.2,
      3,     8,   106.6,  58.7,  165.3,
      3,     9,    96.9,  58.7,  155.6) %>% 
  mutate(trt1 = forcats::fct_relevel(factor(trt1, levels = 1:9, labels = trt_code), 
                                     !! levels(thrombo_net$treatments)),
         trt2 = forcats::fct_relevel(factor(trt2, levels = 1:9, labels = trt_code), 
                                     !! levels(thrombo_net$treatments)))

for (i in 1:nrow(dias2010_nodesplit_dic)) {
  if (as.numeric(dias2010_nodesplit_dic$trt1[i]) > as.numeric(dias2010_nodesplit_dic$trt2[i])) {
    dias2010_nodesplit_dic[i, c("trt1", "trt2")] <- dias2010_nodesplit_dic[i, c("trt2", "trt1")]
  }
}

dias2010_nodesplit_dic <- arrange(dias2010_nodesplit_dic, trt1, trt2)

thrombo_nodesplit_dic <- as_tibble(summary(thrombo_nodesplit), nest = FALSE) %>% 
  distinct(trt1, trt2, resdev, pd, dic)

test_that("Node-splitting DIC", {
  expect_equal(thrombo_nodesplit_dic$dic, dias2010_nodesplit_dic$dic, tolerance = tol_dic)
  expect_equal(thrombo_nodesplit_dic$resdev, dias2010_nodesplit_dic$resdev, tolerance = tol_dic)
  expect_equal(thrombo_nodesplit_dic$pd, dias2010_nodesplit_dic$pd, tolerance = tol_dic)
})