# Generated by vignette example_blocker.Rmd: do not edit by hand # Instead edit example_blocker.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(blocker) ## ------------------------------------------------------------------------------------------------- blocker_net <- set_agd_arm(blocker, study = studyn, trt = trtc, r = r, n = n, trt_ref = "Control") blocker_net ## ------------------------------------------------------------------------------------------------- summary(normal(scale = 100)) ## ------------------------------------------------------------------------------------------------- blocker_fit_FE <- nma(blocker_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100)) ## ------------------------------------------------------------------------------------------------- blocker_fit_FE ## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(blocker_fit_FE, pars = c("d", "mu")) ## ----blocker_FE_pp_plot--------------------------------------------------------------------------- plot_prior_posterior(blocker_fit_FE, prior = "trt") ## ------------------------------------------------------------------------------------------------- summary(normal(scale = 100)) summary(half_normal(scale = 5)) ## ------------------------------------------------------------------------------------------------- blocker_fit_RE <- nma(blocker_net, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5)) ## ------------------------------------------------------------------------------------------------- blocker_fit_RE ## ----eval=FALSE----------------------------------------------------------------------------------- ## # Not run ## print(blocker_fit_RE, pars = c("d", "mu", "delta")) ## ----blocker_RE_pp_plot--------------------------------------------------------------------------- plot_prior_posterior(blocker_fit_RE, prior = c("trt", "het")) ## ------------------------------------------------------------------------------------------------- (dic_FE <- dic(blocker_fit_FE)) ## ------------------------------------------------------------------------------------------------- (dic_RE <- dic(blocker_fit_RE)) ## ----blocker_FE_resdev_plot----------------------------------------------------------------------- plot(dic_FE) ## ----blocker_RE_resdev_plot----------------------------------------------------------------------- plot(dic_RE) ## ----blocker_pred_FE, fig.height = 2-------------------------------------------------------------- pred_FE <- predict(blocker_fit_FE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response") pred_FE plot(pred_FE) ## ----blocker_pred_RE, fig.height = 2-------------------------------------------------------------- pred_RE <- predict(blocker_fit_RE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response") pred_RE plot(pred_RE) ## ----blocker_pred_FE_beta, fig.height = 2--------------------------------------------------------- pred_FE_beta <- predict(blocker_fit_FE, baseline = distr(qbeta, 4, 36-4), baseline_type = "response", type = "response") pred_FE_beta plot(pred_FE_beta) ## ----blocker_pred_RE_beta, fig.height = 2--------------------------------------------------------- pred_RE_beta <- predict(blocker_fit_RE, baseline = distr(qbeta, 4, 36-4), baseline_type = "response", type = "response") pred_RE_beta plot(pred_RE_beta) ## ----blocker_tests, include=FALSE, eval=params$run_tests------------------------------------------ #--- Test against TSD 2 results --- library(testthat) library(dplyr) tol <- 0.05 tol_dic <- 0.1 # Relative effects blocker_FE_releff <- as.data.frame(relative_effects(blocker_fit_FE)) test_that("FE relative effects", { expect_equivalent(blocker_FE_releff$mean, -0.26, tolerance = tol) expect_equivalent(blocker_FE_releff$sd, 0.050, tolerance = tol) expect_equivalent(blocker_FE_releff$`2.5%`, -0.36, tolerance = tol) expect_equivalent(blocker_FE_releff$`50%`, -0.26, tolerance = tol) expect_equivalent(blocker_FE_releff$`97.5%`, -0.16, tolerance = tol) }) blocker_RE_releff <- as.data.frame(relative_effects(blocker_fit_RE)) test_that("RE relative effects", { expect_equivalent(blocker_RE_releff$mean, -0.25, tolerance = tol) expect_equivalent(blocker_RE_releff$sd, 0.066, tolerance = tol) expect_equivalent(blocker_RE_releff$`2.5%`, -0.38, tolerance = tol) expect_equivalent(blocker_RE_releff$`50%`, -0.25, tolerance = tol) expect_equivalent(blocker_RE_releff$`97.5%`, -0.12, tolerance = tol) }) # RE heterogeneity SD blocker_RE_sd <- as.data.frame(summary(blocker_fit_RE, pars = "tau")) test_that("RE heterogeneity SD", { expect_equivalent(blocker_RE_sd$mean, 0.14, tolerance = tol) expect_equivalent(blocker_RE_sd$sd, 0.082, tolerance = tol) expect_equivalent(blocker_RE_sd$`2.5%`, 0.01, tolerance = tol) expect_equivalent(blocker_RE_sd$`50%`, 0.13, tolerance = tol) expect_equivalent(blocker_RE_sd$`97.5%`, 0.32, tolerance = tol) }) # DIC test_that("FE DIC", { expect_equivalent(dic_FE$resdev, 46.8, tolerance = tol_dic) expect_equivalent(dic_FE$pd, 23.0, tolerance = tol_dic) expect_equivalent(dic_FE$dic, 69.8, tolerance = tol_dic) }) test_that("RE DIC", { expect_equivalent(dic_RE$resdev, 41.9, tolerance = tol_dic) expect_equivalent(dic_RE$pd, 28.1, tolerance = tol_dic) expect_equivalent(dic_RE$dic, 70.0, tolerance = tol_dic) }) # Predictions blocker_pred_FE <- as.data.frame(pred_FE) test_that("FE predicted probabilities", { expect_equivalent(blocker_pred_FE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$sd, c(0.055, 0.045), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) blocker_pred_RE <- as.data.frame(pred_RE) test_that("RE predicted probabilities", { expect_equivalent(blocker_pred_RE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$sd, c(0.055, 0.046), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) # Check predictions with Beta distribution on baseline probability blocker_predbeta_FE <- as.data.frame(pred_FE_beta) test_that("FE predicted probabilities (Beta distribution)", { expect_equal(blocker_pred_FE$mean, blocker_predbeta_FE$mean, tolerance = tol) expect_equal(blocker_pred_FE$sd, blocker_predbeta_FE$sd, tolerance = tol) expect_equal(blocker_pred_FE$`2.5%`, blocker_predbeta_FE$`2.5%`, tolerance = tol) expect_equal(blocker_pred_FE$`50%`, blocker_predbeta_FE$`50%`, tolerance = tol) expect_equal(blocker_pred_FE$`97.5%`, blocker_predbeta_FE$`97.5%`, tolerance = tol) }) blocker_predbeta_RE <- as.data.frame(pred_RE_beta) test_that("RE predicted probabilities (Beta distribution)", { expect_equal(blocker_pred_RE$mean, blocker_predbeta_RE$mean, tolerance = tol) expect_equal(blocker_pred_RE$sd, blocker_predbeta_RE$sd, tolerance = tol) expect_equal(blocker_pred_RE$`2.5%`, blocker_predbeta_RE$`2.5%`, tolerance = tol) expect_equal(blocker_pred_RE$`50%`, blocker_predbeta_RE$`50%`, tolerance = tol) expect_equal(blocker_pred_RE$`97.5%`, blocker_predbeta_RE$`97.5%`, tolerance = tol) }) # Test that ordered multinomial model is equivalent blocker_ord_net <- set_agd_arm(blocker, study = studyn, trt = trtc, r = multi(nonevents = n, events = r, inclusive = TRUE), trt_ref = "Control") blocker_ord_fit_FE <- nma(blocker_ord_net, trt_effects = "fixed", link = "logit", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_aux = flat()) blocker_ord_fit_RE <- nma(blocker_ord_net, trt_effects = "random", link = "logit", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), prior_aux = flat()) blocker_ord_FE_releff <- as.data.frame(relative_effects(blocker_ord_fit_FE)) test_that("Equivalent ordered multinomial FE relative effects", { expect_equivalent(blocker_ord_FE_releff$mean, -0.26, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$sd, 0.050, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$`2.5%`, -0.36, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$`50%`, -0.26, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$`97.5%`, -0.16, tolerance = tol) }) blocker_ord_RE_releff <- as.data.frame(relative_effects(blocker_ord_fit_RE)) test_that("Equivalent ordered multinomial RE relative effects", { expect_equivalent(blocker_ord_RE_releff$mean, -0.25, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$sd, 0.066, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$`2.5%`, -0.38, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$`50%`, -0.25, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$`97.5%`, -0.12, tolerance = tol) }) blocker_ord_RE_sd <- as.data.frame(summary(blocker_ord_fit_RE, pars = "tau")) test_that("Equivalent ordered multinomial RE heterogeneity SD", { expect_equivalent(blocker_ord_RE_sd$mean, 0.14, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$sd, 0.082, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$`2.5%`, 0.01, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$`50%`, 0.13, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$`97.5%`, 0.32, tolerance = tol) }) test_that("Equivalent ordered multinomial FE DIC", { expect_equivalent(dic_FE$resdev, 46.8, tolerance = tol_dic) expect_equivalent(dic_FE$pd, 23.0, tolerance = tol_dic) expect_equivalent(dic_FE$dic, 69.8, tolerance = tol_dic) }) test_that("Equivalent ordered multinomial RE DIC", { expect_equivalent(dic_RE$resdev, 41.9, tolerance = tol_dic) expect_equivalent(dic_RE$pd, 28.1, tolerance = tol_dic) expect_equivalent(dic_RE$dic, 70.0, tolerance = tol_dic) }) blocker_ord_pred_FE <- as.data.frame(predict(blocker_ord_fit_FE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response")) test_that("Equivalent ordered multinomial FE predicted probabilities", { expect_equivalent(blocker_ord_pred_FE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$sd, c(0.055, 0.045), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) blocker_ord_pred_RE <- as.data.frame(predict(blocker_ord_fit_RE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response")) test_that("Equivalent ordered multinomial RE predicted probabilities", { expect_equivalent(blocker_ord_pred_RE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$sd, c(0.055, 0.046), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) # Check predictions with Beta distribution on baseline probability blocker_ord_predbeta_FE <- predict(blocker_ord_fit_FE, baseline = distr(qbeta, 4, 36 - 4), #3.66565, 36.74819 - 3.66565), baseline_type = "response", type = "response") %>% as.data.frame() test_that("FE ordered multinomial predicted probabilities (Beta distribution)", { expect_equal(blocker_ord_pred_FE$mean, blocker_ord_predbeta_FE$mean, tolerance = tol) expect_equal(blocker_ord_pred_FE$sd, blocker_ord_predbeta_FE$sd, tolerance = tol) expect_equal(blocker_ord_pred_FE$`2.5%`, blocker_ord_predbeta_FE$`2.5%`, tolerance = tol) expect_equal(blocker_ord_pred_FE$`50%`, blocker_ord_predbeta_FE$`50%`, tolerance = tol) expect_equal(blocker_ord_pred_FE$`97.5%`, blocker_ord_predbeta_FE$`97.5%`, tolerance = tol) }) blocker_ord_predbeta_RE <- predict(blocker_ord_fit_RE, baseline = distr(qbeta, 4, 36 - 4), #3.66565, 36.74819 - 3.66565), baseline_type = "response", type = "response") %>% as.data.frame() test_that("RE ordered multinomial predicted probabilities (Beta distribution)", { expect_equal(blocker_ord_pred_RE$mean, blocker_ord_predbeta_RE$mean, tolerance = tol) expect_equal(blocker_ord_pred_RE$sd, blocker_ord_predbeta_RE$sd, tolerance = tol) expect_equal(blocker_ord_pred_RE$`2.5%`, blocker_ord_predbeta_RE$`2.5%`, tolerance = tol) expect_equal(blocker_ord_pred_RE$`50%`, blocker_ord_predbeta_RE$`50%`, tolerance = tol) expect_equal(blocker_ord_pred_RE$`97.5%`, blocker_ord_predbeta_RE$`97.5%`, tolerance = tol) }) # Force clean up rm(list = ls()) gc()