# 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()