tests/testthat/test-example_hta_psoriasis.R

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

# Generated by vignette example_hta_psoriasis.Rmd: do not edit by hand
# Instead edit example_hta_psoriasis.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)
set.seed(65498431)


## -------------------------------------------------------------------------------------------------
head(hta_psoriasis)


## -------------------------------------------------------------------------------------------------
dplyr::filter(hta_psoriasis, studyc %in% c("Elewski", "Gordon", "ACD2058g", "Altmeyer"))


## -------------------------------------------------------------------------------------------------
pso_net <- set_agd_arm(hta_psoriasis, 
                       study = paste(studyc, year), 
                       trt = trtc, 
                       r = multi(r0 = sample_size - rowSums(cbind(PASI50, PASI75, PASI90), na.rm = TRUE), 
                                 PASI50, PASI75, PASI90,
                                 inclusive = FALSE, 
                                 type = "ordered"))
pso_net


## ----hta_psoriasis_network_plot-------------------------------------------------------------------
plot(pso_net, weight_edges = TRUE, weight_nodes = TRUE) + 
  # Nudge the legend over
  ggplot2::theme(legend.box.spacing = ggplot2::unit(0.75, "in"),
                 plot.margin = ggplot2::margin(0.1, 0, 0.1, 0.75, "in"))


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


## -------------------------------------------------------------------------------------------------
pso_fit_FE <- nma(pso_net, 
                  trt_effects = "fixed",
                  link = "probit",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 10),
                  prior_aux = flat())


## -------------------------------------------------------------------------------------------------
pso_fit_FE


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


## ----pso_FE_pp_plot-------------------------------------------------------------------------------
plot_prior_posterior(pso_fit_FE)


## ----pso_FE_pp_cutpoint_plot----------------------------------------------------------------------
plot_prior_posterior(pso_fit_FE, prior = "aux")


## -------------------------------------------------------------------------------------------------
summary(normal(scale = 10))
summary(normal(scale = 100))
summary(half_normal(scale = 2.5))


## ----eval = FALSE---------------------------------------------------------------------------------
## pso_fit_RE <- nma(pso_net,
##                   trt_effects = "random",
##                   link = "probit",
##                   prior_intercept = normal(scale = 100),
##                   prior_trt = normal(scale = 10),
##                   prior_aux = flat(),
##                   prior_het = half_normal(scale = 2.5),
##                   adapt_delta = 0.99)

## ----echo = FALSE, warning = FALSE----------------------------------------------------------------
pso_fit_RE <- nowarn_on_ci(nma(pso_net, 
                  trt_effects = "random",
                  link = "probit",
                  prior_intercept = normal(scale = 100),
                  prior_trt = normal(scale = 10),
                  prior_aux = flat(),
                  prior_het = half_normal(scale = 2.5),
                  adapt_delta = 0.99,
                  iter = if (isTRUE(as.logical(Sys.getenv("CI")))) 5000 else 2000,
                  seed = 1713435794))


## -------------------------------------------------------------------------------------------------
pso_fit_RE


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


## ----pso_RE_pp_plot-------------------------------------------------------------------------------
plot_prior_posterior(pso_fit_RE, prior = c("trt", "aux", "het"))


## -------------------------------------------------------------------------------------------------
(dic_FE <- dic(pso_fit_FE))

## -------------------------------------------------------------------------------------------------
(dic_RE <- dic(pso_fit_RE))


## ----pso_FE_resdev_plot---------------------------------------------------------------------------
plot(dic_FE)


## ----pso_RE_resdev_plot---------------------------------------------------------------------------
plot(dic_RE)


## ----pso_pred_FE, fig.height = 2------------------------------------------------------------------
pred_FE <- predict(pso_fit_FE, 
                   baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5), 
                   type = "response")
pred_FE
plot(pred_FE)

## ----pso_pred_RE, fig.height = 2------------------------------------------------------------------
pred_RE <- predict(pso_fit_RE, 
                   baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5), 
                   type = "response")
pred_RE
plot(pred_RE)


## ----pso_pred_FE_beta, fig.height = 2-------------------------------------------------------------
pred_FE_beta <- predict(pso_fit_FE, 
                        baseline = distr(qbeta, 56, 408-56),
                        baseline_type = "response",
                        type = "response")
pred_FE_beta
plot(pred_FE_beta)


## ----pso_pred_RE_beta, fig.height = 2-------------------------------------------------------------
pred_RE_beta <- predict(pso_fit_RE, 
                        baseline = distr(qbeta, 56, 408-56),
                        baseline_type = "response",
                        type = "response")
pred_RE_beta
plot(pred_RE_beta)


## ----pso_pred_RE_colour, fig.height = 3-----------------------------------------------------------
library(ggplot2)
plot(pred_RE, position = position_dodge(width = 0.75)) +
  facet_null() +
  aes(colour = Category) +
  scale_colour_brewer(palette = "Blues")


## ----hta_psoriasis_ranks, fig.height=3------------------------------------------------------------
(pso_ranks <- posterior_ranks(pso_fit_RE, lower_better = FALSE))
plot(pso_ranks)

## ----hta_psoriasis_rankprobs----------------------------------------------------------------------
(pso_rankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE))
plot(pso_rankprobs)

## ----hta_psoriasis_cumrankprobs-------------------------------------------------------------------
(pso_cumrankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE, cumulative = TRUE))
plot(pso_cumrankprobs)


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

tol <- 0.05
tol_dic <- 0.1

trt_codes <- c(
  "1"=  "Supportive care",
  "2" = "Etanercept 25 mg",
  "3"= "Etanercept 50 mg",
  "4" =       "Efalizumab",
  "5"=      "Ciclosporin",
  "6" =         "Fumaderm",
  "7"=       "Infliximab",
  "8" =     "Methotrexate")

# FE relative effects
tsd_fe <- tribble(
  ~contrast, ~est, ~sd, ~median, ~lower, ~upper,
  "d12", -1.51, 0.10, -1.51, -1.70, -1.32,
  "d13", -1.92, 0.10, -1.92, -2.12, -1.72,
  "d14", -1.19, 0.06, -1.19, -1.30, -1.08,
  "d15", -1.92, 0.34, -1.90, -2.62, -1.30,
  "d16", -1.49, 0.49, -1.46, -2.55, -0.63,
  "d17", -2.33, 0.26, -2.33, -2.87, -1.84,
  "d18", -1.61, 0.44, -1.60, -2.50, -0.77) %>% 
  mutate(trt = recode(substr(contrast, 3, 3), !!! trt_codes),
         trt = ordered(trt, levels = levels(pso_net$treatments)),
         # Parameter estimates are opposite sign
         .l = lower, .u = upper,
         est = -est, median = -median,
         lower = -.u, upper = -.l) %>% 
  arrange(trt)

pso_releff_FE <- as.data.frame(relative_effects(pso_fit_FE))
test_that("FE relative effects", {
  expect_equivalent(pso_releff_FE$mean, tsd_fe$est, tolerance = tol)
  expect_equivalent(pso_releff_FE$sd, tsd_fe$sd, tolerance = tol)
  expect_equivalent(pso_releff_FE$`50%`, tsd_fe$median, tolerance = tol)
  expect_equivalent(pso_releff_FE$`2.5%`, tsd_fe$lower, tolerance = tol)
  expect_equivalent(pso_releff_FE$`97.5%`, tsd_fe$upper, tolerance = tol)
})

# FE predicted probabilities
tsd_pred_fe <- tribble(
  ~outcome, ~trt              , ~mean, ~sd , ~median, ~lower, ~upper,
  50      , "Supportive care" , 0.14 , 0.02, 0.14   , 0.10  , 0.18  ,
  50      , "Etanercept 25 mg", 0.66 , 0.05, 0.66   , 0.56  , 0.75  ,
  50      , "Etanercept 50 mg", 0.79 , 0.04, 0.79   , 0.71  , 0.86  ,
  50      , "Efalizumab"      , 0.54 , 0.04, 0.54   , 0.45  , 0.62  ,
  50      , "Ciclosporin"     , 0.78 , 0.10, 0.79   , 0.57  , 0.94  ,
  50      , "Fumaderm"        , 0.64 , 0.16, 0.64   , 0.31  , 0.93  ,
  50      , "Infliximab"      , 0.88 , 0.05, 0.89   , 0.76  , 0.96  ,
  50      , "Methotrexate"    , 0.68 , 0.15, 0.69   , 0.37  , 0.92  ,
  75      , "Supportive care" , 0.03 , 0.01, 0.03   , 0.02  , 0.05  ,
  75      , "Etanercept 25 mg", 0.37 , 0.05, 0.37   , 0.28  , 0.47  ,
  75      , "Etanercept 50 mg", 0.53 , 0.05, 0.53   , 0.42  , 0.63  ,
  75      , "Efalizumab"      , 0.25 , 0.03, 0.25   , 0.19  , 0.33  ,
  75      , "Ciclosporin"     , 0.52 , 0.13, 0.52   , 0.28  , 0.79  ,
  75      , "Fumaderm"        , 0.37 , 0.17, 0.35   , 0.11  , 0.76  ,
  75      , "Infliximab"      , 0.68 , 0.10, 0.68   , 0.48  , 0.85  ,
  75      , "Methotrexate"    , 0.41 , 0.16, 0.40   , 0.14  , 0.75  ,
  90      , "Supportive care" , 0.00 , 0.00, 0.00   , 0.00  , 0.01  ,
  90      , "Etanercept 25 mg", 0.13 , 0.03, 0.13   , 0.08  , 0.19  ,
  90      , "Etanercept 50 mg", 0.23 , 0.04, 0.23   , 0.16  , 0.32  ,
  90      , "Efalizumab"      , 0.07 , 0.02, 0.07   , 0.04  , 0.11  ,
  90      , "Ciclosporin"     , 0.24 , 0.11, 0.22   , 0.08  , 0.49  ,
  90      , "Fumaderm"        , 0.15 , 0.12, 0.11   , 0.02  , 0.46  ,
  90      , "Infliximab"      , 0.38 , 0.10, 0.37   , 0.19  , 0.60  ,
  90      , "Methotrexate"    , 0.17 , 0.11, 0.14   , 0.03  , 0.44  ) %>% 
  mutate(trt = ordered(trt, levels = levels(pso_net$treatments))) %>% 
  arrange(trt, outcome)

pred_FE <- as.data.frame(pred_FE)
test_that("FE predicted probabilities", {
  expect_equivalent(pred_FE$mean, tsd_pred_fe$mean, tolerance = tol)
  expect_equivalent(pred_FE$sd, tsd_pred_fe$sd, tolerance = tol)
  expect_equivalent(pred_FE$`50%`, tsd_pred_fe$median, tolerance = tol)
  expect_equivalent(pred_FE$`2.5%`, tsd_pred_fe$lower, tolerance = tol)
  expect_equivalent(pred_FE$`97.5%`, tsd_pred_fe$upper, tolerance = tol)
})

# FE DIC
test_that("FE DIC", {
  expect_equivalent(dic_FE$resdev, 74.9, tolerance = tol_dic)
  expect_equivalent(dic_FE$pd, 25.0, tolerance = tol_dic)
  expect_equivalent(dic_FE$dic, 99.9, tolerance = tol_dic)
})

# RE relative effects
tsd_re <- tribble(
  ~contrast, ~est, ~sd, ~median, ~lower, ~upper,
  "d12", -1.53, 0.24, -1.52, -2.05, -1.03,
  "d13", -1.93, 0.28, -1.92, -2.51, -1.35,
  "d14", -1.19, 0.18, -1.19, -1.56, -0.81,
  "d15", -2.04, 0.43, -2.00, -3.02, -1.30,
  "d16", -1.49, 0.62, -1.46, -2.81, -0.33,
  "d17", -2.32, 0.38, -2.32, -3.06, -1.55,
  "d18", -1.74, 0.64, -1.70, -3.14, -0.59) %>% 
  mutate(trt = recode(substr(contrast, 3, 3), !!! trt_codes),
         trt = ordered(trt, levels = levels(pso_net$treatments)),
         # Parameter estimates are opposite sign
         .l = lower, .u = upper,
         est = -est, median = -median,
         lower = -.u, upper = -.l) %>% 
  arrange(trt)

pso_releff_RE <- as.data.frame(relative_effects(pso_fit_RE))
test_that("RE relative effects", {
  expect_equivalent(pso_releff_RE$mean, tsd_re$est, tolerance = tol)
  expect_equivalent(pso_releff_RE$sd, tsd_re$sd, tolerance = tol)
  expect_equivalent(pso_releff_RE$`50%`, tsd_re$median, tolerance = tol)
  expect_equivalent(pso_releff_RE$`2.5%`, tsd_re$lower, tolerance = tol)
  expect_equivalent(pso_releff_RE$`97.5%`, tsd_re$upper, tolerance = tol)
})

# RE predicted probabilities
tsd_pred_re <- tribble(
  ~outcome, ~trt              , ~mean, ~sd , ~median, ~lower, ~upper,
  50      , "Supportive care" , 0.14 , 0.02, 0.14   , 0.10  , 0.18,
  50      , "Etanercept 25 mg", 0.66 , 0.09, 0.66   , 0.46  , 0.83,
  50      , "Etanercept 50 mg", 0.79 , 0.08, 0.80   , 0.59  , 0.92,
  50      , "Efalizumab"      , 0.54 , 0.08, 0.54   , 0.38  , 0.69,
  50      , "Ciclosporin"     , 0.81 , 0.10, 0.82   , 0.57  , 0.97,
  50      , "Fumaderm"        , 0.63 , 0.20, 0.64   , 0.22  , 0.96,
  50      , "Infliximab"      , 0.87 , 0.08, 0.89   , 0.67  , 0.98,
  50      , "Methotrexate"    , 0.70 , 0.18, 0.73   , 0.30  , 0.98,
  75      , "Supportive care" , 0.03 , 0.01, 0.03   , 0.02  , 0.05,
  75      , "Etanercept 25 mg", 0.38 , 0.09, 0.37   , 0.20  , 0.58,
  75      , "Etanercept 50 mg", 0.53 , 0.11, 0.53   , 0.30  , 0.75,
  75      , "Efalizumab"      , 0.26 , 0.06, 0.25   , 0.14  , 0.40,
  75      , "Ciclosporin"     , 0.56 , 0.15, 0.56   , 0.28  , 0.88,
  75      , "Fumaderm"        , 0.38 , 0.20, 0.35   , 0.06  , 0.83,
  75      , "Infliximab"      , 0.67 , 0.13, 0.68   , 0.37  , 0.89,
  75      , "Methotrexate"    , 0.46 , 0.21, 0.44   , 0.10  , 0.90,
  90      , "Supportive care" , 0.00 , 0.00, 0.00   , 0.00  , 0.01,
  90      , "Etanercept 25 mg", 0.14 , 0.06, 0.13   , 0.05  , 0.28,
  90      , "Etanercept 50 mg", 0.24 , 0.09, 0.23   , 0.09  , 0.45,
  90      , "Efalizumab"      , 0.07 , 0.03, 0.07   , 0.03  , 0.15,
  90      , "Ciclosporin"     , 0.28 , 0.14, 0.26   , 0.08  , 0.65,
  90      , "Fumaderm"        , 0.16 , 0.15, 0.12   , 0.01  , 0.57,
  90      , "Infliximab"      , 0.37 , 0.14, 0.37   , 0.13  , 0.67,
  90      , "Methotrexate"    , 0.21 , 0.17, 0.17   , 0.02  , 0.69) %>% 
  mutate(trt = ordered(trt, levels = levels(pso_net$treatments))) %>% 
  arrange(trt, outcome)

pred_RE <- as.data.frame(pred_RE)
test_that("RE predicted probabilities", {
  expect_equivalent(pred_RE$mean, tsd_pred_re$mean, tolerance = tol)
  expect_equivalent(pred_RE$sd, tsd_pred_re$sd, tolerance = tol)
  expect_equivalent(pred_RE$`50%`, tsd_pred_re$median, tolerance = tol)
  expect_equivalent(pred_RE$`2.5%`, tsd_pred_re$lower, tolerance = tol)
  expect_equivalent(pred_RE$`97.5%`, tsd_pred_re$upper, tolerance = tol)
})

# Heterogeneity SD
pso_tau <- summary(pso_fit_RE, pars = "tau")

test_that("RE heterogeneity SD", {
  expect_equivalent(pso_tau$summary$mean, 0.31, tolerance = tol)
  expect_equivalent(pso_tau$summary$sd, 0.23, tolerance = tol)
  expect_equivalent(pso_tau$summary$`50%`, 0.26, tolerance = tol)
  skip_on_ci()
  expect_equivalent(pso_tau$summary$`2.5%`, 0.01, tolerance = tol)
  expect_equivalent(pso_tau$summary$`97.5%`, 0.88, tolerance = tol)
})

# RE DIC
test_that("RE DIC", {
  expect_equivalent(dic_RE$resdev, 63.0, tolerance = tol_dic)
  expect_equivalent(dic_RE$pd, 33.3, tolerance = tol_dic)
  expect_equivalent(dic_RE$dic, 96.2, tolerance = tol_dic)
})

# Check predictions with Beta distribution on baseline probability
pred_FE_beta <- as.data.frame(pred_FE_beta)
test_that("FE predicted probabilities (Beta distribution)", {
  expect_equal(pred_FE[c("mean", "sd", "2.5%", "50%", "97.5%")],
               pred_FE_beta[c("mean", "sd", "2.5%", "50%", "97.5%")],
               tolerance = tol)
})

pred_RE_beta <- as.data.frame(pred_RE_beta)
test_that("RE predicted probabilities (Beta distribution)", {
  expect_equal(pred_RE[c("mean", "sd", "2.5%", "50%", "97.5%")],
               pred_RE_beta[c("mean", "sd", "2.5%", "50%", "97.5%")],
               tolerance = tol)
})


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