tests/testthat/test-example_smoking.R

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

# Generated by vignette example_smoking.Rmd: do not edit by hand
# Instead edit example_smoking.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(smoking)


## -------------------------------------------------------------------------------------------------
smknet <- set_agd_arm(smoking, 
                      study = studyn,
                      trt = trtc,
                      r = r, 
                      n = n,
                      trt_ref = "No intervention")
smknet


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

## ----smoking_network_plot, echo=FALSE, fig.width=8, fig.height=6, out.width="100%"----------------
plot(smknet, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(plot.margin = ggplot2::unit(c(1, 1, 1, 6), "lines"))


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

## -------------------------------------------------------------------------------------------------
summary(half_normal(scale = 5))


## -------------------------------------------------------------------------------------------------
smkfit <- nma(smknet, 
              trt_effects = "random",
              prior_intercept = normal(scale = 100),
              prior_trt = normal(scale = 100),
              prior_het = normal(scale = 5))


## -------------------------------------------------------------------------------------------------
smkfit


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


## ----smoking_pp_plot, fig.width=8, fig.height=6, out.width="100%"---------------------------------
plot_prior_posterior(smkfit)


## ----smoking_pp_plot_tau--------------------------------------------------------------------------
plot_prior_posterior(smkfit, prior = "het")


## -------------------------------------------------------------------------------------------------
(dic_consistency <- dic(smkfit))


## ----smoking_resdev_plot, fig.width=8-------------------------------------------------------------
plot(dic_consistency)


## -------------------------------------------------------------------------------------------------
smoking[smoking$r == 0, ]


## -------------------------------------------------------------------------------------------------
smkfit_ume <- nma(smknet, 
                  consistency = "ume",
                  trt_effects = "random",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 100),
                  prior_het = normal(scale = 5))
smkfit_ume


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


## ----smoking_devdev_plot--------------------------------------------------------------------------
plot(dic_consistency, dic_ume, point_alpha = 0.5, interval_alpha = 0.2)


## -------------------------------------------------------------------------------------------------
smk_nodesplit <- nma(smknet, 
                     consistency = "nodesplit",
                     trt_effects = "random",
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_het = normal(scale = 5))


## -------------------------------------------------------------------------------------------------
summary(smk_nodesplit)


## ----smk_nodesplit, fig.width = 7-----------------------------------------------------------------
plot(smk_nodesplit) +
  ggplot2::theme(legend.position = "bottom", legend.direction = "horizontal")


## ----smoking_releff, fig.height=4.5---------------------------------------------------------------
(smk_releff <- relative_effects(smkfit, all_contrasts = TRUE))
plot(smk_releff, ref_line = 0)


## ----smoking_ranks, fig.height=3------------------------------------------------------------------
(smk_ranks <- posterior_ranks(smkfit, lower_better = FALSE))
plot(smk_ranks)

## ----smoking_rankprobs----------------------------------------------------------------------------
(smk_rankprobs <- posterior_rank_probs(smkfit, lower_better = FALSE))
plot(smk_rankprobs)

## ----smoking_cumrankprobs-------------------------------------------------------------------------
(smk_cumrankprobs <- posterior_rank_probs(smkfit, lower_better = FALSE, cumulative = TRUE))
plot(smk_cumrankprobs)


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

tol <- 0.05
tol_dic <- 0.1

# A = No intervention
# B = Self-help
# C = Individual counselling
# D = Group counselling

trt_codes <- c(A = "No intervention",
               B = "Self-help",
               C = "Individual counselling",
               D = "Group counselling")

# Relative effects
tsd_re <- tribble(
  ~contrast, ~est, ~sd, ~lower, ~upper,
  "dAB", 0.49, 0.40, -0.29, 1.31,
  "dAC", 0.84, 0.24,  0.39, 1.34,
  "dAD", 1.10, 0.44,  0.26, 2.00,
  "dBC", 0.35, 0.41, -0.46, 1.18,
  "dBD", 0.61, 0.49, -0.34, 1.59,
  "dCD", 0.26, 0.41, -0.55, 1.09) %>% 
  mutate(trt_b = recode(substr(contrast, 2, 2),
                        !!! trt_codes),
         trt = recode(substr(contrast, 3, 3),
                        !!! trt_codes),
         trt_b = ordered(trt_b, levels = levels(smknet$treatments)),
         trt = ordered(trt, levels = levels(smknet$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)) %>% 
  arrange(trt_b, trt) %>% 
  mutate_at(vars(est, lower, upper), ~.*rev)

test_that("RE relative effects", {
  expect_equivalent(smk_releff$summary$mean, tsd_re$est, tolerance = tol)
  expect_equivalent(smk_releff$summary$sd, tsd_re$sd, tolerance = tol)
  expect_equivalent(smk_releff$summary$`2.5%`, tsd_re$lower, tolerance = tol)
  expect_equivalent(smk_releff$summary$`97.5%`, tsd_re$upper, tolerance = tol)
})

# Heterogeneity SD
smk_tau <- summary(smkfit, pars = "tau")

test_that("RE heterogeneity SD", {
  expect_equivalent(smk_tau$summary$`50%`, 0.82, tolerance = tol)
  expect_equivalent(smk_tau$summary$sd, 0.19, tolerance = tol)
  expect_equivalent(smk_tau$summary$`2.5%`, 0.55, tolerance = tol)
  expect_equivalent(smk_tau$summary$`97.5%`, 1.27, tolerance = tol)
})

# DIC
test_that("DIC", {
  expect_equivalent(dic_consistency$resdev, 54.0, tolerance = tol_dic)
  expect_equivalent(dic_consistency$pd, 45.0, tolerance = tol_dic)
  expect_equivalent(dic_consistency$dic, 99.0, tolerance = tol_dic)
})

# Relative effects (UME)

# Treatment ordering is different in TSD 4 - use trtn instead which reflects this
smknet2 <- set_agd_arm(smoking, studyn, trtn, r = r, n = n, trt_ref = 1)
smkfit_ume2 <- nma(smknet2, 
                   consistency = "ume",
                   trt_effects = "random",
                   prior_intercept = normal(scale = 100),
                   prior_trt = normal(scale = 100),
                   prior_het = normal(scale = 5),
                   iter = 4000)

# Results of modified TSD 4 model including multi-arm correction
tsd_ume <- tribble(
  ~contrast, ~est, ~sd, ~lower, ~upper,
  # "dAB",  0.34, 0.58, -0.81, 1.50,
  # "dAC",  0.86, 0.27,  0.34, 1.43,
  # "dAD",  1.43, 0.88, -0.21, 3.29,
  # "dBC", -0.05, 0.74, -1.53, 1.42,
  # "dBD",  0.65, 0.73, -0.80, 2.12,
  # "dCD",  0.20, 0.78, -1.37, 1.73)  %>% 
  "dAB",  0.34, 0.59, -0.82, 1.51,
  "dAC",  0.90, 0.28,  0.37, 1.49,
  "dAD",  1.12, 0.80, -0.36, 2.80,
  "dBC",  0.04, 0.73, -1.39, 1.50,
  "dBD",  0.61, 0.72, -0.80, 2.05,
  "dCD",  0.20, 0.79, -1.40, 1.75)
  # mutate(trt_b = recode(substr(contrast, 2, 2),
  #                       !!! trt_codes),
  #        trt = recode(substr(contrast, 3, 3),
  #                       !!! trt_codes),
  #        trt_b = ordered(trt_b, levels = levels(smknet$treatments)),
  #        trt = ordered(trt, levels = levels(smknet$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)) %>% 
  # arrange(trt_b, trt) %>% 
  # mutate_at(vars(est:upper), ~.*rev)

smk_ume_releff <- summary(smkfit_ume2, pars = "d")

test_that("UME relative effects", {
  # skip("Different model parameterisation")
  expect_equivalent(smk_ume_releff$summary$mean, tsd_ume$est, tolerance = tol)
  expect_equivalent(smk_ume_releff$summary$sd, tsd_ume$sd, tolerance = tol)
  expect_equivalent(smk_ume_releff$summary$`2.5%`, tsd_ume$lower, tolerance = tol)
  expect_equivalent(smk_ume_releff$summary$`97.5%`, tsd_ume$upper, tolerance = tol)
})

# Heterogeneity SD (UME)
smk_ume_tau <- summary(smkfit_ume2, pars = "tau")

test_that("UME heterogeneity SD", {
  # skip("Different model parameterisation")
  # expect_equivalent(smk_ume_tau$summary$`50%`, 0.89, tolerance = tol)
  # expect_equivalent(smk_ume_tau$summary$sd, 0.22, tolerance = tol)
  # expect_equivalent(smk_ume_tau$summary$`2.5%`, 0.58, tolerance = tol)
  # expect_equivalent(smk_ume_tau$summary$`97.5%`, 1.45, tolerance = tol)
  expect_equivalent(smk_ume_tau$summary$`50%`, 0.91, tolerance = tol)
  expect_equivalent(smk_ume_tau$summary$sd, 0.21, tolerance = tol)
  expect_equivalent(smk_ume_tau$summary$`2.5%`, 0.59, tolerance = tol)
  expect_equivalent(smk_ume_tau$summary$`97.5%`, 1.48, tolerance = tol)
})

# DIC (UME)
test_that("UME DIC", {
  expect_equivalent(dic_ume$resdev, 53.4, tolerance = tol_dic)
  expect_equivalent(dic_ume$pd, 46.1, tolerance = tol_dic)
  expect_equivalent(dic_ume$dic, 99.5, tolerance = tol_dic)
})

dic_ume2 <- dic(smkfit_ume2)
test_that("UME DIC", {
  expect_equivalent(dic_ume2$resdev, 53.4, tolerance = tol_dic)
  expect_equivalent(dic_ume2$pd, 46.1, tolerance = tol_dic)
  expect_equivalent(dic_ume2$dic, 99.5, tolerance = tol_dic)
})

# Check that multinomial ordered model produces same results
smknet3 <- set_agd_arm(smoking,
                       studyn,
                       trtc,
                       r = multi(nonevent = n, event = r, inclusive = TRUE),
                       trt_ref = "No intervention")
smkfit_ord <- nma(smknet3,
              trt_effects = "random",
              link = "logit",
              prior_intercept = normal(scale = 100),
              prior_trt = normal(scale = 100),
              prior_het = normal(scale = 5),
              prior_aux = flat())


smk_ord_releff <- relative_effects(smkfit_ord, all_contrasts = TRUE)
test_that("Equivalent ordered multinomial RE relative effects", {
  expect_equivalent(smk_ord_releff$summary$mean, tsd_re$est, tolerance = tol)
  expect_equivalent(smk_ord_releff$summary$sd, tsd_re$sd, tolerance = tol)
  expect_equivalent(smk_ord_releff$summary$`2.5%`, tsd_re$lower, tolerance = tol)
  expect_equivalent(smk_ord_releff$summary$`97.5%`, tsd_re$upper, tolerance = tol)
})

smk_ord_tau <- summary(smkfit_ord, pars = "tau")
test_that("Equivalent ordered multinomial RE heterogeneity SD", {
  expect_equivalent(smk_ord_tau$summary$`50%`, 0.82, tolerance = tol)
  expect_equivalent(smk_ord_tau$summary$sd, 0.19, tolerance = tol)
  expect_equivalent(smk_ord_tau$summary$`2.5%`, 0.55, tolerance = tol)
  expect_equivalent(smk_ord_tau$summary$`97.5%`, 1.27, tolerance = tol)
})

# DIC
dic_ord <- dic(smkfit_ord)
test_that("Equivalent ordered multinomial DIC", {
  expect_equivalent(dic_ord$resdev, 54.0, tolerance = tol_dic)
  expect_equivalent(dic_ord$pd, 45.0, tolerance = tol_dic)
  expect_equivalent(dic_ord$dic, 99.0, tolerance = tol_dic)
})

test_that("Robust to custom options(contrasts) settings", {
  skip_on_cran()
  
  withr::with_options(list(contrasts = c(ordered = "contr.SAS",
                                                       unordered = "contr.SAS")), {
    smkfit_SAS <- nma(smknet,
                      trt_effects = "random",
                      prior_intercept = normal(scale = 100),
                      prior_trt = normal(scale = 100),
                      prior_het = normal(scale = 5))
    
    smkfit_SAS_summary <- as_tibble(summary(smkfit_SAS))[, c("parameter", "mean", "sd")]
    smkfit_SAS_releff <- as_tibble(relative_effects(smkfit_SAS))[, c("parameter", "mean", "sd")]
    smkfit_SAS_pred <- as_tibble(predict(smkfit_SAS))[, c("parameter", "mean", "sd")]
    
    smkfit_ume_SAS <- nma(smknet,
                      trt_effects = "random",
                      consistency = "ume",
                      prior_intercept = normal(scale = 100),
                      prior_trt = normal(scale = 100),
                      prior_het = normal(scale = 5))
    
    smkfit_ume_SAS_summary <- as_tibble(summary(smkfit_ume_SAS))[, c("parameter", "mean", "sd")]
  })

  expect_equal(smkfit_SAS_summary,
               as_tibble(summary(smkfit))[, c("parameter", "mean", "sd")],
               tolerance = tol)
  expect_equal(smkfit_SAS_releff,
               as_tibble(relative_effects(smkfit))[, c("parameter", "mean", "sd")],
               tolerance = tol)
  expect_equal(smkfit_SAS_pred,
               as_tibble(predict(smkfit))[, c("parameter", "mean", "sd")],
               tolerance = tol)
  
  expect_equal(smkfit_ume_SAS_summary,
               as_tibble(summary(smkfit_ume))[, c("parameter", "mean", "sd")],
               tolerance = tol)
})


# Force clean up
rm(list = ls())
gc()