tests/testthat/test-example_plaque_psoriasis.R

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

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

skip_on_cran()



params <-
list(run_tests = FALSE, eval_multinomial = FALSE)

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


## ----setup----------------------------------------------------------------------------------------
library(multinma)
library(dplyr)      # dplyr and tidyr for data manipulation
library(tidyr)
library(ggplot2)    # ggplot2 for plotting covariate distributions

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

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


## -------------------------------------------------------------------------------------------------
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)


## -------------------------------------------------------------------------------------------------
pso_ipd <- pso_ipd %>% 
  mutate(# Variable transformations
         bsa = bsa / 100,
         prevsys = as.numeric(prevsys),
         psa = as.numeric(psa),
         weight = weight / 10,
         durnpso = durnpso / 10,
         # Treatment classes
         trtclass = case_when(trtn == 1 ~ "Placebo",
                              trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                              trtn == 4 ~ "TNFa blocker"),
         # Check complete cases for covariates of interest
         complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>% 
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                              trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                              trtn == 4 ~ "TNFa blocker")
  )


## -------------------------------------------------------------------------------------------------
sum(!pso_ipd$complete)
mean(!pso_ipd$complete)


## -------------------------------------------------------------------------------------------------
pso_ipd <- filter(pso_ipd, complete)


## -------------------------------------------------------------------------------------------------
pso_net <- combine_network(
  set_ipd(pso_ipd, 
          study = studyc, 
          trt = trtc, 
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd, 
              study = studyc, 
              trt = trtc, 
              r = pasi75_r, 
              n = pasi75_n,
              trt_class = trtclass)
)

pso_net


## ----pso_network_plot, fig.width=8, fig.height=6--------------------------------------------------
plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE) + 
  ggplot2::theme(legend.position = "bottom", legend.box = "vertical")


## ----pso_covariate_plot---------------------------------------------------------------------------
# Get mean and sd of covariates in each study
ipd_summary <- pso_ipd %>% 
  group_by(studyc) %>% 
  summarise_at(vars(weight, durnpso, bsa), list(mean = mean, sd = sd, min = min, max = max)) %>% 
  pivot_longer(weight_mean:bsa_max, names_sep = "_", names_to = c("covariate", ".value")) %>% 
  # Assign distributions
  mutate(dist = recode(covariate,
                       bsa = "dlogitnorm",
                       durnpso = "dgamma",
                       weight = "dgamma")) %>% 
  # Compute density curves
  group_by(studyc, covariate) %>% 
  mutate(value = if_else(dist == "dlogitnorm",
                         list(seq(0, 1, length.out = 101)),
                         list(seq(min*0.8, max*1.2, length.out = 101)))) %>% 
  unnest(cols = value) %>% 
  mutate(dens = eval(call(first(dist), x = value, mean = first(mean), sd = first(sd))))

# Plot histograms and assumed densities
pso_ipd %>% 
  pivot_longer(c(weight, durnpso, bsa), names_to = "covariate", values_to = "value") %>% 
ggplot(aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), 
                 binwidth = function(x) diff(range(x)) / nclass.Sturges(x),
                 boundary = 0,
                 fill = "grey50") +
  geom_line(aes(y = dens), data = ipd_summary,
            colour = "darkred", linewidth = 0.5) +
  facet_wrap(~studyc + covariate, scales = "free", ncol = 3) +
  theme_multinma()


## -------------------------------------------------------------------------------------------------
pso_net <- add_integration(pso_net,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  n_int = 1000
)


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


## -------------------------------------------------------------------------------------------------
pso_fit_FE <- nma(pso_net, 
                  trt_effects = "fixed",
                  link = "probit", 
                  likelihood = "bernoulli2",
                  regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
                  class_interactions = "common",
                  prior_intercept = normal(scale = 10),
                  prior_trt = normal(scale = 10),
                  prior_reg = normal(scale = 10),
                  init_r = 0.1,
                  QR = TRUE)


## -------------------------------------------------------------------------------------------------
print(pso_fit_FE)


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


## ----pso_FE_pp_plot, fig.width=8, fig.height=6----------------------------------------------------
plot_prior_posterior(pso_fit_FE, prior = c("intercept", "trt", "reg"))


## ----pso_FE_cumint--------------------------------------------------------------------------------
plot_integration_error(pso_fit_FE)


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


## ---- eval=!params$run_tests----------------------------------------------------------------------
## pso_fit_RE <- nma(pso_net,
##                   trt_effects = "random",
##                   link = "probit",
##                   likelihood = "bernoulli2",
##                   regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
##                   class_interactions = "common",
##                   prior_intercept = normal(scale = 10),
##                   prior_trt = normal(scale = 10),
##                   prior_reg = normal(scale = 10),
##                   prior_het = half_normal(scale = 2.5),
##                   init_r = 0.1,
##                   QR = TRUE)


## ---- eval=!params$run_tests----------------------------------------------------------------------
## print(pso_fit_RE)


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


## ----pso_RE_pairs, eval=!params$run_tests---------------------------------------------------------
## pairs(pso_fit_RE, pars = c("delta[UNCOVER-2: ETN]", "d[ETN]", "tau", "lp__"))


## ----pso_RE_pp_plot, eval=!params$run_tests, fig.width=8, fig.height=6----------------------------
## plot_prior_posterior(pso_fit_RE, prior = c("intercept", "trt", "reg", "het"))


## ----pso_RE_cumint, eval=!params$run_tests--------------------------------------------------------
## plot_integration_error(pso_fit_RE)


## ---- eval=!params$run_tests----------------------------------------------------------------------
## (pso_dic_FE <- dic(pso_fit_FE))
## (pso_dic_RE <- dic(pso_fit_RE))

## ---- eval=params$run_tests, echo=FALSE-----------------------------------------------------------
(pso_dic_FE <- dic(pso_fit_FE))


## -------------------------------------------------------------------------------------------------
plot(pso_fit_FE,
     pars = "beta",
     stat = "halfeye",
     ref_line = 0)


## ----pso_releff_FE--------------------------------------------------------------------------------
(pso_releff_FE <- relative_effects(pso_fit_FE))
plot(pso_releff_FE, ref_line = 0)


## ----pso_pred_FE----------------------------------------------------------------------------------
(pso_pred_FE <- predict(pso_fit_FE, type = "response"))
plot(pso_pred_FE, ref_line = c(0, 1))


## ----pso_ranks_FE---------------------------------------------------------------------------------
(pso_ranks_FE <- posterior_ranks(pso_fit_FE, lower_better = FALSE))
plot(pso_ranks_FE)


## ----pso_rankprobs_FE-----------------------------------------------------------------------------
(pso_rankprobs_FE <- posterior_rank_probs(pso_fit_FE, lower_better = FALSE))
plot(pso_rankprobs_FE)


## ----pso_cumrankprobs_FE--------------------------------------------------------------------------
(pso_cumrankprobs_FE <- posterior_rank_probs(pso_fit_FE, lower_better = FALSE, cumulative = TRUE))
plot(pso_cumrankprobs_FE)


## -------------------------------------------------------------------------------------------------
new_agd_means <- tibble(
  bsa = 0.6,
  prevsys = 0.1,
  psa = 0.2,
  weight = 10,
  durnpso = 3)


## ----pso_releff_FE_new----------------------------------------------------------------------------
(pso_releff_FE_new <- relative_effects(pso_fit_FE, newdata = new_agd_means))
plot(pso_releff_FE_new, ref_line = 0)


## -------------------------------------------------------------------------------------------------
new_agd_int <- tibble(
  bsa_mean = 0.6,
  bsa_sd = 0.3,
  prevsys = 0.1,
  psa = 0.2,
  weight_mean = 10,
  weight_sd = 1,
  durnpso_mean = 3,
  durnpso_sd = 1
)


## -------------------------------------------------------------------------------------------------
new_agd_int <- add_integration(new_agd_int,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  cor = pso_net$int_cor,
  n_int = 1000)


## ----pso_pred_FE_new------------------------------------------------------------------------------
(pso_pred_FE_new <- predict(pso_fit_FE, 
                            type = "response",
                            newdata = new_agd_int,
                            baseline = distr(qnorm, -1.75, 0.08)))
plot(pso_pred_FE_new, ref_line = c(0, 1))


## ---- eval=!params$run_tests----------------------------------------------------------------------
## # IPD studies
## pso_ipd <- plaque_psoriasis_ipd %>%
##   mutate(
##     # Variable transformations
##     bsa = bsa / 100,
##     weight = weight / 10,
##     durnpso = durnpso / 10,
##     prevsys = as.numeric(prevsys),
##     psa = as.numeric(psa),
##     # Treatment classes
##     trtclass = case_when(trtn == 1 ~ "Placebo",
##                          trtn %in% c(2, 3, 5, 6) ~ "IL-17 blocker",
##                          trtn == 4 ~ "TNFa blocker",
##                          trtn == 7 ~ "IL-12/23 blocker"),
##     # Check complete cases for covariates of interest
##     is_complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
##   ) %>%
##   arrange(studyc, trtn)
## 
## # AgD studies
## pso_agd <- plaque_psoriasis_agd %>%
##   mutate(
##     # Variable transformations
##     bsa_mean = bsa_mean / 100,
##     bsa_sd = bsa_sd / 100,
##     weight_mean = weight_mean / 10,
##     weight_sd = weight_sd / 10,
##     durnpso_mean = durnpso_mean / 10,
##     durnpso_sd = durnpso_sd / 10,
##     prevsys = prevsys / 100,
##     psa = psa / 100,
##     # Treatment classes
##     trtclass = case_when(trtn == 1 ~ "Placebo",
##                          trtn %in% c(2, 3, 5, 6) ~ "IL-17 blocker",
##                          trtn == 4 ~ "TNFa blocker",
##                          trtn == 7 ~ "IL-12/23 blocker")
##     ) %>%
##   arrange(studyc, trtn)


## ---- eval=!params$run_tests----------------------------------------------------------------------
## pso_ipd %>%
##   group_by(studyc) %>%
##   summarise(n_total = n(),
##             n_missing = sum(!is_complete),
##             pct_missing = mean(!is_complete) * 100)
## 
## pso_ipd <- filter(pso_ipd, is_complete)


## ---- eval=!params$run_tests----------------------------------------------------------------------
## pso_net <- combine_network(
##   set_ipd(pso_ipd,
##     study = studyc,
##     trt = trtc,
##     r = multi(r0 = 1,
##               PASI75 = pasi75,
##               PASI90 = pasi90,
##               PASI100 = pasi100,
##               type = "ordered", inclusive = TRUE),
##     trt_class = trtclass),
##   set_agd_arm(pso_agd,
##     study = studyc,
##     trt = trtc,
##     r = multi(r0 = pasi75_n,
##               PASI75 = pasi75_r,
##               PASI90 = pasi90_r,
##               PASI100 = pasi100_r,
##               type = "ordered", inclusive = TRUE),
##     trt_class = trtclass)
## )
## 
## pso_net


## ----pso-full-network, eval=!params$run_tests, fig.width=8, fig.height=6--------------------------
## class_pal <- c("#D95F02", "#7570B3", "#E7298A", "#E6AB02")
## 
## plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE, nudge = 0.1) +
##   ggraph::scale_edge_colour_manual("Data",
##                                    values = c(AgD = "#113259", IPD = "#55A480"),
##                                    guide = guide_legend(override.aes = list(edge_width = 2))) +
##   scale_fill_manual("Treatment class",
##                     values = class_pal,
##                     aesthetics = c("fill", "colour"),
##                     guide = guide_legend(override.aes = list(size = 2)))


## ---- eval=!params$run_tests----------------------------------------------------------------------
## pso_net <- add_integration(pso_net,
##   durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
##   prevsys = distr(qbern, prob = prevsys),
##   bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
##   weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
##   psa = distr(qbern, prob = psa),
##   n_int = 1000)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_fit_FE <- nma(pso_net,
##                   trt_effects = "fixed",
##                   link = "probit",
##                   regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
##                   class_interactions = "common",
##                   prior_intercept = normal(scale = 10),
##                   prior_trt = normal(scale = 10),
##                   prior_reg = normal(scale = 10),
##                   prior_aux = flat(),
##                   QR = TRUE,
##                   init_r = 0.5)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_fit_FE


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_fit_RE <- nma(pso_net,
##                   trt_effects = "random",
##                   link = "probit",
##                   regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
##                   class_interactions = "common",
##                   prior_intercept = normal(scale = 10),
##                   prior_trt = normal(scale = 10),
##                   prior_reg = normal(scale = 10),
##                   prior_aux = flat(),
##                   prior_het = half_normal(scale = 2.5),
##                   QR = TRUE,
##                   init_r = 0.5)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_fit_RE


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_dic_FE <- dic(pso_fit_FE))


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_dic_RE <- dic(pso_fit_RE))


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_fit_UME <- nma(pso_net,
##                    trt_effects = "fixed",
##                    consistency = "ume",
##                    link = "probit",
##                    regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
##                    class_interactions = "common",
##                    prior_intercept = normal(scale = 10),
##                    prior_trt = normal(scale = 10),
##                    prior_reg = normal(scale = 10),
##                    prior_aux = flat(),
##                    QR = TRUE,
##                    init_r = 0.5)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_fit_UME


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_dic_FE


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_dic_UME <- dic(pso_fit_UME))


## ----pso-full-dev-dev, eval=!params$run_tests && params$eval_multinomial--------------------------
## plot(pso_dic_FE, pso_dic_UME, show_uncertainty = FALSE) +
##   xlab("Residual deviance - consistency model") +
##   ylab("Residual deviance - inconsistency (UME) model")


## ---- eval=!params$run_tests----------------------------------------------------------------------
## data.frame(classes = pso_net$classes, treatments = pso_net$treatments)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## noSEM_mods <- list(
##   durnpso = ~(prevsys + bsa + weight + psa)*.trtclass + durnpso*.trt,
##   prevsys = ~(durnpso + bsa + weight + psa)*.trtclass + prevsys*.trt,
##   bsa = ~(durnpso + prevsys + weight + psa)*.trtclass + bsa*.trt,
##   weight = ~(durnpso + prevsys + bsa + psa)*.trtclass + weight*.trt,
##   psa = ~(durnpso + prevsys + bsa + weight)*.trtclass + psa*.trt
##   )
## 
## noSEM_fits <- noSEM_mods
## 
## for (m in 1:length(noSEM_mods)) {
##   cat("Fitting model with independent interactions for", names(noSEM_mods)[m], "\n")
## 
##   noSEM_fits[[m]] <-
##     nma(pso_net,
##         trt_effects = "fixed",
##         link = "probit",
##         regression = noSEM_mods[[m]],
##         class_interactions = "independent",
##         prior_intercept = normal(scale = 10),
##         prior_trt = normal(scale = 10),
##         prior_reg = normal(scale = 10),
##         prior_aux = flat(),
##         QR = TRUE,
##         init_r = 0.5,
##         # Using save_warmup = FALSE reduces memory footprint when
##         # fitting many models in one session
##         save_warmup = FALSE)
## }


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## pso_dic_FE


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## lapply(noSEM_fits, dic)


## ----pso-full-relaxing-SEM, eval=!params$run_tests && params$eval_multinomial, fig.height=6, fig.width=8----
## library(purrr)
## library(stringr)
## library(forcats)
## 
## # Extract draws from relaxed models
## imap_dfr(noSEM_fits,
##         ~as_tibble(as.matrix(.x, pars = "beta")) %>%
##            pivot_longer(cols = everything(), names_to = "parameter", values_to = "value") %>%
##            filter(str_detect(parameter, paste0("(IXE|SEC).+:", .y))) %>%
##            mutate(model = .y)) %>%
## 
##   # Add in draws from the original model
##   bind_rows(
##     as_tibble(as.matrix(pso_fit_FE, pars = "beta")) %>%
##     pivot_longer(cols = everything(), names_to = "parameter", values_to = "value") %>%
##     filter(str_detect(parameter, ":.+IL\\-17 blocker")) %>%
##     mutate(model = "all")
##   ) %>%
## 
##   mutate(
##     # Rescale BSA to per 10%
##     value = if_else(str_detect(parameter, "bsa"), value / 10, value),
##     # Create labels
##     covariate = str_extract(parameter, "durnpso|prevsys|bsa|weight|psa"),
##     covariatef = recode_factor(covariate,
##                                durnpso = "Duration of psoriasis, per 10 years",
##                                prevsys = "Previous systemic use",
##                                bsa = "Body surface area, per 10%",
##                                weight = "Weight, per 10 kg",
##                                psa = "Psoriatic arthritis"),
##     treatment = str_remove(str_extract(parameter, "\\.trt(class)?.+?(?=[\\]:])"),
##                            "\\.trt(class)?"),
##     Interactions = fct_collapse(factor(model),
##                                 Common = "all",
##                                 other_level = "Independent")) %>%
## 
## # Plot
## ggplot(aes(x = value, y = fct_rev(treatment), colour = Interactions, fill = Interactions)) +
##   geom_vline(xintercept = 0, colour = "grey70") +
##   ggdist::stat_halfeye(normalize = "panels", slab_alpha = 0.3, .width = c(0, 0.95)) +
##   facet_wrap("covariatef", scales = "free") +
##   xlab("Interaction effect (SMD)") +
##   ylab("Treatment / Class") +
##   scale_colour_manual(values = c(Common = "#7B3294", Independent = "#91D388"),
##                       aesthetics = c("colour", "fill")) +
##   theme_multinma() +
##   theme(legend.position = c(0.85, 0.2))


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_releff_FE <- relative_effects(pso_fit_FE))


## ----pso-full-releff, eval=!params$run_tests && params$eval_multinomial, fig.height=8, fig.width=6----
## plot(pso_releff_FE, ref_line = 0)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_pred_FE <- predict(pso_fit_FE, type = "response"))


## ----pso-full-pred, eval=!params$run_tests && params$eval_multinomial, fig.width=6, fig.height=8----
## plot(pso_pred_FE, ref_line = c(0, 1))


## ---- eval=!params$run_tests----------------------------------------------------------------------
## new_agd_means <- tibble::tribble(
##              ~study, ~covariate,  ~mean,   ~sd,
##           "PsoBest",      "bsa",     24,  20.5,
##           "PsoBest",  "durnpso",   18.2,  14.1,
##           "PsoBest",  "prevsys",   0.54,    NA,
##           "PsoBest",      "psa",  0.207,    NA,
##           "PsoBest",   "weight",     85,  19.1,
##          "PROSPECT",      "bsa",   18.7,  18.4,
##          "PROSPECT",  "durnpso",   19.6,  13.5,
##          "PROSPECT",  "prevsys", 0.9095,    NA,
##          "PROSPECT",      "psa",  0.202,    NA,
##          "PROSPECT",   "weight",   87.5,  20.3,
##   "Chiricozzi 2019",      "bsa",     23, 16.79,
##   "Chiricozzi 2019",  "durnpso",  16.93, 10.82,
##   "Chiricozzi 2019",  "prevsys", 0.9061,    NA,
##   "Chiricozzi 2019",      "psa", 0.2152,    NA,
##   "Chiricozzi 2019",   "weight",   78.3, 15.87
##   ) %>%
##   # Tidy up
##   pivot_wider(id_cols = study,
##               names_from = covariate,
##               values_from = c(mean, sd),
##               names_glue = "{covariate}_{.value}") %>%
##   # Rescale as per analysis
##   transmute(study,
##             bsa_mean = bsa_mean / 100,
##             bsa_sd = bsa_sd / 100,
##             weight_mean = weight_mean / 10,
##             weight_sd = weight_sd / 10,
##             durnpso_mean = durnpso_mean / 10,
##             durnpso_sd = durnpso_sd / 10,
##             prevsys = prevsys_mean,
##             psa = psa_mean)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_releff_FE_new <- relative_effects(pso_fit_FE,
##                                        newdata = transmute(new_agd_means,
##                                                            study,
##                                                            bsa = bsa_mean,
##                                                            weight = weight_mean,
##                                                            durnpso = durnpso_mean,
##                                                            prevsys,
##                                                            psa),
##                                        study = study))


## ----pso-full-releff-new, eval=!params$run_tests && params$eval_multinomial-----------------------
## plot(pso_releff_FE_new, ref_line = 0) + facet_wrap("Study")


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## new_agd_int <- add_integration(filter(new_agd_means, study != "PsoBest"),
##                                durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
##                                prevsys = distr(qbern, prob = prevsys),
##                                bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
##                                weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
##                                psa = distr(qbern, prob = psa),
##                                n_int = 1000,
##                                cor = pso_net$int_cor)


## ---- eval=!params$run_tests && params$eval_multinomial-------------------------------------------
## (pso_pred_FE_new <- predict(pso_fit_FE,
##         type = "response",
##         newdata = new_agd_int,
##         study = study,
##         baseline = list(PROSPECT = distr(qbeta, 1156, 1509-1156),
##                         "Chiricozzi 2019" = distr(qbeta, 243, 330-243)),
##         baseline_type = "response",
##         baseline_level = "aggregate",
##         trt_ref = "SEC_300"))


## ----pso-full-pred-new, eval=!params$run_tests && params$eval_multinomial-------------------------
## plot(pso_pred_FE_new, ref_line = c(0, 1)) +
##   facet_grid(rows = "Study") +
##   aes(colour = Category) +
##   scale_colour_brewer(palette = "Blues")


## ----pso_tests, include=FALSE, eval=params$run_tests----------------------------------------------
library(testthat)
library(dplyr)

tol <- 0.05
tol_dic <- 0.1

# FE model parameters
test_fe <- tribble(
  ~parameter                            , ~mean, ~sd , ~`2.5%`, ~`50%`, ~`97.5%`,
  "beta[durnpso]"                       , 0.05 , 0.06, -0.08  , 0.05  , 0.17    ,
  "beta[prevsys]"                       , -0.13, 0.16, -0.44  , -0.13 , 0.17    ,
  "beta[bsa]"                           , -0.06, 0.45, -0.98  , -0.05 , 0.78    ,
  "beta[weight]"                        , 0.04 , 0.03, -0.02  , 0.04  , 0.10    ,
  "beta[psa]"                           , -0.08, 0.17, -0.42  , -0.08 , 0.25    ,
  "beta[durnpso:.trtclassTNFa blocker]", -0.03, 0.08, -0.18  , -0.03 , 0.12    ,
  "beta[durnpso:.trtclassIL blocker]"  , -0.01, 0.07, -0.15  , -0.01 , 0.12    ,
  "beta[prevsys:.trtclassTNFa blocker]", 0.19 , 0.19, -0.19  , 0.19  , 0.56    ,
  "beta[prevsys:.trtclassIL blocker]"  , 0.06 , 0.17, -0.28  , 0.06  , 0.40    ,
  "beta[bsa:.trtclassTNFa blocker]"    , 0.05 , 0.52, -0.95  , 0.04  , 1.10    ,
  "beta[bsa:.trtclassIL blocker]"      , 0.29 , 0.49, -0.64  , 0.29  , 1.27    ,
  "beta[weight:.trtclassTNFa blocker]" , -0.17, 0.04, -0.23  , -0.17 , -0.10   ,
  "beta[weight:.trtclassIL blocker]"   , -0.10, 0.03, -0.16  , -0.10 , -0.03   ,
  "beta[psa:.trtclassTNFa blocker]"    , -0.06, 0.21, -0.46  , -0.06 , 0.37    ,
  "beta[psa:.trtclassIL blocker]"      , 0.00 , 0.18, -0.35  , 0.00  , 0.37    ,
  "d[ETN]"                              , 1.55 , 0.08, 1.39   , 1.55  , 1.72    ,
  "d[IXE_Q2W]"                          , 2.95 , 0.09, 2.79   , 2.95  , 3.13    ,
  "d[IXE_Q4W]"                          , 2.54 , 0.08, 2.38   , 2.54  , 2.71    ,
  "d[SEC_150]"                          , 2.14 , 0.11, 1.93   , 2.14  , 2.37    ,
  "d[SEC_300]"                          , 2.45 , 0.12, 2.22   , 2.45  , 2.69    )

summary_fe <- summary(pso_fit_FE, pars = c("beta", "d")) %>% 
  as_tibble() %>% 
  select(parameter, mean, sd, `2.5%`, `50%`, `97.5%`)

test_that("FE model parameters", {
  expect_equal(summary_fe, test_fe, tolerance = tol, check.attributes = FALSE)
})

# FE DIC
test_that("FE DIC", {
  expect_equivalent(pso_dic_FE$resdev, 3129.4, tolerance = tol_dic)
  expect_equivalent(pso_dic_FE$pd, 24.1, tolerance = tol_dic)
  expect_equivalent(pso_dic_FE$dic, 3153.5, tolerance = tol_dic)
})

# Population average relative effects in target population
test_releff_fe_new <- tribble(
  ~parameter         , ~mean, ~sd , ~`2.5%`, ~`50%`, ~`97.5%`,
  "d[New 1: ETN]"    , 1.25 , 0.24, 0.81 , 1.24, 1.73  ,
  "d[New 1: IXE_Q2W]", 2.89 , 0.23, 2.46 , 2.88, 3.36  ,
  "d[New 1: IXE_Q4W]", 2.48 , 0.23, 2.05 , 2.47, 2.94  ,
  "d[New 1: SEC_150]", 2.08 , 0.23, 1.64 , 2.07, 2.55  ,
  "d[New 1: SEC_300]", 2.39 , 0.23, 1.95 , 2.38, 2.86  )

summary_releff_fe_new <- pso_releff_FE_new %>% 
  as_tibble() %>% 
  select(parameter, mean, sd, `2.5%`, `50%`, `97.5%`)

test_that("FE PATE in target population", {
  expect_equal(summary_releff_fe_new, test_releff_fe_new, tolerance = tol, check.attributes = FALSE)
})

# Check construction of all contrasts
pso_releff_FE_all_contr <- relative_effects(pso_fit_FE, all_contrasts = TRUE)

# Reconstruct from basic contrasts in each study
dk <- function(study, trt, sims = pso_releff_FE$sims) {
  if (trt == "PBO") return(0)
  else return(sims[ , , paste0("d[", study, ": ", trt, "]"), drop = FALSE])
}

test_all_contr <- tibble(
  contr = pso_releff_FE_all_contr$summary$parameter,
  .study = factor(stringr::str_extract(contr, "(?<=\\[)(.+)(?=:)")),
  .trtb = factor(stringr::str_extract(contr, "(?<=\\: )(.+)(?= vs\\.)"), levels = levels(pso_net$treatments)),
  .trta = factor(stringr::str_extract(contr, "(?<=vs\\. )(.+)(?=\\])"), levels = levels(pso_net$treatments))
) %>% 
  rowwise() %>% 
  mutate(as_tibble(multinma:::summary.mcmc_array(dk(.study, .trtb) - dk(.study, .trta)))) %>% 
  select(.study, .trtb, .trta, parameter = contr, mean:Rhat)

test_that("Construction of all contrasts is correct", {
  ntrt <- nlevels(pso_net$treatments)
  nstudy <- nlevels(test_all_contr$.study)
  expect_equal(nrow(pso_releff_FE_all_contr$summary), nstudy * ntrt * (ntrt - 1) / 2)
  expect_equal(pso_releff_FE_all_contr$summary, test_all_contr, check.attributes = FALSE)
})

# Check all contrasts in target population
pso_releff_FE_all_contr_new <- relative_effects(pso_fit_FE, newdata = new_agd_means, all_contrasts = TRUE)

test_all_contr_new <- tibble(
  contr = pso_releff_FE_all_contr_new$summary$parameter,
  .study = factor(stringr::str_extract(contr, "(?<=\\[)(.+)(?=:)")),
  .trtb = factor(stringr::str_extract(contr, "(?<=\\: )(.+)(?= vs\\.)"), levels = levels(pso_net$treatments)),
  .trta = factor(stringr::str_extract(contr, "(?<=vs\\. )(.+)(?=\\])"), levels = levels(pso_net$treatments))
) %>% 
  rowwise() %>% 
  mutate(as_tibble(multinma:::summary.mcmc_array(dk(.study, .trtb, pso_releff_FE_new$sims) - dk(.study, .trta, pso_releff_FE_new$sims)))) %>% 
  select(.study, .trtb, .trta, parameter = contr, mean:Rhat)

test_that("Construction of all contrasts in target population is correct", {
  ntrt <- nlevels(pso_net$treatments)
  nstudy <- nlevels(test_all_contr_new$.study)
  expect_equal(nrow(pso_releff_FE_all_contr_new$summary), nstudy * ntrt * (ntrt - 1) / 2)
  expect_equal(pso_releff_FE_all_contr_new$summary, test_all_contr_new, check.attributes = FALSE)
})

# Population average response probabilities in target population
test_pred_fe_new <- tribble(
  ~parameter            , ~mean, ~sd , ~`2.5%`, ~`50%`, ~`97.5%`,
  "pred[New 1: PBO]"    , 0.06 , 0.03, 0.02   , 0.06  , 0.12    ,
  "pred[New 1: ETN]"    , 0.37 , 0.06, 0.26   , 0.37  , 0.48    ,
  "pred[New 1: IXE_Q2W]", 0.90 , 0.03, 0.84   , 0.90  , 0.94    , 
  "pred[New 1: IXE_Q4W]", 0.81 , 0.04, 0.72   , 0.81  , 0.88    ,
  "pred[New 1: SEC_150]", 0.68 , 0.06, 0.57   , 0.68  , 0.78    ,
  "pred[New 1: SEC_300]", 0.78 , 0.05, 0.68   , 0.78  , 0.86    )

summary_pred_fe_new <- pso_pred_FE_new %>% 
  as_tibble() %>% 
  select(parameter, mean, sd, `2.5%`, `50%`, `97.5%`)

test_that("FE average response probabilities in target population", {
  expect_equal(summary_pred_fe_new, test_pred_fe_new, tolerance = tol, check.attributes = FALSE)
})

test_that("Robust to custom options(contrasts) settings", {
  withr::with_options(list(contrasts = c(ordered = "contr.SAS", unordered = "contr.SAS")), {
    pso_fit_FE_SAS <- nma(pso_net, 
                          trt_effects = "fixed",
                          link = "probit", 
                          likelihood = "bernoulli2",
                          regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
                          class_interactions = "common",
                          prior_intercept = normal(scale = 10),
                          prior_trt = normal(scale = 10),
                          prior_reg = normal(scale = 10),
                          init_r = 0.1,
                          QR = TRUE)
    
    pso_fit_FE_SAS_releff <- as_tibble(relative_effects(pso_fit_FE_SAS))[, c("parameter", "mean", "sd")]
  })
  
  expect_equal(pso_fit_FE_SAS_releff,
               as_tibble(relative_effects(pso_fit_FE))[, c("parameter", "mean", "sd")],
               tolerance = tol)
})