## ----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_derived.rda")

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

## -----------------------------------------------------------------------------
set.seed(85957)
pop_sub <- pop_syn |> filter(school %in% sample(unique(school), 20)) |>
  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)))

## -----------------------------------------------------------------------------
# Get the point estimate of theta:
th_est <- getME(m0, "theta")
# Get the variance of theta using the bootmlm::vcov_theta() function
th_var <- bootmlm::vcov_theta(m0)
# We can get the asymptotic variance of ICC using delta method with an 
# appropriate transformation. `x1` is the variable name needed for the first
# variable, which is theta in our case
icc_avar <- msm::deltamethod(g = ~ 1 / (1 + x1^(-2)), 
                            mean = th_est, cov = th_var, ses = FALSE)

## -----------------------------------------------------------------------------
icc0 + c(-2, 2) * sqrt(icc_avar)

## -----------------------------------------------------------------------------
icc <- function(x) {
  th_est <- x@theta
  est <- 1 / (1 + th_est^(-2))
  th_var <- vcov_theta(x)
  var <- msm::deltamethod(g = ~ 1 / (1 + x1^(-2)), 
                            mean = th_est, cov = th_var, ses = FALSE)
  c(est, var)
}
icc(m0)

## ----eval=FALSE---------------------------------------------------------------
# # These heavy computations are not run live on CRAN, but you can run them locally.
# boo <- bootstrap_mer(m0, icc, 999L, type = "case")
# 
# # Influence values for BCa
# inf_val <- empinf_mer(m0, icc, index = 1)

## -----------------------------------------------------------------------------
# boo and inf_val are pre-loaded in this vignette
boot::boot.ci(boo, L = inf_val)

## ----out.width="80%", fig.height=4, fig.width=4-------------------------------
hist(boo$t[ , 1])
hist(-log(1 / boo$t[ , 1] - 1))

## -----------------------------------------------------------------------------
boot::boot.ci(boo, L = inf_val, h = qlogis, 
              hdot = function(x) 1 / (x - x^2), 
              hinv = plogis)