# Generated by vignette example_transfusion.Rmd: do not edit by hand
# Instead edit example_transfusion.Rmd and then run precompile.R

skip_on_cran()



params <-
list(run_tests = FALSE)

## ----code=readLines("children/knitr_setup.R"), include=FALSE--------------------------------------

## ----include=FALSE--------------------------------------------------------------------------------
set.seed(2684319)


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


## -------------------------------------------------------------------------------------------------
tr_net <- set_agd_arm(transfusion, 
                           study = studyc,
                           trt = trtc,
                           r = r, 
                           n = n,
                           trt_ref = "Control")
tr_net


## -------------------------------------------------------------------------------------------------
summary(normal(scale = 100))
summary(half_normal(scale = 5))


## ----eval=FALSE, echo=TRUE------------------------------------------------------------------------
## tr_fit_RE_noninf <- nma(tr_net,
##                         trt_effects = "random",
##                         prior_intercept = normal(scale = 100),
##                         prior_trt = normal(scale = 100),
##                         prior_het = half_normal(scale = 5))

## ----echo=FALSE, eval=!params$run_tests-----------------------------------------------------------
## tr_fit_RE_noninf <- nma(tr_net,
##                         seed = 857369814,
##                         trt_effects = "random",
##                         prior_intercept = normal(scale = 100),
##                         prior_trt = normal(scale = 100),
##                         prior_het = half_normal(scale = 5))

## ----echo=FALSE, eval=params$run_tests------------------------------------------------------------
tr_fit_RE_noninf <- suppressWarnings(nma(tr_net, 
                        seed = 857369814,
                        trt_effects = "random",
                        prior_intercept = normal(scale = 100),
                        prior_trt = normal(scale = 100),
                        prior_het = half_normal(scale = 5),
                        iter = 10000,
                        save_warmup = FALSE))


## -------------------------------------------------------------------------------------------------
tr_fit_RE_noninf


## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(tr_fit_RE_noninf, pars = c("d", "mu", "delta"))


## ----tr_RE_noninf_pp_plot-------------------------------------------------------------------------
plot_prior_posterior(tr_fit_RE_noninf, prior = "het")


## -------------------------------------------------------------------------------------------------
noninf_tau <- as.array(tr_fit_RE_noninf, pars = "tau")
noninf_tausq <- noninf_tau^2
names(noninf_tausq) <- "tausq"
summary(noninf_tausq)


## -------------------------------------------------------------------------------------------------
summary(log_normal(-3.93, 1.51))


## ----echo=TRUE, eval=FALSE------------------------------------------------------------------------
## tr_fit_RE_inf <- nma(tr_net,
##                      trt_effects = "random",
##                      prior_intercept = normal(scale = 100),
##                      prior_trt = normal(scale = 100),
##                      prior_het = log_normal(-3.93, 1.51),
##                      prior_het_type = "var")

## ----echo=FALSE, eval=!params$run_tests-----------------------------------------------------------
## tr_fit_RE_inf <- nma(tr_net,
##                      seed = 1803772660,
##                      trt_effects = "random",
##                      prior_intercept = normal(scale = 100),
##                      prior_trt = normal(scale = 100),
##                      prior_het = log_normal(-3.93, 1.51),
##                      prior_het_type = "var")

## ----echo=FALSE, eval=params$run_tests------------------------------------------------------------
tr_fit_RE_inf <- suppressWarnings(nma(tr_net, 
                     seed = 1803772660,
                     trt_effects = "random",
                     prior_intercept = normal(scale = 100),
                     prior_trt = normal(scale = 100),
                     prior_het = log_normal(-3.93, 1.51),
                     prior_het_type = "var",
                     iter = 10000, save_warmup = FALSE))


## -------------------------------------------------------------------------------------------------
tr_fit_RE_inf


## ----eval=FALSE-----------------------------------------------------------------------------------
## # Not run
## print(tr_fit_RE_inf, pars = c("d", "mu", "delta"))


## ----tr_RE_inf_pp_plot----------------------------------------------------------------------------
plot_prior_posterior(tr_fit_RE_inf, prior = "het")


## -------------------------------------------------------------------------------------------------
inf_tau <- as.array(tr_fit_RE_inf, pars = "tau")
inf_tausq <- inf_tau^2
names(inf_tausq) <- "tausq"
summary(inf_tausq)


## ----transfusion_tests, include=FALSE, eval=params$run_tests--------------------------------------
#--- Test against TSD 2 results ---
library(testthat)

tol <- 0.05

# Non-informative prior
tr_RE_noninf_var <- as.data.frame(summary(noninf_tausq))

test_that("Non-informative RE heterogeneity variance", {
  skip("Non-informative priors not identical")
  expect_equivalent(tr_RE_noninf_var$`50%`, 2.74, tolerance = tol)
  expect_equivalent(tr_RE_noninf_var$`2.5%`, 0.34, tolerance = tol)
  expect_equivalent(tr_RE_noninf_var$`97.5%`, 18.1, tolerance = tol)
})

# Informative prior
tr_RE_inf_var <- as.data.frame(summary(inf_tausq))

test_that("Informative RE heterogeneity variance", {
  expect_equivalent(tr_RE_inf_var$`50%`, 0.18, tolerance = tol)
  expect_equivalent(tr_RE_inf_var$`2.5%`, 0.003, tolerance = tol)
  skip_on_ci()
  expect_equivalent(tr_RE_inf_var$`97.5%`, 1.84, tolerance = tol)
})


# Force clean up
rm(list = ls())
gc()