tests/testthat/test-example_thrombolytics.R
In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data
# Generated by vignette example_thrombolytics.Rmd: do not edit by hand
# Instead edit example_thrombolytics.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)
## -------------------------------------------------------------------------------------------------
head(thrombolytics)
## -------------------------------------------------------------------------------------------------
thrombo_net <- set_agd_arm(thrombolytics,
study = studyn,
trt = trtc,
r = r,
n = n)
thrombo_net
## ----include=FALSE, eval=params$run_tests---------------------------------------------------------
# Make trtf factor to order treatments in same way as Dias analysis - needed to
# recreate inconsistency analyses
trts <- dplyr::distinct(thrombolytics, trtn, trtc)
trts <- dplyr::arrange(trts, trtn)
thrombolytics$trtf <- factor(thrombolytics$trtn, levels = trts$trtn, labels = trts$trtc)
thrombo_net2 <- set_agd_arm(thrombolytics,
study = studyn,
trt = trtf,
r = r,
n = n)
## ----eval=FALSE-----------------------------------------------------------------------------------
## plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE)
## ----thrombo_net_plot, echo=FALSE-----------------------------------------------------------------
plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(legend.margin = ggplot2::margin(l = 4, unit = "lines"))
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
## -------------------------------------------------------------------------------------------------
thrombo_fit <- nma(thrombo_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
## -------------------------------------------------------------------------------------------------
thrombo_fit
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(thrombo_fit, pars = c("d", "mu"))
## ----thrombo_pp_plot------------------------------------------------------------------------------
plot_prior_posterior(thrombo_fit, prior = "trt")
## -------------------------------------------------------------------------------------------------
(dic_consistency <- dic(thrombo_fit))
## ----thrombo_resdev_plot, fig.width=12------------------------------------------------------------
plot(dic_consistency)
## -------------------------------------------------------------------------------------------------
thrombo_fit_ume <- nma(thrombo_net,
consistency = "ume",
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
thrombo_fit_ume
## -------------------------------------------------------------------------------------------------
dic_consistency
(dic_ume <- dic(thrombo_fit_ume))
## ----thrombo_devdev_plot--------------------------------------------------------------------------
plot(dic_consistency, dic_ume, show_uncertainty = FALSE)
## ----eval=!params$run_tests-----------------------------------------------------------------------
## thrombo_nodesplit <- nma(thrombo_net,
## consistency = "nodesplit",
## trt_effects = "fixed",
## prior_intercept = normal(scale = 100),
## prior_trt = normal(scale = 100))
## ----include=FALSE, eval=params$run_tests---------------------------------------------------------
# Run node-splits with treatments ordered as per Dias
thrombo_nodesplit <- nma(thrombo_net2,
consistency = "nodesplit",
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
## -------------------------------------------------------------------------------------------------
summary(thrombo_nodesplit)
## ----thrombo_nodesplit, fig.width = 7-------------------------------------------------------------
plot(thrombo_nodesplit)
## ----thrombo_nodesplit_omega, fig.height = 6------------------------------------------------------
plot(thrombo_nodesplit, pars = "omega", stat = "halfeye", ref_line = 0) +
ggplot2::aes(y = comparison) +
ggplot2::facet_null()
## ----thrombo_releff-------------------------------------------------------------------------------
(thrombo_releff <- relative_effects(thrombo_fit, all_contrasts = TRUE))
plot(thrombo_releff, ref_line = 0)
## ----thrombo_ranks--------------------------------------------------------------------------------
(thrombo_ranks <- posterior_ranks(thrombo_fit))
plot(thrombo_ranks)
## ----thrombo_rankprobs----------------------------------------------------------------------------
(thrombo_rankprobs <- posterior_rank_probs(thrombo_fit))
plot(thrombo_rankprobs)
## ----thrombo_cumrankprobs-------------------------------------------------------------------------
(thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE))
plot(thrombo_cumrankprobs)
## ----thrombo_tests, include=FALSE, eval=params$run_tests------------------------------------------
#--- Test against TSD 4 results ---
library(testthat)
library(dplyr)
library(tidyr)
test_that("Reference trt is SK", {
expect_equivalent(levels(thrombo_net$treatments)[1], "SK")
})
tol <- 0.05
tol_dic <- 0.1
# Relative effects
tsd_releff <- tribble(
~trt_b , ~trt , ~mean , ~sd , ~lower, ~upper,
"SK" , "t-PA" , 0.002 , 0.030, -0.06 , 0.06 ,
"SK" , "Acc t-PA" , -0.177, 0.043, -0.26 , -0.09 ,
"SK" , "SK + t-PA", -0.049, 0.046, -0.14 , 0.04 ,
"SK" , "r-PA" , -0.124, 0.060, -0.24 , -0.01 ,
"SK" , "PTCA" , -0.476, 0.101, -0.67 , -0.28 ,
"SK" , "UK" , -0.203, 0.221, -0.64 , 0.23 ,
"SK" , "ASPAC" , 0.016 , 0.037, -0.06 , 0.09 ,
"t-PA" , "PTCA" , -0.478, 0.104, -0.68 , -0.27 ,
"t-PA" , "UK" , -0.206, 0.221, -0.64 , 0.23 ,
"t-PA" , "ASPAC" , 0.013 , 0.037, -0.06 , 0.09 ,
"Acc t-PA", "r-PA" , 0.054 , 0.055, -0.05 , 0.16 ,
"Acc t-PA", "TNK" , 0.005 , 0.064, -0.12 , 0.13 ,
"Acc t-PA", "PTCA" , -0.298, 0.098, -0.49 , -0.11 ,
"Acc t-PA", "UK" , -0.026, 0.221, -0.46 , 0.41 ,
"Acc t-PA", "ASPAC" , 0.193 , 0.056, 0.08 , 0.30 ) %>%
mutate(.trt_b = ordered(trt_b, levels = levels(thrombo_net$treatments)),
.trt = ordered(trt, levels = levels(thrombo_net$treatments)),
rev = if_else(.trt_b > .trt, -1, 1),
.l = lower, .u = upper,
lower = if_else(.trt_b > .trt, .u, .l),
upper = if_else(.trt_b > .trt, .l, .u),
trt_b = if_else(.trt_b > .trt, .trt, .trt_b),
trt = if_else(.trt_b > .trt, .trt_b, .trt),
lab = paste0("d[", trt, " vs. ", trt_b, "]")) %>%
arrange(trt_b, trt) %>%
mutate_at(vars(mean, lower, upper), ~.*rev)
thrombo_releff_summary <- as.data.frame(thrombo_releff) %>%
filter(parameter %in% tsd_releff$lab)
test_that("FE relative effects", {
expect_equivalent(thrombo_releff_summary$mean, tsd_releff$mean, tolerance = tol)
expect_equivalent(thrombo_releff_summary$sd, tsd_releff$sd, tolerance = tol)
expect_equivalent(thrombo_releff_summary$`2.5%`, tsd_releff$lower, tolerance = tol)
expect_equivalent(thrombo_releff_summary$`97.5%`, tsd_releff$upper, tolerance = tol)
})
test_that("SUCRAs", {
thrombo_ranks <- posterior_ranks(thrombo_fit, sucra = TRUE)
thrombo_rankprobs <- posterior_rank_probs(thrombo_fit, sucra = TRUE)
thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE, sucra = TRUE)
expect_equal(thrombo_ranks$summary$sucra, thrombo_rankprobs$summary$sucra)
expect_equal(thrombo_ranks$summary$sucra, thrombo_cumrankprobs$summary$sucra)
})
# DIC
test_that("DIC", {
expect_equivalent(dic_consistency$resdev, 105.9, tolerance = tol_dic)
expect_equivalent(dic_consistency$pd, 58, tolerance = tol_dic)
expect_equivalent(dic_consistency$dic, 163.9, tolerance = tol_dic)
})
# Relative effects (UME)
# FE UME model, so no differences by reference treatment, no multi-arm correction
tsd_ume <- tribble(
~trt_b , ~trt , ~mean , ~sd , ~lower, ~upper,
"SK" , "t-PA" , -0.004, 0.030, -0.06 , 0.06 ,
"SK" , "Acc t-PA" , -0.158, 0.049, -0.25 , -0.06 ,
"SK" , "SK + t-PA", -0.044, 0.047, -0.14 , 0.05 ,
"SK" , "r-PA" , -0.060, 0.089, -0.23 , 0.11 ,
"SK" , "PTCA" , -0.665, 0.185, -1.03 , -0.31 ,
"SK" , "UK" , -0.369, 0.518, -1.41 , 0.63 ,
"SK" , "ASPAC" , 0.005 , 0.037, -0.07 , 0.08 ,
"t-PA" , "PTCA" , -0.544, 0.417, -1.38 , 0.25 ,
"t-PA" , "UK" , -0.294, 0.347, -0.99 , 0.37 ,
"t-PA" , "ASPAC" , -0.290, 0.361, -1.01 , 0.41 ,
"Acc t-PA", "r-PA" , 0.019 , 0.066, -0.11 , 0.15 ,
"Acc t-PA", "TNK" , 0.006 , 0.064, -0.12 , 0.13 ,
"Acc t-PA", "PTCA" , -0.216, 0.119, -0.45 , 0.02 ,
"Acc t-PA", "UK" , 0.146 , 0.358, -0.54 , 0.86 ,
"Acc t-PA", "ASPAC" , 1.405 , 0.417, 0.63 , 2.27 ) %>%
mutate(.trt_b = ordered(trt_b, levels = levels(thrombo_net$treatments)),
.trt = ordered(trt, levels = levels(thrombo_net$treatments)),
rev = if_else(.trt_b > .trt, -1, 1),
.l = lower, .u = upper,
lower = if_else(.trt_b > .trt, .u, .l),
upper = if_else(.trt_b > .trt, .l, .u),
trt_b = if_else(.trt_b > .trt, .trt, .trt_b),
trt = if_else(.trt_b > .trt, .trt_b, .trt),
lab = paste0("d[", trt, " vs. ", trt_b, "]")) %>%
arrange(trt_b, trt) %>%
mutate_at(vars(mean, lower, upper), ~.*rev)
thrombo_ume_releff <- summary(thrombo_fit_ume, pars = "d")
test_that("UME relative effects", {
expect_equivalent(thrombo_ume_releff$summary$mean, tsd_ume$mean, tolerance = tol)
expect_equivalent(thrombo_ume_releff$summary$sd, tsd_ume$sd, tolerance = tol)
expect_equivalent(thrombo_ume_releff$summary$`2.5%`, tsd_ume$lower, tolerance = tol)
expect_equivalent(thrombo_ume_releff$summary$`97.5%`, tsd_ume$upper, tolerance = tol)
})
# DIC (UME)
test_that("UME DIC", {
expect_equivalent(dic_ume$resdev, 99.7, tolerance = tol_dic)
expect_equivalent(dic_ume$pd, 65, tolerance = tol_dic)
expect_equivalent(dic_ume$dic, 164.7, tolerance = tol_dic)
})
# Node-splitting
trt_code <- c("SK", "t-PA", "Acc t-PA", "SK + t-PA", "r-PA",
"TNK", "PTCA", "UK", "ASPAC")
dias2010_nodesplit_est <- tribble(
~trt1, ~trt2, ~net_mean, ~net_sd, ~dir_mean, ~dir_sd, ~ind_mean, ~ind_sd, ~omega_mean, ~omega_sd, ~p_value,
1, 2, 0.002, 0.030, 0.000, 0.030, 0.189, 0.235, -0.190, 0.236, 0.42,
1, 3, -0.177, 0.043, -0.158, 0.048, -0.247, 0.092, 0.088, 0.104, 0.39,
1, 5, -0.124, 0.060, -0.060, 0.089, -0.175, 0.081, 0.115, 0.121, 0.34,
1, 7, -0.475, 0.101, -0.666, 0.185, -0.393, 0.120, -0.272, 0.222, 0.22,
1, 8, -0.203, 0.219, -0.369, 0.518, -0.168, 0.244, -0.207, 0.575, 0.73,
1, 9, 0.016, 0.037, 0.009, 0.037, 0.424, 0.252, -0.413, 0.253, 0.10,
2, 7, -0.477, 0.104, -0.545, 0.417, -0.475, 0.108, -0.073, 0.432, 0.88,
2, 8, -0.205, 0.220, -0.295, 0.347, -0.144, 0.290, -0.155, 0.452, 0.74,
2, 9, 0.014, 0.037, 0.006, 0.037, 0.471, 0.241, -0.468, 0.241, 0.05,
3, 4, 0.128, 0.054, 0.126, 0.054, 0.630, 0.697, -0.506, 0.696, 0.47,
3, 5, 0.053, 0.056, 0.019, 0.066, 0.135, 0.101, -0.116, 0.120, 0.34,
3, 7, -0.298, 0.097, -0.216, 0.118, -0.477, 0.174, 0.260, 0.211, 0.21,
3, 8, -0.026, 0.220, 0.143, 0.356, -0.136, 0.288, 0.277, 0.461, 0.55,
3, 9, 0.193, 0.056, 1.409, 0.415, 0.165, 0.057, 1.239, 0.420, 0.001
) %>%
mutate(trt1 = forcats::fct_relevel(factor(trt1, levels = 1:9, labels = trt_code),
!! levels(thrombo_net2$treatments)),
trt2 = forcats::fct_relevel(factor(trt2, levels = 1:9, labels = trt_code),
!! levels(thrombo_net2$treatments)))
for (i in 1:nrow(dias2010_nodesplit_est)) {
if (as.numeric(dias2010_nodesplit_est$trt1[i]) > as.numeric(dias2010_nodesplit_est$trt2[i])) {
dias2010_nodesplit_est[i, c("trt1", "trt2")] <- dias2010_nodesplit_est[i, c("trt2", "trt1")]
dias2010_nodesplit_est[i, c("net_mean", "dir_mean", "ind_mean")] <- -dias2010_nodesplit_est[i, c("net_mean", "dir_mean", "ind_mean")]
}
}
dias2010_nodesplit_est <- arrange(dias2010_nodesplit_est, trt1, trt2)
thrombo_nodesplit_est <- as_tibble(summary(thrombo_nodesplit), nest = FALSE) %>%
mutate(parameter = stringr::str_replace(parameter, "(^d_)?(.+)\\[.+\\]$", "\\2")) %>%
pivot_wider(names_from = "parameter",
values_from = c("mean", "sd"),
names_glue = "{parameter}_{.value}",
id_cols = c("trt1", "trt2", "p_value")) %>%
mutate(across(where(is.numeric), unname)) %>%
select(!!! colnames(dias2010_nodesplit_est))
test_that("Node-splitting estimates", {
# expect_equal(thrombo_nodesplit_est, dias2010_nodesplit_est, tolerance = tol)
expect_equal(thrombo_nodesplit_est$net_mean, dias2010_nodesplit_est$net_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$net_sd, dias2010_nodesplit_est$net_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$ind_mean, dias2010_nodesplit_est$ind_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$ind_sd, dias2010_nodesplit_est$ind_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$dir_mean, dias2010_nodesplit_est$dir_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$dir_sd, dias2010_nodesplit_est$dir_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$omega_mean, dias2010_nodesplit_est$omega_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$omega_sd, dias2010_nodesplit_est$omega_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$p_value, dias2010_nodesplit_est$p_value, tolerance = tol)
})
dias2010_nodesplit_dic <- tribble(
~trt1, ~trt2, ~resdev, ~pd, ~dic,
1, 2, 106.4, 58.7, 165.1,
1, 3, 106.2, 58.7, 165.0,
1, 5, 106.1, 58.7, 164.8,
1, 7, 105.5, 58.7, 164.2,
1, 8, 106.9, 58.7, 165.6,
1, 9, 104.3, 58.7, 163.0,
2, 7, 106.9, 58.7, 165.6,
2, 8, 106.8, 58.7, 165.5,
2, 9, 103.3, 58.7, 162.0,
3, 4, 106.5, 58.7, 165.2,
3, 5, 106.1, 58.7, 164.8,
3, 7, 105.5, 58.7, 164.2,
3, 8, 106.6, 58.7, 165.3,
3, 9, 96.9, 58.7, 155.6) %>%
mutate(trt1 = forcats::fct_relevel(factor(trt1, levels = 1:9, labels = trt_code),
!! levels(thrombo_net$treatments)),
trt2 = forcats::fct_relevel(factor(trt2, levels = 1:9, labels = trt_code),
!! levels(thrombo_net$treatments)))
for (i in 1:nrow(dias2010_nodesplit_dic)) {
if (as.numeric(dias2010_nodesplit_dic$trt1[i]) > as.numeric(dias2010_nodesplit_dic$trt2[i])) {
dias2010_nodesplit_dic[i, c("trt1", "trt2")] <- dias2010_nodesplit_dic[i, c("trt2", "trt1")]
}
}
dias2010_nodesplit_dic <- arrange(dias2010_nodesplit_dic, trt1, trt2)
thrombo_nodesplit_dic <- as_tibble(summary(thrombo_nodesplit), nest = FALSE) %>%
distinct(trt1, trt2, resdev, pd, dic)
test_that("Node-splitting DIC", {
expect_equal(thrombo_nodesplit_dic$dic, dias2010_nodesplit_dic$dic, tolerance = tol_dic)
expect_equal(thrombo_nodesplit_dic$resdev, dias2010_nodesplit_dic$resdev, tolerance = tol_dic)
expect_equal(thrombo_nodesplit_dic$pd, dias2010_nodesplit_dic$pd, tolerance = tol_dic)
})
# 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_thrombolytics.Rmd: do not edit by hand
# Instead edit example_thrombolytics.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)
## -------------------------------------------------------------------------------------------------
head(thrombolytics)
## -------------------------------------------------------------------------------------------------
thrombo_net <- set_agd_arm(thrombolytics,
study = studyn,
trt = trtc,
r = r,
n = n)
thrombo_net
## ----include=FALSE, eval=params$run_tests---------------------------------------------------------
# Make trtf factor to order treatments in same way as Dias analysis - needed to
# recreate inconsistency analyses
trts <- dplyr::distinct(thrombolytics, trtn, trtc)
trts <- dplyr::arrange(trts, trtn)
thrombolytics$trtf <- factor(thrombolytics$trtn, levels = trts$trtn, labels = trts$trtc)
thrombo_net2 <- set_agd_arm(thrombolytics,
study = studyn,
trt = trtf,
r = r,
n = n)
## ----eval=FALSE-----------------------------------------------------------------------------------
## plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE)
## ----thrombo_net_plot, echo=FALSE-----------------------------------------------------------------
plot(thrombo_net, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(legend.margin = ggplot2::margin(l = 4, unit = "lines"))
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
## -------------------------------------------------------------------------------------------------
thrombo_fit <- nma(thrombo_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
## -------------------------------------------------------------------------------------------------
thrombo_fit
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(thrombo_fit, pars = c("d", "mu"))
## ----thrombo_pp_plot------------------------------------------------------------------------------
plot_prior_posterior(thrombo_fit, prior = "trt")
## -------------------------------------------------------------------------------------------------
(dic_consistency <- dic(thrombo_fit))
## ----thrombo_resdev_plot, fig.width=12------------------------------------------------------------
plot(dic_consistency)
## -------------------------------------------------------------------------------------------------
thrombo_fit_ume <- nma(thrombo_net,
consistency = "ume",
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
thrombo_fit_ume
## -------------------------------------------------------------------------------------------------
dic_consistency
(dic_ume <- dic(thrombo_fit_ume))
## ----thrombo_devdev_plot--------------------------------------------------------------------------
plot(dic_consistency, dic_ume, show_uncertainty = FALSE)
## ----eval=!params$run_tests-----------------------------------------------------------------------
## thrombo_nodesplit <- nma(thrombo_net,
## consistency = "nodesplit",
## trt_effects = "fixed",
## prior_intercept = normal(scale = 100),
## prior_trt = normal(scale = 100))
## ----include=FALSE, eval=params$run_tests---------------------------------------------------------
# Run node-splits with treatments ordered as per Dias
thrombo_nodesplit <- nma(thrombo_net2,
consistency = "nodesplit",
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
## -------------------------------------------------------------------------------------------------
summary(thrombo_nodesplit)
## ----thrombo_nodesplit, fig.width = 7-------------------------------------------------------------
plot(thrombo_nodesplit)
## ----thrombo_nodesplit_omega, fig.height = 6------------------------------------------------------
plot(thrombo_nodesplit, pars = "omega", stat = "halfeye", ref_line = 0) +
ggplot2::aes(y = comparison) +
ggplot2::facet_null()
## ----thrombo_releff-------------------------------------------------------------------------------
(thrombo_releff <- relative_effects(thrombo_fit, all_contrasts = TRUE))
plot(thrombo_releff, ref_line = 0)
## ----thrombo_ranks--------------------------------------------------------------------------------
(thrombo_ranks <- posterior_ranks(thrombo_fit))
plot(thrombo_ranks)
## ----thrombo_rankprobs----------------------------------------------------------------------------
(thrombo_rankprobs <- posterior_rank_probs(thrombo_fit))
plot(thrombo_rankprobs)
## ----thrombo_cumrankprobs-------------------------------------------------------------------------
(thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE))
plot(thrombo_cumrankprobs)
## ----thrombo_tests, include=FALSE, eval=params$run_tests------------------------------------------
#--- Test against TSD 4 results ---
library(testthat)
library(dplyr)
library(tidyr)
test_that("Reference trt is SK", {
expect_equivalent(levels(thrombo_net$treatments)[1], "SK")
})
tol <- 0.05
tol_dic <- 0.1
# Relative effects
tsd_releff <- tribble(
~trt_b , ~trt , ~mean , ~sd , ~lower, ~upper,
"SK" , "t-PA" , 0.002 , 0.030, -0.06 , 0.06 ,
"SK" , "Acc t-PA" , -0.177, 0.043, -0.26 , -0.09 ,
"SK" , "SK + t-PA", -0.049, 0.046, -0.14 , 0.04 ,
"SK" , "r-PA" , -0.124, 0.060, -0.24 , -0.01 ,
"SK" , "PTCA" , -0.476, 0.101, -0.67 , -0.28 ,
"SK" , "UK" , -0.203, 0.221, -0.64 , 0.23 ,
"SK" , "ASPAC" , 0.016 , 0.037, -0.06 , 0.09 ,
"t-PA" , "PTCA" , -0.478, 0.104, -0.68 , -0.27 ,
"t-PA" , "UK" , -0.206, 0.221, -0.64 , 0.23 ,
"t-PA" , "ASPAC" , 0.013 , 0.037, -0.06 , 0.09 ,
"Acc t-PA", "r-PA" , 0.054 , 0.055, -0.05 , 0.16 ,
"Acc t-PA", "TNK" , 0.005 , 0.064, -0.12 , 0.13 ,
"Acc t-PA", "PTCA" , -0.298, 0.098, -0.49 , -0.11 ,
"Acc t-PA", "UK" , -0.026, 0.221, -0.46 , 0.41 ,
"Acc t-PA", "ASPAC" , 0.193 , 0.056, 0.08 , 0.30 ) %>%
mutate(.trt_b = ordered(trt_b, levels = levels(thrombo_net$treatments)),
.trt = ordered(trt, levels = levels(thrombo_net$treatments)),
rev = if_else(.trt_b > .trt, -1, 1),
.l = lower, .u = upper,
lower = if_else(.trt_b > .trt, .u, .l),
upper = if_else(.trt_b > .trt, .l, .u),
trt_b = if_else(.trt_b > .trt, .trt, .trt_b),
trt = if_else(.trt_b > .trt, .trt_b, .trt),
lab = paste0("d[", trt, " vs. ", trt_b, "]")) %>%
arrange(trt_b, trt) %>%
mutate_at(vars(mean, lower, upper), ~.*rev)
thrombo_releff_summary <- as.data.frame(thrombo_releff) %>%
filter(parameter %in% tsd_releff$lab)
test_that("FE relative effects", {
expect_equivalent(thrombo_releff_summary$mean, tsd_releff$mean, tolerance = tol)
expect_equivalent(thrombo_releff_summary$sd, tsd_releff$sd, tolerance = tol)
expect_equivalent(thrombo_releff_summary$`2.5%`, tsd_releff$lower, tolerance = tol)
expect_equivalent(thrombo_releff_summary$`97.5%`, tsd_releff$upper, tolerance = tol)
})
test_that("SUCRAs", {
thrombo_ranks <- posterior_ranks(thrombo_fit, sucra = TRUE)
thrombo_rankprobs <- posterior_rank_probs(thrombo_fit, sucra = TRUE)
thrombo_cumrankprobs <- posterior_rank_probs(thrombo_fit, cumulative = TRUE, sucra = TRUE)
expect_equal(thrombo_ranks$summary$sucra, thrombo_rankprobs$summary$sucra)
expect_equal(thrombo_ranks$summary$sucra, thrombo_cumrankprobs$summary$sucra)
})
# DIC
test_that("DIC", {
expect_equivalent(dic_consistency$resdev, 105.9, tolerance = tol_dic)
expect_equivalent(dic_consistency$pd, 58, tolerance = tol_dic)
expect_equivalent(dic_consistency$dic, 163.9, tolerance = tol_dic)
})
# Relative effects (UME)
# FE UME model, so no differences by reference treatment, no multi-arm correction
tsd_ume <- tribble(
~trt_b , ~trt , ~mean , ~sd , ~lower, ~upper,
"SK" , "t-PA" , -0.004, 0.030, -0.06 , 0.06 ,
"SK" , "Acc t-PA" , -0.158, 0.049, -0.25 , -0.06 ,
"SK" , "SK + t-PA", -0.044, 0.047, -0.14 , 0.05 ,
"SK" , "r-PA" , -0.060, 0.089, -0.23 , 0.11 ,
"SK" , "PTCA" , -0.665, 0.185, -1.03 , -0.31 ,
"SK" , "UK" , -0.369, 0.518, -1.41 , 0.63 ,
"SK" , "ASPAC" , 0.005 , 0.037, -0.07 , 0.08 ,
"t-PA" , "PTCA" , -0.544, 0.417, -1.38 , 0.25 ,
"t-PA" , "UK" , -0.294, 0.347, -0.99 , 0.37 ,
"t-PA" , "ASPAC" , -0.290, 0.361, -1.01 , 0.41 ,
"Acc t-PA", "r-PA" , 0.019 , 0.066, -0.11 , 0.15 ,
"Acc t-PA", "TNK" , 0.006 , 0.064, -0.12 , 0.13 ,
"Acc t-PA", "PTCA" , -0.216, 0.119, -0.45 , 0.02 ,
"Acc t-PA", "UK" , 0.146 , 0.358, -0.54 , 0.86 ,
"Acc t-PA", "ASPAC" , 1.405 , 0.417, 0.63 , 2.27 ) %>%
mutate(.trt_b = ordered(trt_b, levels = levels(thrombo_net$treatments)),
.trt = ordered(trt, levels = levels(thrombo_net$treatments)),
rev = if_else(.trt_b > .trt, -1, 1),
.l = lower, .u = upper,
lower = if_else(.trt_b > .trt, .u, .l),
upper = if_else(.trt_b > .trt, .l, .u),
trt_b = if_else(.trt_b > .trt, .trt, .trt_b),
trt = if_else(.trt_b > .trt, .trt_b, .trt),
lab = paste0("d[", trt, " vs. ", trt_b, "]")) %>%
arrange(trt_b, trt) %>%
mutate_at(vars(mean, lower, upper), ~.*rev)
thrombo_ume_releff <- summary(thrombo_fit_ume, pars = "d")
test_that("UME relative effects", {
expect_equivalent(thrombo_ume_releff$summary$mean, tsd_ume$mean, tolerance = tol)
expect_equivalent(thrombo_ume_releff$summary$sd, tsd_ume$sd, tolerance = tol)
expect_equivalent(thrombo_ume_releff$summary$`2.5%`, tsd_ume$lower, tolerance = tol)
expect_equivalent(thrombo_ume_releff$summary$`97.5%`, tsd_ume$upper, tolerance = tol)
})
# DIC (UME)
test_that("UME DIC", {
expect_equivalent(dic_ume$resdev, 99.7, tolerance = tol_dic)
expect_equivalent(dic_ume$pd, 65, tolerance = tol_dic)
expect_equivalent(dic_ume$dic, 164.7, tolerance = tol_dic)
})
# Node-splitting
trt_code <- c("SK", "t-PA", "Acc t-PA", "SK + t-PA", "r-PA",
"TNK", "PTCA", "UK", "ASPAC")
dias2010_nodesplit_est <- tribble(
~trt1, ~trt2, ~net_mean, ~net_sd, ~dir_mean, ~dir_sd, ~ind_mean, ~ind_sd, ~omega_mean, ~omega_sd, ~p_value,
1, 2, 0.002, 0.030, 0.000, 0.030, 0.189, 0.235, -0.190, 0.236, 0.42,
1, 3, -0.177, 0.043, -0.158, 0.048, -0.247, 0.092, 0.088, 0.104, 0.39,
1, 5, -0.124, 0.060, -0.060, 0.089, -0.175, 0.081, 0.115, 0.121, 0.34,
1, 7, -0.475, 0.101, -0.666, 0.185, -0.393, 0.120, -0.272, 0.222, 0.22,
1, 8, -0.203, 0.219, -0.369, 0.518, -0.168, 0.244, -0.207, 0.575, 0.73,
1, 9, 0.016, 0.037, 0.009, 0.037, 0.424, 0.252, -0.413, 0.253, 0.10,
2, 7, -0.477, 0.104, -0.545, 0.417, -0.475, 0.108, -0.073, 0.432, 0.88,
2, 8, -0.205, 0.220, -0.295, 0.347, -0.144, 0.290, -0.155, 0.452, 0.74,
2, 9, 0.014, 0.037, 0.006, 0.037, 0.471, 0.241, -0.468, 0.241, 0.05,
3, 4, 0.128, 0.054, 0.126, 0.054, 0.630, 0.697, -0.506, 0.696, 0.47,
3, 5, 0.053, 0.056, 0.019, 0.066, 0.135, 0.101, -0.116, 0.120, 0.34,
3, 7, -0.298, 0.097, -0.216, 0.118, -0.477, 0.174, 0.260, 0.211, 0.21,
3, 8, -0.026, 0.220, 0.143, 0.356, -0.136, 0.288, 0.277, 0.461, 0.55,
3, 9, 0.193, 0.056, 1.409, 0.415, 0.165, 0.057, 1.239, 0.420, 0.001
) %>%
mutate(trt1 = forcats::fct_relevel(factor(trt1, levels = 1:9, labels = trt_code),
!! levels(thrombo_net2$treatments)),
trt2 = forcats::fct_relevel(factor(trt2, levels = 1:9, labels = trt_code),
!! levels(thrombo_net2$treatments)))
for (i in 1:nrow(dias2010_nodesplit_est)) {
if (as.numeric(dias2010_nodesplit_est$trt1[i]) > as.numeric(dias2010_nodesplit_est$trt2[i])) {
dias2010_nodesplit_est[i, c("trt1", "trt2")] <- dias2010_nodesplit_est[i, c("trt2", "trt1")]
dias2010_nodesplit_est[i, c("net_mean", "dir_mean", "ind_mean")] <- -dias2010_nodesplit_est[i, c("net_mean", "dir_mean", "ind_mean")]
}
}
dias2010_nodesplit_est <- arrange(dias2010_nodesplit_est, trt1, trt2)
thrombo_nodesplit_est <- as_tibble(summary(thrombo_nodesplit), nest = FALSE) %>%
mutate(parameter = stringr::str_replace(parameter, "(^d_)?(.+)\\[.+\\]$", "\\2")) %>%
pivot_wider(names_from = "parameter",
values_from = c("mean", "sd"),
names_glue = "{parameter}_{.value}",
id_cols = c("trt1", "trt2", "p_value")) %>%
mutate(across(where(is.numeric), unname)) %>%
select(!!! colnames(dias2010_nodesplit_est))
test_that("Node-splitting estimates", {
# expect_equal(thrombo_nodesplit_est, dias2010_nodesplit_est, tolerance = tol)
expect_equal(thrombo_nodesplit_est$net_mean, dias2010_nodesplit_est$net_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$net_sd, dias2010_nodesplit_est$net_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$ind_mean, dias2010_nodesplit_est$ind_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$ind_sd, dias2010_nodesplit_est$ind_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$dir_mean, dias2010_nodesplit_est$dir_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$dir_sd, dias2010_nodesplit_est$dir_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$omega_mean, dias2010_nodesplit_est$omega_mean, tolerance = tol)
expect_equal(thrombo_nodesplit_est$omega_sd, dias2010_nodesplit_est$omega_sd, tolerance = tol)
expect_equal(thrombo_nodesplit_est$p_value, dias2010_nodesplit_est$p_value, tolerance = tol)
})
dias2010_nodesplit_dic <- tribble(
~trt1, ~trt2, ~resdev, ~pd, ~dic,
1, 2, 106.4, 58.7, 165.1,
1, 3, 106.2, 58.7, 165.0,
1, 5, 106.1, 58.7, 164.8,
1, 7, 105.5, 58.7, 164.2,
1, 8, 106.9, 58.7, 165.6,
1, 9, 104.3, 58.7, 163.0,
2, 7, 106.9, 58.7, 165.6,
2, 8, 106.8, 58.7, 165.5,
2, 9, 103.3, 58.7, 162.0,
3, 4, 106.5, 58.7, 165.2,
3, 5, 106.1, 58.7, 164.8,
3, 7, 105.5, 58.7, 164.2,
3, 8, 106.6, 58.7, 165.3,
3, 9, 96.9, 58.7, 155.6) %>%
mutate(trt1 = forcats::fct_relevel(factor(trt1, levels = 1:9, labels = trt_code),
!! levels(thrombo_net$treatments)),
trt2 = forcats::fct_relevel(factor(trt2, levels = 1:9, labels = trt_code),
!! levels(thrombo_net$treatments)))
for (i in 1:nrow(dias2010_nodesplit_dic)) {
if (as.numeric(dias2010_nodesplit_dic$trt1[i]) > as.numeric(dias2010_nodesplit_dic$trt2[i])) {
dias2010_nodesplit_dic[i, c("trt1", "trt2")] <- dias2010_nodesplit_dic[i, c("trt2", "trt1")]
}
}
dias2010_nodesplit_dic <- arrange(dias2010_nodesplit_dic, trt1, trt2)
thrombo_nodesplit_dic <- as_tibble(summary(thrombo_nodesplit), nest = FALSE) %>%
distinct(trt1, trt2, resdev, pd, dic)
test_that("Node-splitting DIC", {
expect_equal(thrombo_nodesplit_dic$dic, dias2010_nodesplit_dic$dic, tolerance = tol_dic)
expect_equal(thrombo_nodesplit_dic$resdev, dias2010_nodesplit_dic$resdev, tolerance = tol_dic)
expect_equal(thrombo_nodesplit_dic$pd, dias2010_nodesplit_dic$pd, tolerance = tol_dic)
})
# 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.