## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, message = FALSE)
comma <- function(x) format(x, digits = 2, big.mark = ",")

# Load the pre-computed heavy bootstrap objects to save CRAN build time.
# These were generated locally and saved to the vignettes directory.
load("precomputed_bootstraps.rda")

## -----------------------------------------------------------------------------
library(dplyr)
library(lme4)
library(boot)
library(bootmlm)

## -----------------------------------------------------------------------------
set.seed(85957)
pop_sub <- pop_syn |> filter(school %in% sample(unique(school), 30)) |>
  group_by(school) |> sample_frac(size = .25) |> ungroup()

## -----------------------------------------------------------------------------
m0 <- lmer(popular ~ (1 | school), data = pop_sub)
summary(m0)

## -----------------------------------------------------------------------------
(icc0 <- 1 / (1 + getME(m0, "theta")^(-2)))

## -----------------------------------------------------------------------------
icc <- function(x) 1 / (1 + x@theta^(-2))
icc(m0)

## ----eval=FALSE---------------------------------------------------------------
# system.time(
#   boo_par <- bootstrap_mer(m0, icc, nsim = 999L, type = "parametric")
# )
# boo_par

## ----echo=FALSE---------------------------------------------------------------
# Outputting the pre-computed object
boo_par

## -----------------------------------------------------------------------------
boo_par_ci <- boot::boot.ci(boo_par, type = c("norm", "basic", "perc"), 
                            index = 1L)
boo_par_ci

## ----eval=FALSE---------------------------------------------------------------
# system.time(
#   boo_res <- bootstrap_mer(m0, icc, nsim = 999L, type = "residual")
# )
# boo_res

## ----echo=FALSE---------------------------------------------------------------
boo_res

## ----eval=FALSE---------------------------------------------------------------
# system.time(
#   boo_cgr <- bootstrap_mer(m0, icc, nsim = 999L, type = "residual_cgr")
# )
# boo_cgr

## ----echo=FALSE---------------------------------------------------------------
boo_cgr

## ----eval=FALSE---------------------------------------------------------------
# # Transformation according to V
# system.time(
#   boo_tra <- bootstrap_mer(m0, icc, nsim = 999L, type = "residual_trans")
# )
# boo_tra
# 
# # Transformation according to the sampling variance of r
# system.time(
#   boo_trac <- bootstrap_mer(m0, icc, nsim = 999L, type = "residual_trans",
#                             corrected_trans = TRUE)
# )
# boo_trac

## ----echo=FALSE---------------------------------------------------------------
boo_tra
boo_trac

## -----------------------------------------------------------------------------
# First need to compute the influence values
inf_val <- empinf_mer(m0, icc, index = 1)

# Residual bootstrap
boo_res_ci <- boot::boot.ci(boo_res, type = c("norm", "basic", "perc", "bca"), 
                            index = 1L, L = inf_val)
boo_res_ci

# CGR
boo_cgr_ci <- boot::boot.ci(boo_cgr, type = c("norm", "basic", "perc", "bca"), 
                            index = 1L, L = inf_val)
boo_cgr_ci

# Transformational (no correction)
boo_tra_ci <- boot::boot.ci(boo_tra, type = c("norm", "basic", "perc", "bca"), 
                            index = 1L, L = inf_val)
boo_tra_ci

# Transformational (with correction)
boo_trac_ci <- boot::boot.ci(boo_trac, type = c("norm", "basic", "perc", "bca"), 
                             index = 1L, L = inf_val)
boo_trac_ci

## ----eval=FALSE---------------------------------------------------------------
# system.time(
#   boo_reb <- bootstrap_mer(m0, icc, nsim = 999L, type = "reb")
# )
# boo_reb
# 
# system.time(
#   boo_rebs <- bootstrap_mer(m0, icc, nsim = 999L, type = "reb",
#                             reb_scale = TRUE)
# )
# boo_rebs

## ----echo=FALSE---------------------------------------------------------------
boo_reb
boo_rebs

## -----------------------------------------------------------------------------
# Only sampling clusters
boo_reb_ci <- boot::boot.ci(boo_reb, type = c("norm", "basic", "perc", "bca"), 
                            index = 1L, L = inf_val)
boo_reb_ci

# Transformational (with correction)
boo_rebs_ci <- boot::boot.ci(boo_rebs, type = c("norm", "basic", "perc", "bca"), 
                             index = 1L, L = inf_val)
boo_rebs_ci

## ----eval=FALSE---------------------------------------------------------------
# system.time(
#   boo_cas <- bootstrap_mer(m0, icc, nsim = 999L, type = "case")
# )
# boo_cas
# 
# system.time(
#   boo_cas1 <- bootstrap_mer(m0, icc, nsim = 999L, type = "case",
#                             lv1_resample = TRUE)
# )
# boo_cas1

## ----echo=FALSE---------------------------------------------------------------
boo_cas
boo_cas1

## -----------------------------------------------------------------------------
# Only sampling clusters
boo_cas_ci <- boot::boot.ci(boo_cas, type = c("norm", "basic", "perc", "bca"), 
                            index = 1L, L = inf_val)
boo_cas_ci

# Transformational (with correction)
boo_cas1_ci <- boot::boot.ci(boo_cas1, type = c("norm", "basic", "perc", "bca"), 
                             index = 1L, L = inf_val)
boo_cas1_ci

## -----------------------------------------------------------------------------
boo_names <- c("parametric", "residual", "cgr", "trans", 
               "trans (cor)", "REB", "REB (scaled)", 
               "case (cluster)", "case (c + i)")
boo_lst <- list(boo_par, boo_res, boo_cgr, boo_tra, boo_trac, 
                boo_reb, boo_rebs, boo_cas, boo_cas1)
boo_ci_lst <- list(boo_par_ci, boo_res_ci, boo_cgr_ci, boo_tra_ci, 
                   boo_trac_ci, boo_reb_ci, boo_rebs_ci, boo_cas_ci, 
                   boo_cas1_ci)

get_ci <- function(boo_ci, type) {
  paste0("(", paste(comma(tail(boo_ci[[type]][1, ], 2L)), collapse = ", "), ")")
}

tab <- tibble(boot_type = boo_names, boo = boo_lst, boo_ci = boo_ci_lst) |>
  mutate(sd = sapply(boo, \(x) comma(sd(x$t))), 
         normal = sapply(boo_ci, \(x) get_ci(x, "normal")), 
         basic = sapply(boo_ci, \(x) get_ci(x, "basic")), 
         percentile = sapply(boo_ci, \(x) get_ci(x, "percent")), 
         bca = sapply(boo_ci, \(x) get_ci(x, "bca"))) |>
  select(-boo, -boo_ci)

knitr::kable(tab)