tests/testthat/test-example_plaque_psoriasis.R
In dmphillippo/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 = 64
)
## -------------------------------------------------------------------------------------------------
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"))
## -------------------------------------------------------------------------------------------------
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"))
## ----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 = 64)
## ----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")) +
## scale_fill_manual("Treatment class",
## values = class_pal,
## aesthetics = c("fill", "colour")) +
## guides(edge_colour = guide_legend(override.aes = list(edge_width = 2)),
## fill = 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 = 64)
## ----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 = 64,
## 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",
## baseline_trt = "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)
})
test_that("Integration checks with int_check throw expected warnings", {
local_edition(3)
pso_net_lownint <-
suppressWarnings(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 = 2))
expect_warning(
expect_warning(
expect_warning(
nma(pso_net_lownint,
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,
),
class = "int_check_rhat"),
class = "int_check_essb"),
class = "int_check_esst")
})
# Force clean up
rm(list = ls())
gc()
dmphillippo/multinma documentation built on April 12, 2025, 11:41 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# 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 = 64
)
## -------------------------------------------------------------------------------------------------
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"))
## -------------------------------------------------------------------------------------------------
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"))
## ----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 = 64)
## ----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")) +
## scale_fill_manual("Treatment class",
## values = class_pal,
## aesthetics = c("fill", "colour")) +
## guides(edge_colour = guide_legend(override.aes = list(edge_width = 2)),
## fill = 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 = 64)
## ----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 = 64,
## 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",
## baseline_trt = "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)
})
test_that("Integration checks with int_check throw expected warnings", {
local_edition(3)
pso_net_lownint <-
suppressWarnings(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 = 2))
expect_warning(
expect_warning(
expect_warning(
nma(pso_net_lownint,
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,
),
class = "int_check_rhat"),
class = "int_check_essb"),
class = "int_check_esst")
})
# Force clean up
rm(list = ls())
gc()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.