tests/testthat/test-example_diabetes.R
In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data
# Generated by vignette example_diabetes.Rmd: do not edit by hand
# Instead edit example_diabetes.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(diabetes)
## -------------------------------------------------------------------------------------------------
db_net <- set_agd_arm(diabetes,
study = studyc,
trt = trtc,
r = r,
n = n)
db_net
## ----eval=FALSE-----------------------------------------------------------------------------------
## plot(db_net, weight_edges = TRUE, weight_nodes = TRUE)
## ----diabetes_network_plot, echo=FALSE------------------------------------------------------------
plot(db_net, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(legend.box.margin = ggplot2::unit(c(0, 0, 0, 4), "lines"))
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
## -------------------------------------------------------------------------------------------------
db_fit_FE <- nma(db_net,
trt_effects = "fixed",
link = "cloglog",
regression = ~offset(log(time)),
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
## -------------------------------------------------------------------------------------------------
db_fit_FE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(db_fit_FE, pars = c("d", "mu"))
## ----db_FE_pp_plot, fig.width=8, fig.height=6, out.width="100%"-----------------------------------
plot_prior_posterior(db_fit_FE)
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
summary(half_normal(scale = 5))
## -------------------------------------------------------------------------------------------------
db_fit_RE <- nma(db_net,
trt_effects = "random",
link = "cloglog",
regression = ~offset(log(time)),
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_het = half_normal(scale = 5),
init_r = 0.5)
## -------------------------------------------------------------------------------------------------
db_fit_RE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(db_fit_RE, pars = c("d", "mu", "delta"))
## ----db_RE_pp_plot--------------------------------------------------------------------------------
plot_prior_posterior(db_fit_RE, prior = c("trt", "het"))
## -------------------------------------------------------------------------------------------------
(dic_FE <- dic(db_fit_FE))
## -------------------------------------------------------------------------------------------------
(dic_RE <- dic(db_fit_RE))
## ----db_FE_resdev_plot----------------------------------------------------------------------------
plot(dic_FE)
## ----db_RE_resdev_plot----------------------------------------------------------------------------
plot(dic_RE)
## ----diabetes_releff_FE, fig.height=3-------------------------------------------------------------
(db_releff_FE <- relative_effects(db_fit_FE, trt_ref = "Diuretic"))
plot(db_releff_FE, ref_line = 0)
## ----diabetes_releff_RE, fig.height=3-------------------------------------------------------------
(db_releff_RE <- relative_effects(db_fit_RE, trt_ref = "Diuretic"))
plot(db_releff_RE, ref_line = 0)
## ----db_pred_FE, fig.height = 2-------------------------------------------------------------------
db_pred_FE <- predict(db_fit_FE,
newdata = data.frame(time = 3),
baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
baseline_trt = "Diuretic",
type = "response")
db_pred_FE
plot(db_pred_FE)
## ----db_pred_RE, fig.height = 2-------------------------------------------------------------------
db_pred_RE <- predict(db_fit_RE,
newdata = data.frame(time = 3),
baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
baseline_trt = "Diuretic",
type = "response")
db_pred_RE
plot(db_pred_RE)
## ----db_pred_RE_all, fig.height=16----------------------------------------------------------------
db_pred_RE_studies <- predict(db_fit_RE, type = "response")
db_pred_RE_studies
plot(db_pred_RE_studies)
## ----diabetes_ranks-------------------------------------------------------------------------------
(db_ranks <- posterior_ranks(db_fit_RE))
plot(db_ranks)
## ----diabetes_rankprobs---------------------------------------------------------------------------
(db_rankprobs <- posterior_rank_probs(db_fit_RE))
plot(db_rankprobs)
## ----diabetes_cumrankprobs------------------------------------------------------------------------
(db_cumrankprobs <- posterior_rank_probs(db_fit_RE, cumulative = TRUE))
plot(db_cumrankprobs)
## ----diabetes_tests, include=FALSE, eval=params$run_tests-----------------------------------------
#--- Test against TSD 2 results ---
library(testthat)
library(dplyr)
tol <- 0.05
tol_dic <- 0.1
# TSD treatment codes
# trt_codes <- c(1 = "Diuretic",
# 2 = "Placebo",
# 3 = "Beta Blocker",
# 4 = "CCB",
# 5 = "ACE Inhibitor",
# 6 = "ARB")
# FE Relative effects
tsd_FE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Placebo" , -0.25, 0.06, -0.25 , -0.36 , -0.14 ,
"Beta Blocker" , -0.06, 0.06, -0.06 , -0.17 , 0.05 ,
"CCB" , -0.25, 0.05, -0.25 , -0.36 , -0.15 ,
"ACE Inhibitor", -0.36, 0.05, -0.36 , -0.46 , -0.25 ,
"ARB" , -0.45, 0.06, -0.45 , -0.58 , -0.33 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_releff_FE <- as.data.frame(db_releff_FE)
test_that("FE relative effects", {
expect_equivalent(db_releff_FE$mean, tsd_FE$mean, tolerance = tol)
expect_equivalent(db_releff_FE$sd, tsd_FE$sd, tolerance = tol)
expect_equivalent(db_releff_FE$`50%`, tsd_FE$median, tolerance = tol)
expect_equivalent(db_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol)
expect_equivalent(db_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol)
})
# RE Relative effects
tsd_RE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Placebo" , -0.29, 0.09, -0.29 , -0.47 , -0.12 ,
"Beta Blocker" , -0.07, 0.09, -0.07 , -0.25 , 0.10 ,
"CCB" , -0.24, 0.08, -0.24 , -0.41 , -0.08 ,
"ACE Inhibitor", -0.40, 0.09, -0.40 , -0.58 , -0.24 ,
"ARB" , -0.47, 0.11, -0.47 , -0.70 , -0.27 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_releff_RE <- as.data.frame(db_releff_RE)
test_that("RE relative effects", {
expect_equivalent(db_releff_RE$mean, tsd_RE$mean, tolerance = tol)
expect_equivalent(db_releff_RE$sd, tsd_RE$sd, tolerance = tol)
expect_equivalent(db_releff_RE$`50%`, tsd_RE$median, tolerance = tol)
expect_equivalent(db_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol)
expect_equivalent(db_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol)
})
# Heterogeneity SD
db_tau <- as.data.frame(summary(db_fit_RE, pars = "tau"))
test_that("RE heterogeneity SD", {
expect_equivalent(db_tau$mean, 0.13, tolerance = tol)
expect_equivalent(db_tau$sd, 0.04, tolerance = tol)
expect_equivalent(db_tau$`50%`, 0.12, tolerance = tol)
expect_equivalent(db_tau$`2.5%`, 0.05, tolerance = tol)
expect_equivalent(db_tau$`97.5%`, 0.23, tolerance = tol)
})
# FE probabilities
tsd_pred_FE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Diuretic" , 0.065, 0.067, 0.044 , 0.01 ,0.25 ,
"Placebo" , 0.052, 0.055, 0.034 , 0.01 ,0.20 ,
"Beta Blocker" , 0.062, 0.064, 0.042 , 0.01 ,0.24 ,
"CCB" , 0.051, 0.055, 0.034 , 0.01 ,0.20 ,
"ACE Inhibitor", 0.047, 0.050, 0.031 , 0.00 ,0.18 ,
"ARB" , 0.043, 0.046, 0.028 , 0.00 ,0.17 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_pred_FE <- as.data.frame(db_pred_FE)
test_that("FE predicted probabilities at 3 years", {
expect_equivalent(db_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol)
expect_equivalent(db_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol)
expect_equivalent(db_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol)
expect_equivalent(db_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol)
expect_equivalent(db_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol)
})
# RE probabilities
tsd_pred_RE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Diuretic" , 0.065, 0.067, 0.044 , 0.01 ,0.25 ,
"Placebo" , 0.050, 0.053, 0.033 , 0.01 ,0.20 ,
"Beta Blocker" , 0.061, 0.064, 0.041 , 0.01 ,0.24 ,
"CCB" , 0.052, 0.056, 0.035 , 0.01 ,0.20 ,
"ACE Inhibitor", 0.045, 0.048, 0.030 , 0.00 ,0.18 ,
"ARB" , 0.042, 0.046, 0.028 , 0.00 ,0.17 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_pred_RE <- as.data.frame(db_pred_RE)
test_that("RE predicted probabilities at 3 years", {
expect_equivalent(db_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol)
expect_equivalent(db_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol)
expect_equivalent(db_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol)
expect_equivalent(db_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol)
expect_equivalent(db_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(dic_RE$resdev, 53.7, tolerance = tol_dic)
expect_equivalent(dic_RE$pd, 38.0, tolerance = tol_dic)
expect_equivalent(dic_RE$dic, 91.7, tolerance = tol_dic)
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(dic_FE$resdev, 78.25, tolerance = tol_dic)
expect_equivalent(dic_FE$pd, 27.0, tolerance = tol_dic)
expect_equivalent(dic_FE$dic, 105.2, tolerance = tol_dic)
})
# Check that predictions for multiple studies works in any order
times <- 1:3
# Baselines named by time
bls <- list("1" = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
"2" = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
"3" = distr(qnorm, mean = -4.2, sd = 1.11^-0.5))
db_pred_FE_multi1 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = times),
baseline = unname(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_RE_multi1 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = times),
baseline = unname(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_FE_multi2 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = times),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_RE_multi2 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = times),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_FE_multi3 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = rev(times)),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response")) %>%
arrange(.study)
db_pred_RE_multi3 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = rev(times)),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response")) %>%
arrange(.study)
db_pred_FE_multi4 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = times),
baseline = rev(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_RE_multi4 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = times),
baseline = rev(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
test_that("Predictions for reordered newdata/baselines works", {
expect_equivalent(db_pred_FE_multi1$mean, db_pred_FE_multi2$mean, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$sd, db_pred_FE_multi2$sd, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`50%`, db_pred_FE_multi2$`50%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`2.5%`, db_pred_FE_multi2$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`97.5%`, db_pred_FE_multi2$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$mean, db_pred_FE_multi3$mean, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$sd, db_pred_FE_multi3$sd, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`50%`, db_pred_FE_multi3$`50%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`2.5%`, db_pred_FE_multi3$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`97.5%`, db_pred_FE_multi3$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$mean, db_pred_FE_multi4$mean, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$sd, db_pred_FE_multi4$sd, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`50%`, db_pred_FE_multi4$`50%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`2.5%`, db_pred_FE_multi4$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`97.5%`, db_pred_FE_multi4$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$mean, db_pred_RE_multi2$mean, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$sd, db_pred_RE_multi2$sd, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`50%`, db_pred_RE_multi2$`50%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`2.5%`, db_pred_RE_multi2$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`97.5%`, db_pred_RE_multi2$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$mean, db_pred_RE_multi3$mean, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$sd, db_pred_RE_multi3$sd, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`50%`, db_pred_RE_multi3$`50%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`2.5%`, db_pred_RE_multi3$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`97.5%`, db_pred_RE_multi3$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$mean, db_pred_RE_multi4$mean, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$sd, db_pred_RE_multi4$sd, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`50%`, db_pred_RE_multi4$`50%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`2.5%`, db_pred_RE_multi4$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`97.5%`, db_pred_RE_multi4$`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_diabetes.Rmd: do not edit by hand
# Instead edit example_diabetes.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(diabetes)
## -------------------------------------------------------------------------------------------------
db_net <- set_agd_arm(diabetes,
study = studyc,
trt = trtc,
r = r,
n = n)
db_net
## ----eval=FALSE-----------------------------------------------------------------------------------
## plot(db_net, weight_edges = TRUE, weight_nodes = TRUE)
## ----diabetes_network_plot, echo=FALSE------------------------------------------------------------
plot(db_net, weight_edges = TRUE, weight_nodes = TRUE) + ggplot2::theme(legend.box.margin = ggplot2::unit(c(0, 0, 0, 4), "lines"))
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
## -------------------------------------------------------------------------------------------------
db_fit_FE <- nma(db_net,
trt_effects = "fixed",
link = "cloglog",
regression = ~offset(log(time)),
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
## -------------------------------------------------------------------------------------------------
db_fit_FE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(db_fit_FE, pars = c("d", "mu"))
## ----db_FE_pp_plot, fig.width=8, fig.height=6, out.width="100%"-----------------------------------
plot_prior_posterior(db_fit_FE)
## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
summary(half_normal(scale = 5))
## -------------------------------------------------------------------------------------------------
db_fit_RE <- nma(db_net,
trt_effects = "random",
link = "cloglog",
regression = ~offset(log(time)),
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_het = half_normal(scale = 5),
init_r = 0.5)
## -------------------------------------------------------------------------------------------------
db_fit_RE
## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(db_fit_RE, pars = c("d", "mu", "delta"))
## ----db_RE_pp_plot--------------------------------------------------------------------------------
plot_prior_posterior(db_fit_RE, prior = c("trt", "het"))
## -------------------------------------------------------------------------------------------------
(dic_FE <- dic(db_fit_FE))
## -------------------------------------------------------------------------------------------------
(dic_RE <- dic(db_fit_RE))
## ----db_FE_resdev_plot----------------------------------------------------------------------------
plot(dic_FE)
## ----db_RE_resdev_plot----------------------------------------------------------------------------
plot(dic_RE)
## ----diabetes_releff_FE, fig.height=3-------------------------------------------------------------
(db_releff_FE <- relative_effects(db_fit_FE, trt_ref = "Diuretic"))
plot(db_releff_FE, ref_line = 0)
## ----diabetes_releff_RE, fig.height=3-------------------------------------------------------------
(db_releff_RE <- relative_effects(db_fit_RE, trt_ref = "Diuretic"))
plot(db_releff_RE, ref_line = 0)
## ----db_pred_FE, fig.height = 2-------------------------------------------------------------------
db_pred_FE <- predict(db_fit_FE,
newdata = data.frame(time = 3),
baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
baseline_trt = "Diuretic",
type = "response")
db_pred_FE
plot(db_pred_FE)
## ----db_pred_RE, fig.height = 2-------------------------------------------------------------------
db_pred_RE <- predict(db_fit_RE,
newdata = data.frame(time = 3),
baseline = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
baseline_trt = "Diuretic",
type = "response")
db_pred_RE
plot(db_pred_RE)
## ----db_pred_RE_all, fig.height=16----------------------------------------------------------------
db_pred_RE_studies <- predict(db_fit_RE, type = "response")
db_pred_RE_studies
plot(db_pred_RE_studies)
## ----diabetes_ranks-------------------------------------------------------------------------------
(db_ranks <- posterior_ranks(db_fit_RE))
plot(db_ranks)
## ----diabetes_rankprobs---------------------------------------------------------------------------
(db_rankprobs <- posterior_rank_probs(db_fit_RE))
plot(db_rankprobs)
## ----diabetes_cumrankprobs------------------------------------------------------------------------
(db_cumrankprobs <- posterior_rank_probs(db_fit_RE, cumulative = TRUE))
plot(db_cumrankprobs)
## ----diabetes_tests, include=FALSE, eval=params$run_tests-----------------------------------------
#--- Test against TSD 2 results ---
library(testthat)
library(dplyr)
tol <- 0.05
tol_dic <- 0.1
# TSD treatment codes
# trt_codes <- c(1 = "Diuretic",
# 2 = "Placebo",
# 3 = "Beta Blocker",
# 4 = "CCB",
# 5 = "ACE Inhibitor",
# 6 = "ARB")
# FE Relative effects
tsd_FE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Placebo" , -0.25, 0.06, -0.25 , -0.36 , -0.14 ,
"Beta Blocker" , -0.06, 0.06, -0.06 , -0.17 , 0.05 ,
"CCB" , -0.25, 0.05, -0.25 , -0.36 , -0.15 ,
"ACE Inhibitor", -0.36, 0.05, -0.36 , -0.46 , -0.25 ,
"ARB" , -0.45, 0.06, -0.45 , -0.58 , -0.33 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_releff_FE <- as.data.frame(db_releff_FE)
test_that("FE relative effects", {
expect_equivalent(db_releff_FE$mean, tsd_FE$mean, tolerance = tol)
expect_equivalent(db_releff_FE$sd, tsd_FE$sd, tolerance = tol)
expect_equivalent(db_releff_FE$`50%`, tsd_FE$median, tolerance = tol)
expect_equivalent(db_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol)
expect_equivalent(db_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol)
})
# RE Relative effects
tsd_RE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Placebo" , -0.29, 0.09, -0.29 , -0.47 , -0.12 ,
"Beta Blocker" , -0.07, 0.09, -0.07 , -0.25 , 0.10 ,
"CCB" , -0.24, 0.08, -0.24 , -0.41 , -0.08 ,
"ACE Inhibitor", -0.40, 0.09, -0.40 , -0.58 , -0.24 ,
"ARB" , -0.47, 0.11, -0.47 , -0.70 , -0.27 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_releff_RE <- as.data.frame(db_releff_RE)
test_that("RE relative effects", {
expect_equivalent(db_releff_RE$mean, tsd_RE$mean, tolerance = tol)
expect_equivalent(db_releff_RE$sd, tsd_RE$sd, tolerance = tol)
expect_equivalent(db_releff_RE$`50%`, tsd_RE$median, tolerance = tol)
expect_equivalent(db_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol)
expect_equivalent(db_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol)
})
# Heterogeneity SD
db_tau <- as.data.frame(summary(db_fit_RE, pars = "tau"))
test_that("RE heterogeneity SD", {
expect_equivalent(db_tau$mean, 0.13, tolerance = tol)
expect_equivalent(db_tau$sd, 0.04, tolerance = tol)
expect_equivalent(db_tau$`50%`, 0.12, tolerance = tol)
expect_equivalent(db_tau$`2.5%`, 0.05, tolerance = tol)
expect_equivalent(db_tau$`97.5%`, 0.23, tolerance = tol)
})
# FE probabilities
tsd_pred_FE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Diuretic" , 0.065, 0.067, 0.044 , 0.01 ,0.25 ,
"Placebo" , 0.052, 0.055, 0.034 , 0.01 ,0.20 ,
"Beta Blocker" , 0.062, 0.064, 0.042 , 0.01 ,0.24 ,
"CCB" , 0.051, 0.055, 0.034 , 0.01 ,0.20 ,
"ACE Inhibitor", 0.047, 0.050, 0.031 , 0.00 ,0.18 ,
"ARB" , 0.043, 0.046, 0.028 , 0.00 ,0.17 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_pred_FE <- as.data.frame(db_pred_FE)
test_that("FE predicted probabilities at 3 years", {
expect_equivalent(db_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol)
expect_equivalent(db_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol)
expect_equivalent(db_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol)
expect_equivalent(db_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol)
expect_equivalent(db_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol)
})
# RE probabilities
tsd_pred_RE <- tribble(
~trt , ~mean, ~sd , ~median, ~lower, ~upper,
"Diuretic" , 0.065, 0.067, 0.044 , 0.01 ,0.25 ,
"Placebo" , 0.050, 0.053, 0.033 , 0.01 ,0.20 ,
"Beta Blocker" , 0.061, 0.064, 0.041 , 0.01 ,0.24 ,
"CCB" , 0.052, 0.056, 0.035 , 0.01 ,0.20 ,
"ACE Inhibitor", 0.045, 0.048, 0.030 , 0.00 ,0.18 ,
"ARB" , 0.042, 0.046, 0.028 , 0.00 ,0.17 ,
) %>%
mutate(trt = ordered(trt, levels = levels(db_net$treatments))) %>%
arrange(trt)
db_pred_RE <- as.data.frame(db_pred_RE)
test_that("RE predicted probabilities at 3 years", {
expect_equivalent(db_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol)
expect_equivalent(db_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol)
expect_equivalent(db_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol)
expect_equivalent(db_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol)
expect_equivalent(db_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(dic_RE$resdev, 53.7, tolerance = tol_dic)
expect_equivalent(dic_RE$pd, 38.0, tolerance = tol_dic)
expect_equivalent(dic_RE$dic, 91.7, tolerance = tol_dic)
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(dic_FE$resdev, 78.25, tolerance = tol_dic)
expect_equivalent(dic_FE$pd, 27.0, tolerance = tol_dic)
expect_equivalent(dic_FE$dic, 105.2, tolerance = tol_dic)
})
# Check that predictions for multiple studies works in any order
times <- 1:3
# Baselines named by time
bls <- list("1" = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
"2" = distr(qnorm, mean = -4.2, sd = 1.11^-0.5),
"3" = distr(qnorm, mean = -4.2, sd = 1.11^-0.5))
db_pred_FE_multi1 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = times),
baseline = unname(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_RE_multi1 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = times),
baseline = unname(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_FE_multi2 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = times),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_RE_multi2 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = times),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_FE_multi3 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = rev(times)),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response")) %>%
arrange(.study)
db_pred_RE_multi3 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = rev(times)),
baseline = bls,
study = time,
baseline_trt = "Diuretic",
type = "response")) %>%
arrange(.study)
db_pred_FE_multi4 <- as.data.frame(predict(db_fit_FE,
newdata = data.frame(time = times),
baseline = rev(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
db_pred_RE_multi4 <- as.data.frame(predict(db_fit_RE,
newdata = data.frame(time = times),
baseline = rev(bls),
study = time,
baseline_trt = "Diuretic",
type = "response"))
test_that("Predictions for reordered newdata/baselines works", {
expect_equivalent(db_pred_FE_multi1$mean, db_pred_FE_multi2$mean, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$sd, db_pred_FE_multi2$sd, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`50%`, db_pred_FE_multi2$`50%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`2.5%`, db_pred_FE_multi2$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`97.5%`, db_pred_FE_multi2$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$mean, db_pred_FE_multi3$mean, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$sd, db_pred_FE_multi3$sd, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`50%`, db_pred_FE_multi3$`50%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`2.5%`, db_pred_FE_multi3$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`97.5%`, db_pred_FE_multi3$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$mean, db_pred_FE_multi4$mean, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$sd, db_pred_FE_multi4$sd, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`50%`, db_pred_FE_multi4$`50%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`2.5%`, db_pred_FE_multi4$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_FE_multi1$`97.5%`, db_pred_FE_multi4$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$mean, db_pred_RE_multi2$mean, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$sd, db_pred_RE_multi2$sd, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`50%`, db_pred_RE_multi2$`50%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`2.5%`, db_pred_RE_multi2$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`97.5%`, db_pred_RE_multi2$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$mean, db_pred_RE_multi3$mean, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$sd, db_pred_RE_multi3$sd, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`50%`, db_pred_RE_multi3$`50%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`2.5%`, db_pred_RE_multi3$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`97.5%`, db_pred_RE_multi3$`97.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$mean, db_pred_RE_multi4$mean, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$sd, db_pred_RE_multi4$sd, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`50%`, db_pred_RE_multi4$`50%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`2.5%`, db_pred_RE_multi4$`2.5%`, tolerance = tol)
expect_equivalent(db_pred_RE_multi1$`97.5%`, db_pred_RE_multi4$`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.