tests/testthat/test-example_blocker.R
In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data
# 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())
# library(ggplot2)
## ----setup, echo = FALSE--------------------------------------------------------------------------
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
"true" =, "warn" = 2,
parallel::detectCores())
options(mc.cores = nc)
library(ggplot2)
## -------------------------------------------------------------------------------------------------
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_leverage_FE--------------------------------------------------------------------------
plot(dic_FE, type = "leverage") +
# Add labels for points outside DIC=3
geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
## ----blocker_leverage_RE--------------------------------------------------------------------------
plot(dic_RE, type = "leverage") +
# Add labels for points outside DIC=3
geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
## ----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()
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_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())
# library(ggplot2)
## ----setup, echo = FALSE--------------------------------------------------------------------------
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
"true" =, "warn" = 2,
parallel::detectCores())
options(mc.cores = nc)
library(ggplot2)
## -------------------------------------------------------------------------------------------------
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_leverage_FE--------------------------------------------------------------------------
plot(dic_FE, type = "leverage") +
# Add labels for points outside DIC=3
geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
## ----blocker_leverage_RE--------------------------------------------------------------------------
plot(dic_RE, type = "leverage") +
# Add labels for points outside DIC=3
geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
## ----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()
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.