tests/testthat/test-example_hta_psoriasis.R
In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data
# Generated by vignette example_hta_psoriasis.Rmd: do not edit by hand
# Instead edit example_hta_psoriasis.Rmd and then run precompile.R
skip_on_cran()
params <-
list(run_tests = FALSE)
## ----code=readLines("children/knitr_setup.R"), include=FALSE--------------------------------------
## ----eval = FALSE---------------------------------------------------------------------------------
## library(multinma)
## options(mc.cores = parallel::detectCores())
## ----setup, echo = FALSE--------------------------------------------------------------------------
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
"true" =, "warn" = 2,
parallel::detectCores())
options(mc.cores = nc)
set.seed(65498431)
## -------------------------------------------------------------------------------------------------
head(hta_psoriasis)
## -------------------------------------------------------------------------------------------------
dplyr::filter(hta_psoriasis, studyc %in% c("Elewski", "Gordon", "ACD2058g", "Altmeyer"))
## -------------------------------------------------------------------------------------------------
pso_net <- set_agd_arm(hta_psoriasis,
study = paste(studyc, year),
trt = trtc,
r = multi(r0 = sample_size - rowSums(cbind(PASI50, PASI75, PASI90), na.rm = TRUE),
PASI50, PASI75, PASI90,
inclusive = FALSE,
type = "ordered"))
pso_net
## ----hta_psoriasis_network_plot-------------------------------------------------------------------
plot(pso_net, weight_edges = TRUE, weight_nodes = TRUE) +
# Nudge the legend over
ggplot2::theme(legend.box.spacing = ggplot2::unit(0.75, "in"),
plot.margin = ggplot2::margin(0.1, 0, 0.1, 0.75, "in"))
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 10))
summary(normal(scale = 100))
## -------------------------------------------------------------------------------------------------
pso_fit_FE <- nma(pso_net,
trt_effects = "fixed",
link = "probit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_aux = flat())
## -------------------------------------------------------------------------------------------------
pso_fit_FE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(pso_fit_FE, pars = c("d", "mu", "cc"))
## ----pso_FE_pp_plot-------------------------------------------------------------------------------
plot_prior_posterior(pso_fit_FE)
## ----pso_FE_pp_cutpoint_plot----------------------------------------------------------------------
plot_prior_posterior(pso_fit_FE, prior = "aux")
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 10))
summary(normal(scale = 100))
summary(half_normal(scale = 2.5))
## ----eval = FALSE---------------------------------------------------------------------------------
## pso_fit_RE <- nma(pso_net,
## trt_effects = "random",
## link = "probit",
## prior_intercept = normal(scale = 100),
## prior_trt = normal(scale = 10),
## prior_aux = flat(),
## prior_het = half_normal(scale = 2.5),
## adapt_delta = 0.99)
## ----echo = FALSE, warning = FALSE----------------------------------------------------------------
pso_fit_RE <- nowarn_on_ci(nma(pso_net,
trt_effects = "random",
link = "probit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_aux = flat(),
prior_het = half_normal(scale = 2.5),
adapt_delta = 0.99,
iter = if (isTRUE(as.logical(Sys.getenv("CI")))) 5000 else 2000,
seed = 1713435794))
## -------------------------------------------------------------------------------------------------
pso_fit_RE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(pso_fit_RE, pars = c("d", "cc", "mu", "delta"))
## ----pso_RE_pp_plot-------------------------------------------------------------------------------
plot_prior_posterior(pso_fit_RE, prior = c("trt", "aux", "het"))
## -------------------------------------------------------------------------------------------------
(dic_FE <- dic(pso_fit_FE))
## -------------------------------------------------------------------------------------------------
(dic_RE <- dic(pso_fit_RE))
## ----pso_FE_resdev_plot---------------------------------------------------------------------------
plot(dic_FE)
## ----pso_RE_resdev_plot---------------------------------------------------------------------------
plot(dic_RE)
## ----pso_pred_FE, fig.height = 2------------------------------------------------------------------
pred_FE <- predict(pso_fit_FE,
baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5),
type = "response")
pred_FE
plot(pred_FE)
## ----pso_pred_RE, fig.height = 2------------------------------------------------------------------
pred_RE <- predict(pso_fit_RE,
baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5),
type = "response")
pred_RE
plot(pred_RE)
## ----pso_pred_FE_beta, fig.height = 2-------------------------------------------------------------
pred_FE_beta <- predict(pso_fit_FE,
baseline = distr(qbeta, 56, 408-56),
baseline_type = "response",
type = "response")
pred_FE_beta
plot(pred_FE_beta)
## ----pso_pred_RE_beta, fig.height = 2-------------------------------------------------------------
pred_RE_beta <- predict(pso_fit_RE,
baseline = distr(qbeta, 56, 408-56),
baseline_type = "response",
type = "response")
pred_RE_beta
plot(pred_RE_beta)
## ----pso_pred_RE_colour, fig.height = 3-----------------------------------------------------------
library(ggplot2)
plot(pred_RE, position = position_dodge(width = 0.75)) +
facet_null() +
aes(colour = Category) +
scale_colour_brewer(palette = "Blues")
## ----hta_psoriasis_ranks, fig.height=3------------------------------------------------------------
(pso_ranks <- posterior_ranks(pso_fit_RE, lower_better = FALSE))
plot(pso_ranks)
## ----hta_psoriasis_rankprobs----------------------------------------------------------------------
(pso_rankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE))
plot(pso_rankprobs)
## ----hta_psoriasis_cumrankprobs-------------------------------------------------------------------
(pso_cumrankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE, cumulative = TRUE))
plot(pso_cumrankprobs)
## ----hta_psoriasis_tests, include=FALSE, eval=params$run_tests------------------------------------
#--- Test against TSD 2 results ---
library(testthat)
library(dplyr)
tol <- 0.05
tol_dic <- 0.1
trt_codes <- c(
"1"= "Supportive care",
"2" = "Etanercept 25 mg",
"3"= "Etanercept 50 mg",
"4" = "Efalizumab",
"5"= "Ciclosporin",
"6" = "Fumaderm",
"7"= "Infliximab",
"8" = "Methotrexate")
# FE relative effects
tsd_fe <- tribble(
~contrast, ~est, ~sd, ~median, ~lower, ~upper,
"d12", -1.51, 0.10, -1.51, -1.70, -1.32,
"d13", -1.92, 0.10, -1.92, -2.12, -1.72,
"d14", -1.19, 0.06, -1.19, -1.30, -1.08,
"d15", -1.92, 0.34, -1.90, -2.62, -1.30,
"d16", -1.49, 0.49, -1.46, -2.55, -0.63,
"d17", -2.33, 0.26, -2.33, -2.87, -1.84,
"d18", -1.61, 0.44, -1.60, -2.50, -0.77) %>%
mutate(trt = recode(substr(contrast, 3, 3), !!! trt_codes),
trt = ordered(trt, levels = levels(pso_net$treatments)),
# Parameter estimates are opposite sign
.l = lower, .u = upper,
est = -est, median = -median,
lower = -.u, upper = -.l) %>%
arrange(trt)
pso_releff_FE <- as.data.frame(relative_effects(pso_fit_FE))
test_that("FE relative effects", {
expect_equivalent(pso_releff_FE$mean, tsd_fe$est, tolerance = tol)
expect_equivalent(pso_releff_FE$sd, tsd_fe$sd, tolerance = tol)
expect_equivalent(pso_releff_FE$`50%`, tsd_fe$median, tolerance = tol)
expect_equivalent(pso_releff_FE$`2.5%`, tsd_fe$lower, tolerance = tol)
expect_equivalent(pso_releff_FE$`97.5%`, tsd_fe$upper, tolerance = tol)
})
# FE predicted probabilities
tsd_pred_fe <- tribble(
~outcome, ~trt , ~mean, ~sd , ~median, ~lower, ~upper,
50 , "Supportive care" , 0.14 , 0.02, 0.14 , 0.10 , 0.18 ,
50 , "Etanercept 25 mg", 0.66 , 0.05, 0.66 , 0.56 , 0.75 ,
50 , "Etanercept 50 mg", 0.79 , 0.04, 0.79 , 0.71 , 0.86 ,
50 , "Efalizumab" , 0.54 , 0.04, 0.54 , 0.45 , 0.62 ,
50 , "Ciclosporin" , 0.78 , 0.10, 0.79 , 0.57 , 0.94 ,
50 , "Fumaderm" , 0.64 , 0.16, 0.64 , 0.31 , 0.93 ,
50 , "Infliximab" , 0.88 , 0.05, 0.89 , 0.76 , 0.96 ,
50 , "Methotrexate" , 0.68 , 0.15, 0.69 , 0.37 , 0.92 ,
75 , "Supportive care" , 0.03 , 0.01, 0.03 , 0.02 , 0.05 ,
75 , "Etanercept 25 mg", 0.37 , 0.05, 0.37 , 0.28 , 0.47 ,
75 , "Etanercept 50 mg", 0.53 , 0.05, 0.53 , 0.42 , 0.63 ,
75 , "Efalizumab" , 0.25 , 0.03, 0.25 , 0.19 , 0.33 ,
75 , "Ciclosporin" , 0.52 , 0.13, 0.52 , 0.28 , 0.79 ,
75 , "Fumaderm" , 0.37 , 0.17, 0.35 , 0.11 , 0.76 ,
75 , "Infliximab" , 0.68 , 0.10, 0.68 , 0.48 , 0.85 ,
75 , "Methotrexate" , 0.41 , 0.16, 0.40 , 0.14 , 0.75 ,
90 , "Supportive care" , 0.00 , 0.00, 0.00 , 0.00 , 0.01 ,
90 , "Etanercept 25 mg", 0.13 , 0.03, 0.13 , 0.08 , 0.19 ,
90 , "Etanercept 50 mg", 0.23 , 0.04, 0.23 , 0.16 , 0.32 ,
90 , "Efalizumab" , 0.07 , 0.02, 0.07 , 0.04 , 0.11 ,
90 , "Ciclosporin" , 0.24 , 0.11, 0.22 , 0.08 , 0.49 ,
90 , "Fumaderm" , 0.15 , 0.12, 0.11 , 0.02 , 0.46 ,
90 , "Infliximab" , 0.38 , 0.10, 0.37 , 0.19 , 0.60 ,
90 , "Methotrexate" , 0.17 , 0.11, 0.14 , 0.03 , 0.44 ) %>%
mutate(trt = ordered(trt, levels = levels(pso_net$treatments))) %>%
arrange(trt, outcome)
pred_FE <- as.data.frame(pred_FE)
test_that("FE predicted probabilities", {
expect_equivalent(pred_FE$mean, tsd_pred_fe$mean, tolerance = tol)
expect_equivalent(pred_FE$sd, tsd_pred_fe$sd, tolerance = tol)
expect_equivalent(pred_FE$`50%`, tsd_pred_fe$median, tolerance = tol)
expect_equivalent(pred_FE$`2.5%`, tsd_pred_fe$lower, tolerance = tol)
expect_equivalent(pred_FE$`97.5%`, tsd_pred_fe$upper, tolerance = tol)
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(dic_FE$resdev, 74.9, tolerance = tol_dic)
expect_equivalent(dic_FE$pd, 25.0, tolerance = tol_dic)
expect_equivalent(dic_FE$dic, 99.9, tolerance = tol_dic)
})
# RE relative effects
tsd_re <- tribble(
~contrast, ~est, ~sd, ~median, ~lower, ~upper,
"d12", -1.53, 0.24, -1.52, -2.05, -1.03,
"d13", -1.93, 0.28, -1.92, -2.51, -1.35,
"d14", -1.19, 0.18, -1.19, -1.56, -0.81,
"d15", -2.04, 0.43, -2.00, -3.02, -1.30,
"d16", -1.49, 0.62, -1.46, -2.81, -0.33,
"d17", -2.32, 0.38, -2.32, -3.06, -1.55,
"d18", -1.74, 0.64, -1.70, -3.14, -0.59) %>%
mutate(trt = recode(substr(contrast, 3, 3), !!! trt_codes),
trt = ordered(trt, levels = levels(pso_net$treatments)),
# Parameter estimates are opposite sign
.l = lower, .u = upper,
est = -est, median = -median,
lower = -.u, upper = -.l) %>%
arrange(trt)
pso_releff_RE <- as.data.frame(relative_effects(pso_fit_RE))
test_that("RE relative effects", {
expect_equivalent(pso_releff_RE$mean, tsd_re$est, tolerance = tol)
expect_equivalent(pso_releff_RE$sd, tsd_re$sd, tolerance = tol)
expect_equivalent(pso_releff_RE$`50%`, tsd_re$median, tolerance = tol)
expect_equivalent(pso_releff_RE$`2.5%`, tsd_re$lower, tolerance = tol)
expect_equivalent(pso_releff_RE$`97.5%`, tsd_re$upper, tolerance = tol)
})
# RE predicted probabilities
tsd_pred_re <- tribble(
~outcome, ~trt , ~mean, ~sd , ~median, ~lower, ~upper,
50 , "Supportive care" , 0.14 , 0.02, 0.14 , 0.10 , 0.18,
50 , "Etanercept 25 mg", 0.66 , 0.09, 0.66 , 0.46 , 0.83,
50 , "Etanercept 50 mg", 0.79 , 0.08, 0.80 , 0.59 , 0.92,
50 , "Efalizumab" , 0.54 , 0.08, 0.54 , 0.38 , 0.69,
50 , "Ciclosporin" , 0.81 , 0.10, 0.82 , 0.57 , 0.97,
50 , "Fumaderm" , 0.63 , 0.20, 0.64 , 0.22 , 0.96,
50 , "Infliximab" , 0.87 , 0.08, 0.89 , 0.67 , 0.98,
50 , "Methotrexate" , 0.70 , 0.18, 0.73 , 0.30 , 0.98,
75 , "Supportive care" , 0.03 , 0.01, 0.03 , 0.02 , 0.05,
75 , "Etanercept 25 mg", 0.38 , 0.09, 0.37 , 0.20 , 0.58,
75 , "Etanercept 50 mg", 0.53 , 0.11, 0.53 , 0.30 , 0.75,
75 , "Efalizumab" , 0.26 , 0.06, 0.25 , 0.14 , 0.40,
75 , "Ciclosporin" , 0.56 , 0.15, 0.56 , 0.28 , 0.88,
75 , "Fumaderm" , 0.38 , 0.20, 0.35 , 0.06 , 0.83,
75 , "Infliximab" , 0.67 , 0.13, 0.68 , 0.37 , 0.89,
75 , "Methotrexate" , 0.46 , 0.21, 0.44 , 0.10 , 0.90,
90 , "Supportive care" , 0.00 , 0.00, 0.00 , 0.00 , 0.01,
90 , "Etanercept 25 mg", 0.14 , 0.06, 0.13 , 0.05 , 0.28,
90 , "Etanercept 50 mg", 0.24 , 0.09, 0.23 , 0.09 , 0.45,
90 , "Efalizumab" , 0.07 , 0.03, 0.07 , 0.03 , 0.15,
90 , "Ciclosporin" , 0.28 , 0.14, 0.26 , 0.08 , 0.65,
90 , "Fumaderm" , 0.16 , 0.15, 0.12 , 0.01 , 0.57,
90 , "Infliximab" , 0.37 , 0.14, 0.37 , 0.13 , 0.67,
90 , "Methotrexate" , 0.21 , 0.17, 0.17 , 0.02 , 0.69) %>%
mutate(trt = ordered(trt, levels = levels(pso_net$treatments))) %>%
arrange(trt, outcome)
pred_RE <- as.data.frame(pred_RE)
test_that("RE predicted probabilities", {
expect_equivalent(pred_RE$mean, tsd_pred_re$mean, tolerance = tol)
expect_equivalent(pred_RE$sd, tsd_pred_re$sd, tolerance = tol)
expect_equivalent(pred_RE$`50%`, tsd_pred_re$median, tolerance = tol)
expect_equivalent(pred_RE$`2.5%`, tsd_pred_re$lower, tolerance = tol)
expect_equivalent(pred_RE$`97.5%`, tsd_pred_re$upper, tolerance = tol)
})
# Heterogeneity SD
pso_tau <- summary(pso_fit_RE, pars = "tau")
test_that("RE heterogeneity SD", {
expect_equivalent(pso_tau$summary$mean, 0.31, tolerance = tol)
expect_equivalent(pso_tau$summary$sd, 0.23, tolerance = tol)
expect_equivalent(pso_tau$summary$`50%`, 0.26, tolerance = tol)
skip_on_ci()
expect_equivalent(pso_tau$summary$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(pso_tau$summary$`97.5%`, 0.88, tolerance = tol)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(dic_RE$resdev, 63.0, tolerance = tol_dic)
expect_equivalent(dic_RE$pd, 33.3, tolerance = tol_dic)
expect_equivalent(dic_RE$dic, 96.2, tolerance = tol_dic)
})
# Check predictions with Beta distribution on baseline probability
pred_FE_beta <- as.data.frame(pred_FE_beta)
test_that("FE predicted probabilities (Beta distribution)", {
expect_equal(pred_FE[c("mean", "sd", "2.5%", "50%", "97.5%")],
pred_FE_beta[c("mean", "sd", "2.5%", "50%", "97.5%")],
tolerance = tol)
})
pred_RE_beta <- as.data.frame(pred_RE_beta)
test_that("RE predicted probabilities (Beta distribution)", {
expect_equal(pred_RE[c("mean", "sd", "2.5%", "50%", "97.5%")],
pred_RE_beta[c("mean", "sd", "2.5%", "50%", "97.5%")],
tolerance = tol)
})
# Force clean up
rm(list = ls())
gc()
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_hta_psoriasis.Rmd: do not edit by hand
# Instead edit example_hta_psoriasis.Rmd and then run precompile.R
skip_on_cran()
params <-
list(run_tests = FALSE)
## ----code=readLines("children/knitr_setup.R"), include=FALSE--------------------------------------
## ----eval = FALSE---------------------------------------------------------------------------------
## library(multinma)
## options(mc.cores = parallel::detectCores())
## ----setup, echo = FALSE--------------------------------------------------------------------------
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
"true" =, "warn" = 2,
parallel::detectCores())
options(mc.cores = nc)
set.seed(65498431)
## -------------------------------------------------------------------------------------------------
head(hta_psoriasis)
## -------------------------------------------------------------------------------------------------
dplyr::filter(hta_psoriasis, studyc %in% c("Elewski", "Gordon", "ACD2058g", "Altmeyer"))
## -------------------------------------------------------------------------------------------------
pso_net <- set_agd_arm(hta_psoriasis,
study = paste(studyc, year),
trt = trtc,
r = multi(r0 = sample_size - rowSums(cbind(PASI50, PASI75, PASI90), na.rm = TRUE),
PASI50, PASI75, PASI90,
inclusive = FALSE,
type = "ordered"))
pso_net
## ----hta_psoriasis_network_plot-------------------------------------------------------------------
plot(pso_net, weight_edges = TRUE, weight_nodes = TRUE) +
# Nudge the legend over
ggplot2::theme(legend.box.spacing = ggplot2::unit(0.75, "in"),
plot.margin = ggplot2::margin(0.1, 0, 0.1, 0.75, "in"))
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 10))
summary(normal(scale = 100))
## -------------------------------------------------------------------------------------------------
pso_fit_FE <- nma(pso_net,
trt_effects = "fixed",
link = "probit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_aux = flat())
## -------------------------------------------------------------------------------------------------
pso_fit_FE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(pso_fit_FE, pars = c("d", "mu", "cc"))
## ----pso_FE_pp_plot-------------------------------------------------------------------------------
plot_prior_posterior(pso_fit_FE)
## ----pso_FE_pp_cutpoint_plot----------------------------------------------------------------------
plot_prior_posterior(pso_fit_FE, prior = "aux")
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 10))
summary(normal(scale = 100))
summary(half_normal(scale = 2.5))
## ----eval = FALSE---------------------------------------------------------------------------------
## pso_fit_RE <- nma(pso_net,
## trt_effects = "random",
## link = "probit",
## prior_intercept = normal(scale = 100),
## prior_trt = normal(scale = 10),
## prior_aux = flat(),
## prior_het = half_normal(scale = 2.5),
## adapt_delta = 0.99)
## ----echo = FALSE, warning = FALSE----------------------------------------------------------------
pso_fit_RE <- nowarn_on_ci(nma(pso_net,
trt_effects = "random",
link = "probit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_aux = flat(),
prior_het = half_normal(scale = 2.5),
adapt_delta = 0.99,
iter = if (isTRUE(as.logical(Sys.getenv("CI")))) 5000 else 2000,
seed = 1713435794))
## -------------------------------------------------------------------------------------------------
pso_fit_RE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(pso_fit_RE, pars = c("d", "cc", "mu", "delta"))
## ----pso_RE_pp_plot-------------------------------------------------------------------------------
plot_prior_posterior(pso_fit_RE, prior = c("trt", "aux", "het"))
## -------------------------------------------------------------------------------------------------
(dic_FE <- dic(pso_fit_FE))
## -------------------------------------------------------------------------------------------------
(dic_RE <- dic(pso_fit_RE))
## ----pso_FE_resdev_plot---------------------------------------------------------------------------
plot(dic_FE)
## ----pso_RE_resdev_plot---------------------------------------------------------------------------
plot(dic_RE)
## ----pso_pred_FE, fig.height = 2------------------------------------------------------------------
pred_FE <- predict(pso_fit_FE,
baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5),
type = "response")
pred_FE
plot(pred_FE)
## ----pso_pred_RE, fig.height = 2------------------------------------------------------------------
pred_RE <- predict(pso_fit_RE,
baseline = distr(qnorm, mean = -1.097, sd = 123^-0.5),
type = "response")
pred_RE
plot(pred_RE)
## ----pso_pred_FE_beta, fig.height = 2-------------------------------------------------------------
pred_FE_beta <- predict(pso_fit_FE,
baseline = distr(qbeta, 56, 408-56),
baseline_type = "response",
type = "response")
pred_FE_beta
plot(pred_FE_beta)
## ----pso_pred_RE_beta, fig.height = 2-------------------------------------------------------------
pred_RE_beta <- predict(pso_fit_RE,
baseline = distr(qbeta, 56, 408-56),
baseline_type = "response",
type = "response")
pred_RE_beta
plot(pred_RE_beta)
## ----pso_pred_RE_colour, fig.height = 3-----------------------------------------------------------
library(ggplot2)
plot(pred_RE, position = position_dodge(width = 0.75)) +
facet_null() +
aes(colour = Category) +
scale_colour_brewer(palette = "Blues")
## ----hta_psoriasis_ranks, fig.height=3------------------------------------------------------------
(pso_ranks <- posterior_ranks(pso_fit_RE, lower_better = FALSE))
plot(pso_ranks)
## ----hta_psoriasis_rankprobs----------------------------------------------------------------------
(pso_rankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE))
plot(pso_rankprobs)
## ----hta_psoriasis_cumrankprobs-------------------------------------------------------------------
(pso_cumrankprobs <- posterior_rank_probs(pso_fit_RE, lower_better = FALSE, cumulative = TRUE))
plot(pso_cumrankprobs)
## ----hta_psoriasis_tests, include=FALSE, eval=params$run_tests------------------------------------
#--- Test against TSD 2 results ---
library(testthat)
library(dplyr)
tol <- 0.05
tol_dic <- 0.1
trt_codes <- c(
"1"= "Supportive care",
"2" = "Etanercept 25 mg",
"3"= "Etanercept 50 mg",
"4" = "Efalizumab",
"5"= "Ciclosporin",
"6" = "Fumaderm",
"7"= "Infliximab",
"8" = "Methotrexate")
# FE relative effects
tsd_fe <- tribble(
~contrast, ~est, ~sd, ~median, ~lower, ~upper,
"d12", -1.51, 0.10, -1.51, -1.70, -1.32,
"d13", -1.92, 0.10, -1.92, -2.12, -1.72,
"d14", -1.19, 0.06, -1.19, -1.30, -1.08,
"d15", -1.92, 0.34, -1.90, -2.62, -1.30,
"d16", -1.49, 0.49, -1.46, -2.55, -0.63,
"d17", -2.33, 0.26, -2.33, -2.87, -1.84,
"d18", -1.61, 0.44, -1.60, -2.50, -0.77) %>%
mutate(trt = recode(substr(contrast, 3, 3), !!! trt_codes),
trt = ordered(trt, levels = levels(pso_net$treatments)),
# Parameter estimates are opposite sign
.l = lower, .u = upper,
est = -est, median = -median,
lower = -.u, upper = -.l) %>%
arrange(trt)
pso_releff_FE <- as.data.frame(relative_effects(pso_fit_FE))
test_that("FE relative effects", {
expect_equivalent(pso_releff_FE$mean, tsd_fe$est, tolerance = tol)
expect_equivalent(pso_releff_FE$sd, tsd_fe$sd, tolerance = tol)
expect_equivalent(pso_releff_FE$`50%`, tsd_fe$median, tolerance = tol)
expect_equivalent(pso_releff_FE$`2.5%`, tsd_fe$lower, tolerance = tol)
expect_equivalent(pso_releff_FE$`97.5%`, tsd_fe$upper, tolerance = tol)
})
# FE predicted probabilities
tsd_pred_fe <- tribble(
~outcome, ~trt , ~mean, ~sd , ~median, ~lower, ~upper,
50 , "Supportive care" , 0.14 , 0.02, 0.14 , 0.10 , 0.18 ,
50 , "Etanercept 25 mg", 0.66 , 0.05, 0.66 , 0.56 , 0.75 ,
50 , "Etanercept 50 mg", 0.79 , 0.04, 0.79 , 0.71 , 0.86 ,
50 , "Efalizumab" , 0.54 , 0.04, 0.54 , 0.45 , 0.62 ,
50 , "Ciclosporin" , 0.78 , 0.10, 0.79 , 0.57 , 0.94 ,
50 , "Fumaderm" , 0.64 , 0.16, 0.64 , 0.31 , 0.93 ,
50 , "Infliximab" , 0.88 , 0.05, 0.89 , 0.76 , 0.96 ,
50 , "Methotrexate" , 0.68 , 0.15, 0.69 , 0.37 , 0.92 ,
75 , "Supportive care" , 0.03 , 0.01, 0.03 , 0.02 , 0.05 ,
75 , "Etanercept 25 mg", 0.37 , 0.05, 0.37 , 0.28 , 0.47 ,
75 , "Etanercept 50 mg", 0.53 , 0.05, 0.53 , 0.42 , 0.63 ,
75 , "Efalizumab" , 0.25 , 0.03, 0.25 , 0.19 , 0.33 ,
75 , "Ciclosporin" , 0.52 , 0.13, 0.52 , 0.28 , 0.79 ,
75 , "Fumaderm" , 0.37 , 0.17, 0.35 , 0.11 , 0.76 ,
75 , "Infliximab" , 0.68 , 0.10, 0.68 , 0.48 , 0.85 ,
75 , "Methotrexate" , 0.41 , 0.16, 0.40 , 0.14 , 0.75 ,
90 , "Supportive care" , 0.00 , 0.00, 0.00 , 0.00 , 0.01 ,
90 , "Etanercept 25 mg", 0.13 , 0.03, 0.13 , 0.08 , 0.19 ,
90 , "Etanercept 50 mg", 0.23 , 0.04, 0.23 , 0.16 , 0.32 ,
90 , "Efalizumab" , 0.07 , 0.02, 0.07 , 0.04 , 0.11 ,
90 , "Ciclosporin" , 0.24 , 0.11, 0.22 , 0.08 , 0.49 ,
90 , "Fumaderm" , 0.15 , 0.12, 0.11 , 0.02 , 0.46 ,
90 , "Infliximab" , 0.38 , 0.10, 0.37 , 0.19 , 0.60 ,
90 , "Methotrexate" , 0.17 , 0.11, 0.14 , 0.03 , 0.44 ) %>%
mutate(trt = ordered(trt, levels = levels(pso_net$treatments))) %>%
arrange(trt, outcome)
pred_FE <- as.data.frame(pred_FE)
test_that("FE predicted probabilities", {
expect_equivalent(pred_FE$mean, tsd_pred_fe$mean, tolerance = tol)
expect_equivalent(pred_FE$sd, tsd_pred_fe$sd, tolerance = tol)
expect_equivalent(pred_FE$`50%`, tsd_pred_fe$median, tolerance = tol)
expect_equivalent(pred_FE$`2.5%`, tsd_pred_fe$lower, tolerance = tol)
expect_equivalent(pred_FE$`97.5%`, tsd_pred_fe$upper, tolerance = tol)
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(dic_FE$resdev, 74.9, tolerance = tol_dic)
expect_equivalent(dic_FE$pd, 25.0, tolerance = tol_dic)
expect_equivalent(dic_FE$dic, 99.9, tolerance = tol_dic)
})
# RE relative effects
tsd_re <- tribble(
~contrast, ~est, ~sd, ~median, ~lower, ~upper,
"d12", -1.53, 0.24, -1.52, -2.05, -1.03,
"d13", -1.93, 0.28, -1.92, -2.51, -1.35,
"d14", -1.19, 0.18, -1.19, -1.56, -0.81,
"d15", -2.04, 0.43, -2.00, -3.02, -1.30,
"d16", -1.49, 0.62, -1.46, -2.81, -0.33,
"d17", -2.32, 0.38, -2.32, -3.06, -1.55,
"d18", -1.74, 0.64, -1.70, -3.14, -0.59) %>%
mutate(trt = recode(substr(contrast, 3, 3), !!! trt_codes),
trt = ordered(trt, levels = levels(pso_net$treatments)),
# Parameter estimates are opposite sign
.l = lower, .u = upper,
est = -est, median = -median,
lower = -.u, upper = -.l) %>%
arrange(trt)
pso_releff_RE <- as.data.frame(relative_effects(pso_fit_RE))
test_that("RE relative effects", {
expect_equivalent(pso_releff_RE$mean, tsd_re$est, tolerance = tol)
expect_equivalent(pso_releff_RE$sd, tsd_re$sd, tolerance = tol)
expect_equivalent(pso_releff_RE$`50%`, tsd_re$median, tolerance = tol)
expect_equivalent(pso_releff_RE$`2.5%`, tsd_re$lower, tolerance = tol)
expect_equivalent(pso_releff_RE$`97.5%`, tsd_re$upper, tolerance = tol)
})
# RE predicted probabilities
tsd_pred_re <- tribble(
~outcome, ~trt , ~mean, ~sd , ~median, ~lower, ~upper,
50 , "Supportive care" , 0.14 , 0.02, 0.14 , 0.10 , 0.18,
50 , "Etanercept 25 mg", 0.66 , 0.09, 0.66 , 0.46 , 0.83,
50 , "Etanercept 50 mg", 0.79 , 0.08, 0.80 , 0.59 , 0.92,
50 , "Efalizumab" , 0.54 , 0.08, 0.54 , 0.38 , 0.69,
50 , "Ciclosporin" , 0.81 , 0.10, 0.82 , 0.57 , 0.97,
50 , "Fumaderm" , 0.63 , 0.20, 0.64 , 0.22 , 0.96,
50 , "Infliximab" , 0.87 , 0.08, 0.89 , 0.67 , 0.98,
50 , "Methotrexate" , 0.70 , 0.18, 0.73 , 0.30 , 0.98,
75 , "Supportive care" , 0.03 , 0.01, 0.03 , 0.02 , 0.05,
75 , "Etanercept 25 mg", 0.38 , 0.09, 0.37 , 0.20 , 0.58,
75 , "Etanercept 50 mg", 0.53 , 0.11, 0.53 , 0.30 , 0.75,
75 , "Efalizumab" , 0.26 , 0.06, 0.25 , 0.14 , 0.40,
75 , "Ciclosporin" , 0.56 , 0.15, 0.56 , 0.28 , 0.88,
75 , "Fumaderm" , 0.38 , 0.20, 0.35 , 0.06 , 0.83,
75 , "Infliximab" , 0.67 , 0.13, 0.68 , 0.37 , 0.89,
75 , "Methotrexate" , 0.46 , 0.21, 0.44 , 0.10 , 0.90,
90 , "Supportive care" , 0.00 , 0.00, 0.00 , 0.00 , 0.01,
90 , "Etanercept 25 mg", 0.14 , 0.06, 0.13 , 0.05 , 0.28,
90 , "Etanercept 50 mg", 0.24 , 0.09, 0.23 , 0.09 , 0.45,
90 , "Efalizumab" , 0.07 , 0.03, 0.07 , 0.03 , 0.15,
90 , "Ciclosporin" , 0.28 , 0.14, 0.26 , 0.08 , 0.65,
90 , "Fumaderm" , 0.16 , 0.15, 0.12 , 0.01 , 0.57,
90 , "Infliximab" , 0.37 , 0.14, 0.37 , 0.13 , 0.67,
90 , "Methotrexate" , 0.21 , 0.17, 0.17 , 0.02 , 0.69) %>%
mutate(trt = ordered(trt, levels = levels(pso_net$treatments))) %>%
arrange(trt, outcome)
pred_RE <- as.data.frame(pred_RE)
test_that("RE predicted probabilities", {
expect_equivalent(pred_RE$mean, tsd_pred_re$mean, tolerance = tol)
expect_equivalent(pred_RE$sd, tsd_pred_re$sd, tolerance = tol)
expect_equivalent(pred_RE$`50%`, tsd_pred_re$median, tolerance = tol)
expect_equivalent(pred_RE$`2.5%`, tsd_pred_re$lower, tolerance = tol)
expect_equivalent(pred_RE$`97.5%`, tsd_pred_re$upper, tolerance = tol)
})
# Heterogeneity SD
pso_tau <- summary(pso_fit_RE, pars = "tau")
test_that("RE heterogeneity SD", {
expect_equivalent(pso_tau$summary$mean, 0.31, tolerance = tol)
expect_equivalent(pso_tau$summary$sd, 0.23, tolerance = tol)
expect_equivalent(pso_tau$summary$`50%`, 0.26, tolerance = tol)
skip_on_ci()
expect_equivalent(pso_tau$summary$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(pso_tau$summary$`97.5%`, 0.88, tolerance = tol)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(dic_RE$resdev, 63.0, tolerance = tol_dic)
expect_equivalent(dic_RE$pd, 33.3, tolerance = tol_dic)
expect_equivalent(dic_RE$dic, 96.2, tolerance = tol_dic)
})
# Check predictions with Beta distribution on baseline probability
pred_FE_beta <- as.data.frame(pred_FE_beta)
test_that("FE predicted probabilities (Beta distribution)", {
expect_equal(pred_FE[c("mean", "sd", "2.5%", "50%", "97.5%")],
pred_FE_beta[c("mean", "sd", "2.5%", "50%", "97.5%")],
tolerance = tol)
})
pred_RE_beta <- as.data.frame(pred_RE_beta)
test_that("RE predicted probabilities (Beta distribution)", {
expect_equal(pred_RE[c("mean", "sd", "2.5%", "50%", "97.5%")],
pred_RE_beta[c("mean", "sd", "2.5%", "50%", "97.5%")],
tolerance = tol)
})
# Force clean up
rm(list = ls())
gc()
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.