Tutorial 3: Diagnostics, Influence and Troubleshooting

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.

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

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.

full$diagnostics$chain_se
#>          mu
#> A 0.2291499
#> B 0.2221072
full$diagnostics$chain_relative_range
#>         mu 
#> 0.03122537
full$diagnostics$chain_ids
#>   A   B 
#> "A" "B"
c(Rhat = unname(full$diagnostics$rhat),
  Bulk_ESS = unname(full$diagnostics$ess_bulk),
  Tail_ESS = unname(full$diagnostics$ess_tail))
#>        Rhat    Bulk_ESS    Tail_ESS 
#>    1.002042 2074.844953 2044.291017

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.

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

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]))
#>      Cluster   Derivative 
#>         "g4" "0.09505291"

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.

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.

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)
#>          Covariance Adjusted_difference 
#>          0.09505291          0.09505291
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.

bad_ids <- group
bad_ids[3] <- NA_character_
show_error(ijse_draws(draws, ll, cluster_var = bad_ids))
#> [1] "'cluster_var' must be a non-missing atomic vector of length 100."
factor_ids <- factor(bad_ids, exclude = NULL)
c(AnyNA_factor = anyNA(factor_ids),
  AnyNA_values = anyNA(as.character(factor_ids)))
#> AnyNA_factor AnyNA_values 
#>        FALSE         TRUE
show_error(ijse_draws(draws, ll, cluster_var = factor_ids))
#> [1] "'cluster_var' must be a non-missing atomic vector of length 100."

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.

show_error(ijse_draws(draws, function(ii) ll[, ii],
  n_obs = length(y), block_size = 1))
#> [1] "'log_lik block' must be a numeric matrix with at least 2 rows and one column."
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.

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))
#> [1] "Row names of 'draws' and 'log_lik' must agree in order, including across likelihood blocks."
show_error(ijse_draws(draws, ll, variables = "mu_typo"))
#> [1] "'variables' must contain unique, existing parameter names."
show_error(influence(pooled))
#> [1] "Influence was not retained. Recompute with 'save_influence = TRUE'."

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.

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))
#> [1] "Row names of 'draws' and 'log_lik' must agree in order, including across likelihood blocks."

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

observation_result <- ijse_draws(draws, ll,
  diagnostics = "none")
c(Observation_units = observation_result$n_units,
  Group_units = full$n_units)
#> Observation_units       Group_units 
#>               100                20
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.

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.