## ----brms-setup, include=FALSE------------------------------
library(IJSE)
stopifnot(requireNamespace("brms", quietly = TRUE))
options(width = 62)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
se_table <- function(x) {
  round(data.frame(Estimate = coef(x),
    Posterior_SD = x$posterior_sd, IJ_SE = x$se), 3)
}

## ----brms-load----------------------------------------------
path <- system.file("extdata", "tutorial",
  "median-examples.rds", package = "IJSE")
stopifnot(nzchar(path))
examples <- readRDS(path)
ind <- examples$independent
cl <- examples$clustered
stopifnot(inherits(ind$fit, "brmsfit"),
          inherits(cl$fit, "brmsfit"))

## ----brms-check-rows----------------------------------------
fit <- cl$fit
ids <- cl$cluster
stopifnot(length(ids) == nrow(cl$data),
  identical(ids, cl$data$group),
  isTRUE(all.equal(fit$data$x, cl$data$x)),
  isTRUE(all.equal(fit$data$y, cl$data$y)),
  !anyNA(ids), length(unique(ids)) == 30L)

## ----brms-independent---------------------------------------
independent <- ijse(ind$fit, diagnostics = "posterior")
se_table(independent)
confint(independent)
stopifnot(!independent$clustered,
          independent$n_units == 200L)
stopifnot(isTRUE(all.equal(
  IJ_se(ind$fit), independent$se)))

## ----brms-clustered-----------------------------------------
observation_ij <- ijse(fit, diagnostics = "none")
cluster_ij <- ijse(fit, cluster_var = ids,
  save_influence = TRUE, diagnostics = "all")
se_table(cluster_ij)
data.frame(Observation_IJ = observation_ij$se,
           Cluster_IJ = cluster_ij$se)
stopifnot(cluster_ij$clustered,
  cluster_ij$n_obs == 300L, cluster_ij$n_units == 30L,
  isTRUE(all.equal(coef(observation_ij), coef(cluster_ij))))

## ----brms-transform-----------------------------------------
predictions <- ijse(fit, cluster_var = ids,
  variables = c("b_Intercept", "b_x"),
  transform = function(d) {
    cbind(at_minus1 = d[, "b_Intercept"] - d[, "b_x"],
          at_plus1 = d[, "b_Intercept"] + d[, "b_x"],
          difference = 2 * d[, "b_x"])
  }, diagnostics = "posterior")
se_table(predictions)
confint(predictions, parm = "difference", level = 0.90)

## ----brms-contrast-check------------------------------------
a <- c(b_Intercept = 0, b_x = 2)
v <- vcov(cluster_ij)[names(a), names(a), drop = FALSE]
contrast_se <- sqrt(drop(t(a) %*% v %*% a))
stopifnot(isTRUE(all.equal(
  unname(predictions$se["difference"]), contrast_se)))

## ----brms-subset--------------------------------------------
arr <- posterior::as_draws_array(fit,
  variable = c("b_Intercept", "b_x"), regex = FALSE)
n_iter <- posterior::niterations(arr)
n_chain <- posterior::nchains(arr)
keep <- unlist(lapply(seq_len(n_chain), function(j) {
  (j - 1L) * n_iter + seq.int(1L, n_iter, by = 2L)
}))
subset_result <- ijse(fit, cluster_var = ids,
  draw_ids = keep, diagnostics = "posterior")
stopifnot(subset_result$n_draws == length(keep))
data.frame(All_draws = cluster_ij$se,
           Selected_draws = subset_result$se)

## ----brms-stream--------------------------------------------
streamed <- ijse(fit, cluster_var = ids,
  block_size = 40, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(streamed),
                         vcov(cluster_ij))))

## ----brms-diagnostics---------------------------------------
data.frame(Rhat = cluster_ij$diagnostics$rhat,
  Bulk_ESS = cluster_ij$diagnostics$ess_bulk)
cluster_ij$diagnostics$chain_relative_range