tests/testthat/test-example_bcg_vaccine.R

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

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

skip_on_cran()



params <-
list(run_tests = FALSE)

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

## ---- include=FALSE-------------------------------------------------------------------------------
set.seed(18284729)


## ---- 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(bcg_vaccine)


## -------------------------------------------------------------------------------------------------
bcg_net <- set_agd_arm(bcg_vaccine, 
                       study = studyn,
                       trt = trtc,
                       r = r, 
                       n = n,
                       trt_ref = "Unvaccinated")
bcg_net


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


## ---- eval = FALSE--------------------------------------------------------------------------------
## bcg_fit_unadj <- nma(bcg_net,
##                      trt_effects = "random",
##                      prior_intercept = normal(scale = 100),
##                      prior_trt = normal(scale = 100),
##                      prior_het = half_normal(scale = 5))

## ---- echo = FALSE--------------------------------------------------------------------------------
bcg_fit_unadj <- nma(bcg_net, 
                     seed = 14308133,
                     trt_effects = "random",
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_het = half_normal(scale = 5))


## -------------------------------------------------------------------------------------------------
bcg_fit_unadj


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


## ----bcg_unadj_pp_plot----------------------------------------------------------------------------
plot_prior_posterior(bcg_fit_unadj, prior = c("trt", "het"))


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


## ---- eval = FALSE--------------------------------------------------------------------------------
## bcg_fit_lat <- nma(bcg_net,
##                    trt_effects = "random",
##                    regression = ~.trt:latitude,
##                    prior_intercept = normal(scale = 100),
##                    prior_trt = normal(scale = 100),
##                    prior_reg = normal(scale = 100),
##                    prior_het = half_normal(scale = 5),
##                    adapt_delta = 0.99)

## ---- echo = FALSE--------------------------------------------------------------------------------
bcg_fit_lat <- nowarn_on_ci(
                 nma(bcg_net, 
                     seed = 1932599147,
                     trt_effects = "random",
                     regression = ~.trt:latitude,
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_reg = normal(scale = 100),
                     prior_het = half_normal(scale = 5),
                     adapt_delta = 0.99)
                 )


## -------------------------------------------------------------------------------------------------
bcg_fit_lat


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


## ----bcg_lat_pp_plot------------------------------------------------------------------------------
plot_prior_posterior(bcg_fit_lat, prior = c("trt", "reg", "het"))


## -------------------------------------------------------------------------------------------------
(bcg_dic_unadj <- dic(bcg_fit_unadj))

## -------------------------------------------------------------------------------------------------
(bcg_dic_lat <- dic(bcg_fit_lat))


## -------------------------------------------------------------------------------------------------
summary(bcg_fit_unadj, pars = "tau")
summary(bcg_fit_lat, pars = "tau")


## ----bcg_vaccine_beta_lat, fig.height = 4---------------------------------------------------------
summary(bcg_fit_lat, pars = "beta")

plot(bcg_fit_lat, 
     pars = "beta", 
     ref_line = 0,
     stat = "halfeye")


## -------------------------------------------------------------------------------------------------
bcg_releff_lat <- relative_effects(bcg_fit_lat,
                                   newdata = tibble::tibble(latitude = seq(10, 50, by = 10),
                                                            label = paste0(latitude, "\u00B0 latitude")),
                                   study = label)

bcg_releff_lat


## ----bcg_vaccine_releff_lat, fig.height = 5-------------------------------------------------------
plot(bcg_releff_lat, 
     ref_line = 0)


## ----bcg_vaccine_reg_plot-------------------------------------------------------------------------
library(dplyr)
library(ggplot2)

# Get data for regression line
lat_range <- range(bcg_vaccine$latitude)
lat_dat <- tibble(latitude = seq(lat_range[1], lat_range[2], by = 1))

bcg_lat_reg <- relative_effects(bcg_fit_lat, 
                                newdata = lat_dat) %>% 
  as_tibble() %>% 
  bind_cols(lat_dat)

# Get study log odds ratios
bcg_lor <- bcg_vaccine %>% 
  group_by(studyn) %>% 
  mutate(lor = log(r / (n - r)) - log(first(r) / (first(n) - first(r))),
         sample_size = sum(n)) %>% 
  slice(-1)

# Plot
ggplot(aes(x = latitude), data = bcg_lor) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), data = bcg_lat_reg,
              fill = "darkred", alpha = 0.3) +
  geom_line(aes(y = mean), data = bcg_lat_reg,
            colour = "darkred") +
  geom_point(aes(y = lor, size = sample_size), alpha = 0.6) +
  coord_cartesian(xlim = c(0, 60)) +
  xlab("Degrees Latitude") + ylab("log Odds Ratio") +
  scale_size("Sample Size") +
  theme_multinma()


## ----bcg_vaccine_predictive_unadj-----------------------------------------------------------------
(bcg_predeff_unadj <- relative_effects(bcg_fit_unadj, predictive_distribution = TRUE))


## -------------------------------------------------------------------------------------------------
mean(as.matrix(bcg_predeff_unadj) > 0)


## -------------------------------------------------------------------------------------------------
bcg_predeff_lat <- relative_effects(bcg_fit_lat,
                                   newdata = tibble::tibble(latitude = seq(0, 50, by = 10),
                                                            label = paste0(latitude, "\u00B0 latitude")),
                                   study = label,
                                   predictive_distribution = TRUE)

bcg_predeff_lat


## -------------------------------------------------------------------------------------------------
colMeans(as.matrix(bcg_predeff_lat) > 0)


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

tol <- 0.05
tol_dic <- 0.1

# Relative effects
bcg_unadj_releff <- as.data.frame(summary(bcg_fit_unadj, pars = "d"))

test_that("Unadjusted relative effects", {
  expect_equivalent(bcg_unadj_releff$mean, -0.762, tolerance = tol)
  expect_equivalent(bcg_unadj_releff$sd, 0.22, tolerance = tol)
  expect_equivalent(bcg_unadj_releff$`2.5%`, -1.21, tolerance = tol)
  expect_equivalent(bcg_unadj_releff$`97.5%`, -0.34, tolerance = tol)
})

test_that("Unadjusted predictive distribution", {
  bcg_unadj_releff_pred <- as.data.frame(relative_effects(bcg_fit_unadj, predictive_distribution = TRUE))
  expect_equivalent(bcg_unadj_releff_pred$mean, -0.762, tolerance = tol)
  expect_equivalent(bcg_unadj_releff_pred$`2.5%`, -2.27, tolerance = tol)
  expect_equivalent(bcg_unadj_releff_pred$`97.5%`, 0.72, tolerance = tol)
})

bcg_lat_releff <- as.data.frame(summary(bcg_fit_lat, pars = "d"))

test_that("Regression relative effects", {
  expect_equivalent(bcg_lat_releff$mean, -0.763, tolerance = tol)
  expect_equivalent(bcg_lat_releff$sd, 0.126, tolerance = tol)
  expect_equivalent(bcg_lat_releff$`2.5%`, -1.04, tolerance = tol)
  expect_equivalent(bcg_lat_releff$`97.5%`, -0.52, tolerance = tol)
})


test_that("Regression predictive distribution", {
  bcg_lat_releff_pred <- relative_effects(bcg_fit_lat,
                                          newdata = data.frame(latitude = c(0, 13, 50)),
                                          predictive_distribution = TRUE)
  expect_equivalent(colMeans(as.matrix(bcg_lat_releff_pred) > 0), c(0.8, 0.35, 0.006), tolerance = tol)
})

# Regression coefficients
bcg_lat_beta <- as.data.frame(summary(bcg_fit_lat, pars = "beta"))

test_that("Regression beta", {
  expect_equivalent(bcg_lat_beta$mean, -0.032, tolerance = tol)
  expect_equivalent(bcg_lat_beta$sd, 0.009, tolerance = tol)
  expect_equivalent(bcg_lat_beta$`2.5%`, -0.05, tolerance = tol)
  expect_equivalent(bcg_lat_beta$`97.5%`, -0.01, tolerance = tol)
})

# RE heterogeneity SD
bcg_unadj_sd <- as.data.frame(summary(bcg_fit_unadj, pars = "tau"))

test_that("Unadjusted heterogeneity SD", {
  expect_equivalent(bcg_unadj_sd$`50%`, 0.649, tolerance = tol)
  expect_equivalent(bcg_unadj_sd$sd, 0.202, tolerance = tol)
  expect_equivalent(bcg_unadj_sd$`2.5%`, 0.39, tolerance = tol)
  expect_equivalent(bcg_unadj_sd$`97.5%`, 1.17, tolerance = tol)
})

bcg_lat_sd <- as.data.frame(summary(bcg_fit_lat, pars = "tau"))

test_that("Regression heterogeneity SD", {
  expect_equivalent(bcg_lat_sd$`50%`, 0.272, tolerance = tol)
  expect_equivalent(bcg_lat_sd$sd, 0.188, tolerance = tol)
  expect_equivalent(bcg_lat_sd$`2.5%`, 0.03, tolerance = tol)
  expect_equivalent(bcg_lat_sd$`97.5%`, 0.75, tolerance = tol)
})

# DIC
test_that("Unadjusted DIC", {
  expect_equivalent(bcg_dic_unadj$resdev, 26.1, tolerance = tol_dic)
  expect_equivalent(bcg_dic_unadj$pd, 23.5, tolerance = tol_dic)
  expect_equivalent(bcg_dic_unadj$dic, 49.6, tolerance = tol_dic)
})

test_that("Regression DIC", {
  expect_equivalent(bcg_dic_lat$resdev, 30.4, tolerance = tol_dic)
  expect_equivalent(bcg_dic_lat$pd, 21.1, tolerance = tol_dic)
  expect_equivalent(bcg_dic_lat$dic, 51.5, tolerance = tol_dic)
})

test_that("Relative effects and predict work with data.frame", {
  new <- tibble::tibble(latitude = seq(10, 50, by = 10), label = paste0(latitude, "\u00B0 latitude"))
  expect_identical(relative_effects(bcg_fit_lat, newdata = new, study = label),
                   relative_effects(bcg_fit_lat, newdata = as.data.frame(new), study = label))
  # For predict() we need to account for the random baseline sample
  # expect_identical(withr::with_seed(predict(bcg_fit_lat, newdata = new, study = label,
  #                                           baseline = distr(qnorm, mean = -2, sd = 0.1)), seed = 1234),
  #                  withr::with_seed(predict(bcg_fit_lat, newdata = as.data.frame(new), study = label,
  #                                           baseline = distr(qnorm, mean = -2, sd = 0.1)), seed = 1234))
  qcons <- function(p, cons = 0) {cons}
  expect_identical(predict(bcg_fit_lat, newdata = new, study = label,
                           baseline = distr(qcons, -2)),
                   predict(bcg_fit_lat, newdata = as.data.frame(new), study = label,
                           baseline = distr(qcons, -2)))
})