## ----diag-setup, include=FALSE------------------------------
library(IJSE)
options(width = 62)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
show_error <- function(expr) {
  message <- tryCatch({force(expr); NULL},
    error = function(e) conditionMessage(e))
  if (is.null(message)) stop("Expected an error in this example.")
  message
}

## ----diag-data----------------------------------------------
set.seed(452)
group <- rep(paste0("g", seq_len(20)), each = 5)
u <- rnorm(20)
y <- 0.5 + sqrt(0.6) * rep(u, each = 5) +
     sqrt(0.4) * rnorm(length(group))
v <- 1 / (length(y) + 1 / 4)
m <- v * sum(y)
mu <- rnorm(2000, m, sqrt(v))
draws <- cbind(mu = mu)
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE),
             numeric(length(mu)))
chain <- rep(c("A", "B"), each = 1000)

## ----diag-modes---------------------------------------------
full <- ijse_draws(draws, ll, cluster_var = group,
  chain_id = chain, save_influence = TRUE)
posterior_only <- ijse_draws(draws, ll, cluster_var = group,
  chain_id = chain, diagnostics = "posterior")
pooled <- ijse_draws(draws, ll, cluster_var = group,
  chain_id = chain, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(full), vcov(posterior_only))),
          isTRUE(all.equal(vcov(full), vcov(pooled))),
          is.null(posterior_only$diagnostics$chain_se),
          is.null(pooled$diagnostics$rhat))

## ----diag-values--------------------------------------------
full$diagnostics$chain_se
full$diagnostics$chain_relative_range
full$diagnostics$chain_ids
c(Rhat = unname(full$diagnostics$rhat),
  Bulk_ESS = unname(full$diagnostics$ess_bulk),
  Tail_ESS = unname(full$diagnostics$ess_tail))

## ----diag-unequal-------------------------------------------
arr <- posterior::as_draws_array(array(mu, c(1000, 2, 1),
  dimnames = list(NULL, NULL, "mu")))
df <- posterior::as_draws_df(arr)
keep <- !(df$.chain == 2 & df$.iteration == 1000)
unequal <- ijse_draws(df[keep, ], ll[keep, , drop = FALSE],
                     cluster_var = group)
stopifnot(unequal$n_draws == 1999L,
  is.null(unequal$diagnostics$rhat),
  nrow(unequal$diagnostics$chain_se) == 2L)

## ----diag-influence-----------------------------------------
case_derivatives <- influence(full) / full$n_units
id <- colnames(case_derivatives)[
  which.max(abs(case_derivatives["mu", ]))]
c(Cluster = id,
  Derivative = format(case_derivatives["mu", id]))

## ----diag-covariance----------------------------------------
stopifnot(isTRUE(all.equal(
  cov(t(influence(full))) / full$n_units, vcov(full))))

## ----diag-derivative----------------------------------------
contribution <- rowSums(ll[, group == id, drop = FALSE])
weighted_mean <- function(delta) {
  log_weight <- delta * contribution
  weight <- exp(log_weight - max(log_weight))
  sum(weight * mu) / sum(weight)
}
h <- 1e-5
finite_difference <- (weighted_mean(h) - weighted_mean(-h)) /
                     (2 * h)
s <- nrow(draws)
adjusted_difference <- finite_difference * s / (s - 1)
c(Covariance = case_derivatives["mu", id],
  Adjusted_difference = adjusted_difference)
stopifnot(isTRUE(all.equal(
  unname(case_derivatives["mu", id]), adjusted_difference,
  tolerance = 1e-7)))

## ----diag-missing-------------------------------------------
bad_ids <- group
bad_ids[3] <- NA_character_
show_error(ijse_draws(draws, ll, cluster_var = bad_ids))
factor_ids <- factor(bad_ids, exclude = NULL)
c(AnyNA_factor = anyNA(factor_ids),
  AnyNA_values = anyNA(as.character(factor_ids)))
show_error(ijse_draws(draws, ll, cluster_var = factor_ids))

## ----diag-dimensions----------------------------------------
show_error(ijse_draws(draws, function(ii) ll[, ii],
  n_obs = length(y), block_size = 1))
fixed <- ijse_draws(draws,
  function(ii) ll[, ii, drop = FALSE], n_obs = length(y),
  block_size = 1, cluster_var = group, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(fixed), vcov(full))))

## ----diag-order---------------------------------------------
named_draws <- draws
named_ll <- ll
rownames(named_draws) <- paste0("s", seq_len(nrow(draws)))
rownames(named_ll) <- rev(rownames(named_draws))
show_error(ijse_draws(named_draws, named_ll))
show_error(ijse_draws(draws, ll, variables = "mu_typo"))
show_error(influence(pooled))

## ----diag-cross-block-order---------------------------------
block_ll <- ll
rownames(block_ll) <- paste0("s", seq_len(nrow(draws)))
misordered_reader <- function(ii) {
  block <- block_ll[, ii, drop = FALSE]
  if (min(ii) > 2) block <- block[nrow(block):1, , drop = FALSE]
  block
}
show_error(ijse_draws(draws, misordered_reader,
  n_obs = ncol(ll), block_size = 2))

## ----diag-units---------------------------------------------
observation_result <- ijse_draws(draws, ll,
  diagnostics = "none")
c(Observation_units = observation_result$n_units,
  Group_units = full$n_units)
stopifnot(isTRUE(all.equal(coef(observation_result),
                         coef(full))), full$n_units == 20L)

## ----diag-permutation---------------------------------------
perm <- rev(seq_along(y))
reordered <- ijse_draws(draws, ll[, perm, drop = FALSE],
  cluster_var = group[perm], diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(reordered), vcov(full))))