# Generated by vignette example_statins.Rmd: do not edit by hand # Instead edit example_statins.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(statins) ## ------------------------------------------------------------------------------------------------- statin_net <- set_agd_arm(statins, study = studyc, trt = trtc, r = r, n = n, trt_ref = "Placebo") statin_net ## ------------------------------------------------------------------------------------------------- summary(normal(scale = 100)) ## ----eval=!params$run_tests----------------------------------------------------------------------- ## statin_fit_FE <- nma(statin_net, ## trt_effects = "fixed", ## regression = ~.trt:prevention, ## prior_intercept = normal(scale = 100), ## prior_trt = normal(scale = 100), ## prior_reg = normal(scale = 100)) ## ----eval=params$run_tests, echo=FALSE------------------------------------------------------------ statin_fit_FE <- nma(statin_net, trt_effects = "fixed", regression = ~.trt:prevention, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_reg = normal(scale = 100), iter = 5000) ## ------------------------------------------------------------------------------------------------- statin_fit_FE ## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(statin_fit_FE, pars = c("d", "beta", "mu")) ## ----statin_FE_pp_plot---------------------------------------------------------------------------- plot_prior_posterior(statin_fit_FE, prior = c("trt", "reg")) ## ------------------------------------------------------------------------------------------------- summary(normal(scale = 100)) summary(half_normal(scale = 5)) ## ----eval=!params$run_tests----------------------------------------------------------------------- ## statin_fit_RE <- nma(statin_net, ## trt_effects = "random", ## regression = ~.trt:prevention, ## prior_intercept = normal(scale = 100), ## prior_trt = normal(scale = 100), ## prior_reg = normal(scale = 100), ## prior_het = half_normal(scale = 5), ## adapt_delta = 0.99) ## ----eval=params$run_tests, echo=FALSE------------------------------------------------------------ statin_fit_RE <- nowarn_on_ci(nma(statin_net, trt_effects = "random", regression = ~.trt:prevention, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_reg = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99, iter = 5000)) ## ------------------------------------------------------------------------------------------------- statin_fit_RE ## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(statin_fit_RE, pars = c("d", "beta", "mu", "delta")) ## ----statin_RE_pp_plot---------------------------------------------------------------------------- plot_prior_posterior(statin_fit_RE, prior = c("trt", "reg", "het")) ## ------------------------------------------------------------------------------------------------- (statin_dic_FE <- dic(statin_fit_FE)) ## ------------------------------------------------------------------------------------------------- (statin_dic_RE <- dic(statin_fit_RE)) ## ----statin_FE_resdev_plot------------------------------------------------------------------------ plot(statin_dic_FE) ## ----statin_RE_resdev_plot------------------------------------------------------------------------ plot(statin_dic_RE) ## ------------------------------------------------------------------------------------------------- statin_releff_FE <- relative_effects(statin_fit_FE, newdata = data.frame(prevention = c("Primary", "Secondary")), study = prevention) statin_releff_FE ## ----statins_releff_FE, fig.height = 2------------------------------------------------------------ plot(statin_releff_FE, ref_line = 0) ## ----statins_beta_FE, fig.height = 4-------------------------------------------------------------- plot(statin_fit_FE, pars = "beta", ref_line = 0, stat = "halfeye") ## ----statins_tests, include=FALSE, eval=params$run_tests------------------------------------------ #--- Test against TSD 3 results --- library(testthat) library(dplyr) tol <- 0.05 tol_dic <- 0.1 # Relative effects statin_FE_releff <- as.data.frame(summary(statin_fit_FE, pars = "d")) test_that("FE relative effects", { expect_equivalent(statin_FE_releff$mean, -0.11, tolerance = tol) expect_equivalent(statin_FE_releff$sd, 0.10, tolerance = tol) expect_equivalent(statin_FE_releff$`2.5%`, -0.30, tolerance = tol) expect_equivalent(statin_FE_releff$`97.5%`, 0.09, tolerance = tol) }) statin_RE_releff <- as.data.frame(summary(statin_fit_RE, pars = "d")) test_that("RE relative effects", { expect_equivalent(statin_RE_releff$mean, -0.07, tolerance = tol) expect_equivalent(statin_RE_releff$sd, 0.20, tolerance = tol) skip_on_ci() expect_equivalent(statin_RE_releff$`2.5%`, -0.48, tolerance = tol) expect_equivalent(statin_RE_releff$`97.5%`, 0.36, tolerance = tol) }) # Regression coefficients statin_FE_beta <- as.data.frame(summary(statin_fit_FE, pars = "beta")) test_that("FE regression beta", { expect_equivalent(statin_FE_beta$mean, -0.21, tolerance = tol) expect_equivalent(statin_FE_beta$sd, 0.11, tolerance = tol) expect_equivalent(statin_FE_beta$`2.5%`, -0.42, tolerance = tol) expect_equivalent(statin_FE_beta$`97.5%`, 0.01, tolerance = tol) }) statin_RE_beta <- as.data.frame(summary(statin_fit_RE, pars = "beta")) test_that("RE regression beta", { expect_equivalent(statin_RE_beta$mean, -0.29, tolerance = tol) expect_equivalent(statin_RE_beta$sd, 0.26, tolerance = tol) skip_on_ci() expect_equivalent(statin_RE_beta$`2.5%`, -0.86, tolerance = tol) expect_equivalent(statin_RE_beta$`97.5%`, 0.20, tolerance = tol) }) # RE heterogeneity SD statin_RE_sd <- as.data.frame(summary(statin_fit_RE, pars = "tau")) test_that("RE heterogeneity SD", { expect_equivalent(statin_RE_sd$`50%`, 0.19, tolerance = tol) expect_equivalent(statin_RE_sd$sd, 0.20, tolerance = tol) skip_on_ci() expect_equivalent(statin_RE_sd$`2.5%`, 0.01, tolerance = tol) expect_equivalent(statin_RE_sd$`97.5%`, 0.76, tolerance = tol) }) # DIC test_that("FE DIC", { expect_equivalent(statin_dic_FE$resdev, 45.9, tolerance = tol_dic) expect_equivalent(statin_dic_FE$pd, 21.0, tolerance = tol_dic) expect_equivalent(statin_dic_FE$dic, 66.9, tolerance = tol_dic) }) test_that("RE DIC", { expect_equivalent(statin_dic_RE$resdev, 42.6, tolerance = tol_dic) expect_equivalent(statin_dic_RE$pd, 24.2, tolerance = tol_dic) expect_equivalent(statin_dic_RE$dic, 66.8, tolerance = tol_dic) }) test_that("Robust to custom options(contrasts) settings", { withr::with_options(list(contrasts = c(ordered = "contr.SAS", unordered = "contr.SAS")), { statin_fit_FE_SAS <- nma(statin_net, trt_effects = "fixed", regression = ~.trt:prevention, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_reg = normal(scale = 100), iter = 5000) # Model pars are different (reference level of prevention is different) but # relative effects should still be calculated correctly statin_fit_FE_SAS_releff <- as_tibble(relative_effects(statin_fit_FE_SAS))[, c("parameter", "mean", "sd")] }) expect_equal(statin_fit_FE_SAS_releff, as_tibble(relative_effects(statin_fit_FE))[, c("parameter", "mean", "sd")], tolerance = tol) }) # Force clean up rm(list = ls()) gc()