---
title: "Tutorial 3: Diagnostics, Influence and Troubleshooting"
output:
  rmarkdown::html_vignette:
    css: ijse.css
vignette: >
  %\VignetteIndexEntry{Tutorial 3: Diagnostics, Influence and Troubleshooting}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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
}
```

## What you will learn

This tutorial runs without brms or Stan. You will choose a diagnostic mode,
preserve unequal chain metadata, interpret saved influence values, and reproduce
common input errors. Every error demonstration catches its intended failure, so
the complete script can run successfully.

## Set up clustered observations and a working posterior

Generate 20 independent groups of five with within-group correlation 0.6. Use a
normal working-independence likelihood with observation variance one and prior
variance four. As in Tutorial 1, its posterior is available exactly.

```{r 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)
```

The two chains here are independently generated normal draws, assigned labels
to demonstrate the interface. They are not outputs from an MCMC sampler.

## Choose how much diagnostic computation to do

```{r 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))
```

The default `"all"` calculates chain-specific IJ SEs and quantity-level posterior
diagnostics. `"posterior"` retains R-hat and ESS without additional chain IJ
calculations. `"none"` computes pooled inference only. The modes do not change
the pooled estimator or covariance. For a plain matrix, supply `chain_id` when
you want chain diagnostics; an unlabeled matrix does not contain chain identity.

```{r 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))
```

The relative range is the maximum minus minimum chain IJ SE divided by the pooled
IJ SE. It has no universal pass threshold and is not a calibrated Monte Carlo
SE. R-hat and ESS diagnose the selected quantities, rather than every
parameter-likelihood covariance. For an actual brms fit, inspect sampler
convergence, divergences and nuisance parameters as well.

## Preserve unequal retained chains

A posterior draws data frame retains `.chain`, `.iteration` and `.draw`. Remove
one draw while preserving that metadata and use the same row selection for the
likelihood. Within each chain, retain iteration order.

```{r 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)
```

The pooled result and separate chain IJ calculations are available. The current
R-hat/ESS implementation requires equal chain lengths with at least four draws
each, so these entries are omitted. A missing diagnostic is not a failed test or
proof of convergence. A matrix whose metadata were lost cannot reconstruct
unequal chain identities reliably; supply correct `chain_id` explicitly.

## Interpret influence as local sensitivity

```{r 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]))
```

Exact original cluster IDs are retained in `full$unit_ids`, named by the
matrix's display labels. Distinct numeric IDs can otherwise print alike;
display names are made unique and should not be parsed back into IDs.
`full$diagnostics$chain_ids` provides the same mapping for enabled chain diagnostics.

`influence(full)` is a quantities-by-units matrix. Divide by the number of units
to obtain sample-covariance estimates of derivatives with respect to likelihood
weights. A positive derivative means that upweighting this group locally raises
the posterior expectation. Downweighting has the opposite sign.

These returned values are uncentered case-weight sensitivities multiplied by
the unit count; their row means need not be zero. IJ covariance centers across
units when constructing the sampling variance. The following reconstructs it.

```{r diag-covariance}
stopifnot(isTRUE(all.equal(
  cov(t(influence(full))) / full$n_units, vcov(full))))
```

A small weight change $\Delta w$ gives the local approximation
$\Delta T\approx C_j\Delta w$. Deleting a whole group is a finite perturbation
and may differ substantially. Large absolute sensitivity identifies a case to
inspect, rather than a reason to remove it.

## Independently check the sign and scale

Perturb a group's likelihood weight by a small amount and recompute the
normalized weighted draw mean. Subtract the maximum log weight before
exponentiating to keep the calculation stable.

```{r 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)))
```

The normalization matters when checking exact finite-draw equality. IJSE uses
sample covariance with denominator $S-1$; differentiating a normalized empirical
weighted mean uses denominator $S$. The adjustment above reconciles them. The
reweighted mean is a numerical check on these draws, not a replacement for
refitting after a substantial data change.

## Reproduce common input failures

Missing cluster IDs must be repaired using the original sampling records. Do
not silently omit affected observations after fitting the posterior.

```{r 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))
```

A used explicit NA factor level is still a missing sampling ID, even when
`anyNA()` on the factor returns FALSE. IJSE checks this case explicitly.
The literal character value `"NA"` is a different, nonmissing ID.

A one-column callback must return a matrix. R's default subsetting can drop that
dimension, so preserve it explicitly.

```{r 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))))
```

Mismatched supplied row names are rejected. Unnamed inputs require equal care;
software cannot establish their semantic correspondence from dimensions alone.

```{r 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))
```

The last error occurs because influence retention is opt-in. Recompute with
`save_influence = TRUE` before calling `influence()`.

When draw rows are unlabeled, labels supplied by different likelihood blocks
must still agree with each other. This can expose a provider that changes draw
order partway through streaming.

```{r 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))
```

The first labeled block establishes the reference across blocks. Agreement is
necessary but does not prove correspondence to unlabeled original draws.

## Check the sampling units, not only the dimensions

```{r 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)
```

Omitting cluster IDs is a valid independent-observation request to the API. It
cannot know that your scientific design required groups. IDs of the correct
length can also be in the wrong order. Validate IDs against fitting-row keys and
inspect `clustered`, `n_obs`, `n_units` and the original records explicitly.
Renaming groups does not alter the estimator, and correctly permuting likelihood
columns together with their group IDs must leave it unchanged.

```{r 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))))
```

For an unexpected result, reproduce the smallest case, inspect original extraction
errors, check draw order and units, and compare a dense calculation with a small
streamed calculation. Preserve the input-generating code and session information.
Only then investigate numerical precision or methodological calibration. The
`ijse-study` vignette gives examples where correct computation still leaves
undercoverage from few clusters or prior bias.