---
title: "Infinitesimal Jackknife Inference with IJSE: A Reproducible Vignette Study"
date: "October 21, 2026"
output:
  rmarkdown::html_vignette:
    css: ijse.css
params:
  data_dir: null
bibliography: study-references.bib
link-citations: true
fontsize: 12pt
mainfont: Times New Roman
sansfont: Times New Roman
monofont: Courier New
geometry: margin=1in
colorlinks: false
header-includes:
  - \usepackage{longtable,booktabs,array,float}
  - \usepackage{titlesec}
  - \titleformat{\section}{\normalfont\fontsize{13}{15}\bfseries}{\thesection}{0.6em}{}
  - \titleformat{\subsection}{\normalfont\normalsize\bfseries}{\thesubsection}{0.6em}{}
  - \makeatletter\renewcommand{\maketitle}{\begin{center}{\fontsize{16}{18}\selectfont\bfseries\@title\par}\vspace{0.5em}\begingroup\let\vspace\@gobble\normalsize\@date\par\endgroup\end{center}\vspace{0.5em}}\makeatother
vignette: >
  %\VignetteIndexEntry{Infinitesimal Jackknife Inference with IJSE: A Reproducible Vignette Study}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
library(IJSE)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dpi = 300,
  fig.align = "center", out.width = "100%",
  dev = if (knitr::is_latex_output()) "cairo_pdf" else "png",
  dev.args = if (knitr::is_latex_output()) list() else list(type = "cairo"))
options(knitr.kable.NA = "-")
data_dir <- params$data_dir
if (is.null(data_dir)) data_dir <- system.file("extdata", "study", package = "IJSE")
stopifnot(dir.exists(data_dir))
read_study <- function(name) read.csv(file.path(data_dir, name), stringsAsFactors = FALSE)
normal <- read_study("normal-summary.csv")
design <- read_study("normal-design.csv")
audit <- read_study("draw-audit-summary.csv")
quantile <- read_study("quantile-summary.csv")
accounting <- read_study("quantile-accounting.csv")
paired <- read_study("quantile-paired-coverage.csv")
source(system.file("study", "study-figures.R", package = "IJSE"), local = TRUE)
table_out <- function(x, caption, digits = 3) {
  knitr::kable(x, caption = caption, digits = digits, row.names = FALSE,
    format = if (knitr::is_latex_output()) "latex" else "html", booktabs = TRUE,
    longtable = FALSE, position = "H")
}
percent <- function(x, digits = 1) sprintf(paste0("%.", digits, "f"), 100 * x)
value <- function(scenario, method, column = "coverage") {
  normal[normal$scenario == scenario & normal$method == method, column]
}
qvalue <- function(scenario, parameter, method, column = "coverage") {
  quantile[quantile$scenario == scenario & quantile$parameter == parameter &
             quantile$method == method, column]
}
short_names <- c("Matched (50)", "Matched (200)", "Variance x4 (50)",
 "Variance x4 (200)", "t5 variance x4 (200)", "8 equal clusters",
 "30 equal clusters", "100 equal clusters", "30 unequal clusters", "Strong prior (20)")
```

## Purpose and scope

IJSE computes infinitesimal jackknife (IJ) standard errors for posterior
expectations using posterior draws and pointwise log likelihoods. This vignette
examines three questions: whether IJ uncertainty tracks repeated-sampling
variation under likelihood misspecification, whether cluster aggregation matters,
and how numerical variation from a finite number of posterior draws differs
from statistical miscalibration. It also demonstrates the package interface.

The study has three parts. Part A uses 4,000 replications in each of ten
conjugate-normal scenarios, where posterior expectations and IJ standard errors
are available exactly. Part B evaluates the package on genuine independent
posterior draws at four simulation sizes. Part C fits 200 actual brms median
regressions in each of two scenarios. These are simulation examples with known
targets; they are not analyses of observed educational or economic data.

The IJ construction follows @giordano2024 and the clustered-data formulation of
@ji2026. The design and reporting separate aims, data-generating mechanisms,
estimands, methods, and performance measures, following @morris2019. Numerical
equivalence and empirical coverage answer different questions, so they are
reported separately.

## What is being estimated?

Write $T(y)=E\{g(\theta)\mid y\}$ for a posterior expectation and
$\ell_i(\theta)$ for observation $i$'s log likelihood. The sensitivity to its
likelihood weight is

$$C_i=\operatorname{Cov}_{\theta\mid y}\{g(\theta),\ell_i(\theta)\}.$$

IJSE uses the finite-sample covariance convention

$$\widehat V_{\mathrm{IJ}}
=\frac{N}{N-1}\sum_{i=1}^{N}(C_i-\bar C)(C_i-\bar C)' .$$

For clustered data, first form $C_j^{\mathrm{cl}}=\sum_{i\in j}C_i$, then use
the same expression over $J$ clusters, replacing $N$ with $J$. The posterior
draw covariance and the sampling covariance of a posterior expectation are
different objects. `vcov(result)` returns the latter;
`vcov(result, type = "posterior")` returns the former.

IJSE's `confint()` method centers normal intervals on the estimated posterior
expectations. Their frequentist coverage depends on estimator bias, the
variance approximation, the reference distribution, and the sampling units.
Correcting a standard error cannot remove point-estimator bias.

For every simulation comparison, the interval labeled "Posterior" is the
posterior mean plus or minus 1.96 posterior SDs. It is a symmetric normal
approximation, rather than an equal-tailed posterior-quantile interval.
The IJ intervals use the same center and replace the posterior SD with the
applicable IJ SE; the cluster t-reference sensitivity uses its stated critical
value instead of 1.96.

## A short executable example

Consider a working model $y_i\mid\mu\sim N(\mu,1)$ and prior
$\mu\sim N(0,100)$. Its posterior variance and mean are
$v=(N+0.01)^{-1}$ and $m=v\sum_i y_i$. Here the data variance is four, creating
a deliberate mismatch with the working likelihood.

```{r short-example}
set.seed(417)
y <- rnorm(200, sd = 2)
v <- 1 / (length(y) + 0.01)
m <- v * sum(y)
mu <- rnorm(6000, m, sqrt(v))
draws <- cbind(mu = mu)
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE),
             numeric(length(mu)))
result <- ijse_draws(draws, ll, diagnostics = "none")
summary(result)$coefficients[, 1:3]
vcov(result)
vcov(result, type = "posterior")
```

The input has draws in rows and observations in likelihood columns. Their row
order must match. Matrix row names are checked when both inputs supply them;
dimension checks alone cannot verify that unnamed rows were paired correctly.

## Part A: repeated sampling with an exact posterior

### Data-generating mechanisms

The working observation variance is one in every scenario. For independent
data, errors are either normal with variance one or four, or Student $t_5$
scaled to variance four. The weak prior has mean zero and SD 10, and the
data-generating mean is zero. The strong-prior example instead has true mean
one, prior mean zero and SD 0.25, and $N=20$.

Clustered observations are generated as

$$y_{ij}=\mu_0+\sqrt{0.6}\,u_j+\sqrt{0.4}\,e_{ij},$$

where all $u_j$ and $e_{ij}$ are independent standard normals. Thus the marginal
variance is one and the intraclass correlation is 0.6. Three designs have equal
cluster size ten. The unequal design has sizes 2, 4, 8, 16, and 32, each repeated
six times. Point estimation deliberately uses the independent working likelihood;
the random effects are in the data-generating mechanism, not the fitted model.

```{r design-table, results='asis', echo=FALSE}
d <- data.frame(Scenario = short_names, N = design$n,
                J = ifelse(design$clusters == 0, "-", design$clusters),
                Variance = design$true_variance, `Prior SD` = design$prior_sd,
                check.names = FALSE)
table_out(d, "Part A design. Every row uses 4,000 independently generated datasets.")
```

### Exact calculations and evaluation criteria

The estimand for sampling-variance evaluation is the posterior mean $m(y)$.
Coverage is evaluated against the data-generating mean $\mu_0$. These targets
must be distinguished in the strong-prior case, where $m(y)$ is biased for
$\mu_0$. Under the working likelihood,

$$C_i=v(y_i-m),\qquad
\mathrm{SE}_{\mathrm{IJ}}=v\sqrt{N s_y^2}.$$

For clusters, the score contribution is
$C_j^{\mathrm{cl}}=v\{\sum_{i\in j}y_i-n_jm\}$. The cluster size $n_j$ is
retained; unequal clusters are not replaced by equally weighted cluster means.
These closed-form formulas remove posterior simulation error from Part A.
They are independently checked against the package's covariance calculation;
Part A does not claim to perform MCMC fits.

For independent errors with variance $\sigma_0^2$, the exact sampling SD of
$m$ is $v\sqrt{N\sigma_0^2}$. For the clustered mechanism it is

$$v\sqrt{0.6\sum_j n_j^2+0.4N}.$$

We report empirical bias, empirical SD, mean estimated SE, and coverage of
nominal 95% intervals. Each coverage estimate has Monte Carlo SE
$\sqrt{\hat p(1-\hat p)/B}$, where $B=4,000$. Figures show Wilson 95%
intervals for the coverage probabilities. These describe uncertainty in the
simulation result, not intervals for a model coefficient. Approximate Monte
Carlo SEs for bias, empirical SD and mean SE are in the supplied CSV tables.

The primary IJ intervals use a standard-normal reference, matching `confint()`.
A $t_{J-1}$ reference is included only as a clustered sensitivity analysis.
For equal-sized normal clusters with true and prior means both zero, the
corresponding statistic has an exact cluster-level t distribution. That result
does not justify applying the same reference universally to unequal clusters
or other likelihoods.

### Coverage and calibration

```{r normal-coverage-table, results='asis', echo=FALSE}
rows <- lapply(seq_len(nrow(design)), function(i) {
  s <- design$scenario[i]
  entry <- function(method) {
    z <- normal[normal$scenario == s & normal$method == method, ]
    if (!nrow(z)) return("-")
    paste0(percent(z$coverage), " (", percent(z$coverage_mcse, 2), ")")
  }
  data.frame(Scenario = short_names[i], Posterior = entry("Posterior"),
               `IJ (ind.)` = entry("IJ independent"), `IJ (cl.)` = entry("IJ cluster"),
               check.names = FALSE)
})
table_out(do.call(rbind, rows),
 "Coverage percentages, with Monte Carlo SE in percentage points in parentheses. Cluster IJ uses normal intervals.")
```

With matched variance, posterior intervals have coverage close to 95%.
Independent IJ intervals show some finite-sample undercoverage at $N=50$.
When the true variance is four but the working variance is one, posterior
coverage is `r percent(value("variance_50", "Posterior"))`% at $N=50$ and
`r percent(value("variance_200", "Posterior"))`% at $N=200$.
The corresponding independent IJ coverages are
`r percent(value("variance_50", "IJ independent"))`% and
`r percent(value("variance_200", "IJ independent"))`%.
The heavy-tailed case has independent IJ coverage of
`r percent(value("heavy_t5_200", "IJ independent"))`%.

Ignoring within-cluster dependence produces substantial undercoverage in all
clustered designs. Cluster aggregation improves calibration, but its normal
intervals cover only `r percent(value("clusters_8", "IJ cluster"))`% with eight
clusters. The t reference gives
`r percent(value("clusters_8", "IJ cluster: t reference"))`% there. In the unequal
design, the normal and t-reference coverages are
`r percent(value("clusters_unequal", "IJ cluster"))`% and
`r percent(value("clusters_unequal", "IJ cluster: t reference"))`%, respectively.
That remaining shortfall is evidence against treating a simple reference change
as a general small-cluster solution.

```{r normal-coverage-figure, echo=FALSE, fig.width=7, fig.height=6, fig.cap="Part A coverage. Points summarize 4,000 replications per scenario; horizontal bars are Wilson 95% intervals for simulation coverage. The dashed vertical line marks 95%. The t reference is a sensitivity analysis and is not the package default."}
normal_coverage_plot(data_dir)
```

```{r calibration-table, echo=FALSE, results='asis'}
rows <- lapply(seq_len(nrow(design)), function(i) {
  s <- design$scenario[i]
  base <- normal[normal$scenario == s & normal$method == "Posterior", ]
  data.frame(Scenario = short_names[i], Bias = base$bias, `Emp. SD` = base$empirical_sd,
    `Post/SD` = base$se_ratio, `IJ ind/SD` = value(s, "IJ independent", "se_ratio"),
    `IJ cl/SD` = if (design$clusters[i]) value(s, "IJ cluster", "se_ratio") else NA_real_,
    check.names = FALSE)
})
table_out(do.call(rbind, rows), "Bias, empirical sampling SD, and mean estimated uncertainty divided by empirical SD.")
```

The strong-prior example isolates the distinction between precision and bias.
Its empirical bias is `r sprintf("%.3f", value("strong_prior_20", "Posterior", "bias"))`,
while independent IJ SEs are about
`r sprintf("%.3f", value("strong_prior_20", "IJ independent", "se_ratio"))` times
the empirical SD. Nevertheless, IJ interval coverage is only
`r percent(value("strong_prior_20", "IJ independent"))`%. The SE is close to the
sampling SD of the biased estimator; this does not make the estimator close to
the true mean.

## Part B: numerical precision of the package calculation

For one fixed variance-misspecified dataset from Part A, we generated 80
independent ensembles at each of 200, 1,000, 5,000 and 20,000 exact normal
posterior draws. Each ensemble was passed through `ijse_draws()` with a
likelihood callback. No draws are shared across the ensembles.

For this intercept-only example, a useful exact identity is available. Because
the observation-specific constant in the covariance disappears after centering,

$$\frac{\widehat{\mathrm{SE}}_{\mathrm{IJ}}}
 {\mathrm{SE}_{\mathrm{IJ,exact}}}=\frac{s_\mu^2}{v}.$$

For IID normal posterior draws, $(S-1)s_\mu^2/v$ has a chi-squared distribution
with $S-1$ degrees of freedom. Hence the relative numerical SD is
$\sqrt{2/(S-1)}$. This identity is checked in every ensemble. It is specific
to this example; autocorrelated MCMC draws and general transformed quantities
require separate precision assessment.

```{r draw-audit-table, echo=FALSE, results='asis'}
d <- data.frame(Draws = audit$draws, Ensembles = audit$ensembles,
  `Mean error (%)` = 100 * audit$mean_relative_error,
  `RMSE (%)` = 100 * audit$relative_rmse,
  `Theory SD (%)` = 100 * audit$theoretical_relative_sd, check.names = FALSE)
table_out(d, "Numerical error relative to the exact IJ SE. Each row contains 80 independent posterior-draw ensembles.", digits = 2)
```

The relative RMSE decreases from
`r sprintf("%.1f", 100 * audit$relative_rmse[1])`% at 200 draws to
`r sprintf("%.1f", 100 * tail(audit$relative_rmse, 1))`% at 20,000 draws.
Increasing the number of draws addresses this numerical component; it cannot
repair an incorrect sampling unit, finite-cluster approximation, or prior bias.

```{r draw-precision-figure, echo=FALSE, fig.width=7, fig.height=3.8, fig.cap="Part B numerical precision. Each dot is one of 80 independent posterior-draw ensembles. Boxes show the median and interquartile range; whiskers reach the most extreme observations within 1.5 interquartile ranges of the box. Dashed curves mark the exact IID-normal 2.5th and 97.5th percentiles of relative SE error, not confidence limits for its mean."}
draw_precision_plot(data_dir)
```

## Part C: actual brms median regression

### Design and fitting protocol

Using brms [@burkner2017], we next fit a working asymmetric-Laplace median regression with estimated
scale. Both scenarios have target conditional marginal median
$0.5+0.7x$, with $x\sim N(0,1)$. For 200 independent observations, errors are
$\exp(0.3x)e$, where $e\sim N(0,1)$, inducing heteroscedasticity while retaining
zero conditional median. For clustered data, there are 30 clusters of ten
observations and errors are
$\sqrt{0.6}u_j+\sqrt{0.4}e_{ij}$. The fitted model includes an intercept and
slope but no random effects. The median is conditional on x and marginal over
the generated cluster effects. Thus cluster-robust IJ targets a marginal
working-independence fit; this is not an evaluation of a hierarchical AL model.

Each scenario has 200 generated datasets. Normal priors with SD 2.5 are used
for the slope and brms's centered intercept, and an exponential(1) prior is
used for the AL scale. The initial fit uses two chains, 1,500 iterations per
chain, 750 warmup iterations, and `adapt_delta = 0.99`. The maximum tree depth
is 12. A single prespecified retry doubles the iterations and uses
`adapt_delta = 0.995` if any divergence occurs, maximum R-hat exceeds 1.05,
minimum bulk ESS is below 200, or fitting fails. Every attempt is retained in
the audit CSV; the protocol does not inspect coverage when deciding on a retry.

Fixed quantities with zero posterior SD are excluded from convergence summaries.
Maximum R-hat and minimum ESS otherwise refer to all saved varying model
quantities. Numerical completion and the stated diagnostic criterion are
reported separately. Summaries use all numerically completed final fits, with
an additional adequate-fit sensitivity table. Unavailable fits are not silently
counted as covered or omitted from availability reporting.

```{r quantile-accounting, echo=FALSE, results='asis'}
d <- accounting[c("scenario", "datasets", "completed", "adequate", "retries")]
names(d) <- c("Scenario", "Datasets", "Completed", "Adequate", "Retried")
table_out(d, "Actual brms fit accounting. Retried datasets contribute their final fit.", digits = 0)
```

### Results

All `r sum(accounting$datasets)` datasets produced numerical results and met
the stated diagnostic criterion on their first attempt. The retained fits had
zero divergent transitions, maximum R-hat
`r sprintf("%.4f", max(accounting$maximum_retained_rhat))`, and minimum bulk ESS
`r sprintf("%.1f", min(accounting$minimum_retained_bulk_ess))`.

At nominal 95% coverage and 200 completed datasets, a coverage estimate near
95% has Monte Carlo SE about 1.54 percentage points. The Wilson intervals in
the figure make that uncertainty visible. These modest simulation sizes can
distinguish large failures of calibration, but do not establish small differences
between methods precisely.

```{r quantile-calibration, echo=FALSE, results='asis'}
keys <- unique(quantile[c("scenario", "parameter")])
rows <- lapply(seq_len(nrow(keys)), function(i) {
  key <- keys[i, ]
  z <- quantile[quantile$scenario == key$scenario & quantile$parameter == key$parameter, ]
  take <- function(method, column) {
    x <- z[z$method == method, column]
    if (!length(x)) NA_real_ else x
  }
  data.frame(Scenario = key$scenario, Term = key$parameter,
    `Emp. SD` = z$empirical_sd[1], `Post/SD` = take("Posterior", "se_ratio"),
    `IJ ind/SD` = take("IJ independent", "se_ratio"),
    `IJ cl/SD` = take("IJ cluster", "se_ratio"), check.names = FALSE)
})
table_out(do.call(rbind, rows), "Actual brms fits: empirical SD and mean SE-to-SD ratios for both regression coefficients.")
```

```{r quantile-coverage-table, echo=FALSE, results='asis'}
rows <- lapply(seq_len(nrow(keys)), function(i) {
  key <- keys[i, ]
  z <- quantile[quantile$scenario == key$scenario & quantile$parameter == key$parameter, ]
  entry <- function(method) {
    x <- z[z$method == method, ]
    if (!nrow(x)) return("-")
    paste0(percent(x$coverage), " (", percent(x$coverage_mcse, 2), ")")
  }
  data.frame(Scenario = key$scenario, Term = key$parameter,
    Posterior = entry("Posterior"), `IJ (ind.)` = entry("IJ independent"),
    `IJ (cl.)` = entry("IJ cluster"), check.names = FALSE)
})
table_out(do.call(rbind, rows), "Actual brms coverage percentages, with Monte Carlo SE in percentage points in parentheses.")
```

```{r quantile-coverage-figure, echo=FALSE, fig.width=7, fig.height=3.8, fig.cap="Part C: actual brms median-regression fits, with 200 datasets per scenario. Points show empirical coverage and bars show Wilson 95% intervals for coverage. The dashed line marks 95%. All IJ intervals use a normal reference. These error bars quantify simulation uncertainty, not MCMC uncertainty within a fit."}
quantile_coverage_plot(data_dir)
```

The slope and intercept should be examined separately: within-cluster dependence
does not have the same effect on every coefficient. The tables provide both
variance calibration and coverage, so a change in interval width is not itself
treated as proof of improved inference. The complete replicate-level estimates,
SEs, diagnostics and attempt records accompany the study.

In the independent scenario, posterior coverage is
`r percent(qvalue("independent", "Intercept", "Posterior"))`% for the intercept and
`r percent(qvalue("independent", "Slope", "Posterior"))`% for the slope, compared with
`r percent(qvalue("independent", "Intercept", "IJ independent"))`% and
`r percent(qvalue("independent", "Slope", "IJ independent"))`% using independent IJ.
In the clustered scenario, intercept coverage is
`r percent(qvalue("clustered", "Intercept", "Posterior"))`% using posterior intervals,
`r percent(qvalue("clustered", "Intercept", "IJ independent"))`% using independent IJ,
and `r percent(qvalue("clustered", "Intercept", "IJ cluster"))`% using clustered IJ.
Clustered IJ slope coverage is
`r percent(qvalue("clustered", "Slope", "IJ cluster"))`%, compared with
`r percent(qvalue("clustered", "Slope", "IJ independent"))`% using independent IJ.
The clustered IJ mean
SE is closer to the empirical SD, but that does not establish nominal coverage.

Coverage comparisons share the same datasets, so their simulation errors are
paired. For the clustered intercept, clustered IJ minus independent IJ improves
coverage by `r sprintf("%.1f", 100 * paired$coverage_difference[paired$scenario == "clustered" & paired$parameter == "Intercept" & paired$comparison == "IJ cluster minus IJ independent"])`
percentage points, with paired Monte Carlo SE
`r sprintf("%.2f", 100 * paired$paired_mcse[paired$scenario == "clustered" & paired$parameter == "Intercept" & paired$comparison == "IJ cluster minus IJ independent"])`
points. Paired differences for every planned comparison are supplied in the
CSV results; these are descriptive simulation comparisons, not claims that
all remaining coverage differences have been resolved.

## Worked brms usage and transformed quantities

The following code uses the first clustered fit saved by the simulation script.
The displayed results below were computed from that saved fit. Re-running the
study produces the file; no model is fitted while building this vignette.

```{r brms-example-code, eval=FALSE}
fit <- readRDS("study-output/quantile-clustered-example.rds")
stopifnot("group" %in% names(fit$data))
result <- ijse(fit, cluster_var = fit$data$group,
               save_influence = TRUE, diagnostics = "all")
stopifnot(result$clustered, result$n_units == 30L)
summary(result)
vcov(result)
vcov(result, type = "posterior")
case_weight_derivatives <- influence(result) / result$n_units
```

```{r worked-example-table, echo=FALSE, results='asis'}
d <- read_study("worked-example.csv")
d <- d[c("parameter", "Estimate", "Posterior_SD", "IJ_SE")]
names(d) <- c("Parameter", "Estimate", "Posterior SD", "IJ SE")
table_out(d,
 "The first clustered brms dataset: posterior expectations, posterior SDs and clustered IJ SEs.")
```

`influence(result)` retains the empirical influence values. Dividing by the
number of sampling units gives derivatives with respect to their likelihood
weights. For a small weight change $\Delta w$, the approximate mean change is
$C_j^{\mathrm{cl}}\Delta w$. Downweighting has $\Delta w<0$, so the mean change
has the opposite sign to the derivative. Deletion is a finite perturbation and
need not equal this local approximation. Here the units are clusters.

For a fixed covariate value $x=1$, compute the median predictor on every draw:

```{r transformed-example-code, eval=FALSE}
prediction <- ijse(
  fit, cluster_var = fit$data$group,
  transform = function(d) {
    cbind(median_at_x1 = d[, "b_Intercept"] + d[, "b_x"])
  }, diagnostics = "posterior"
)
coef(prediction)
confint(prediction)
```

```{r worked-transform-table, echo=FALSE, results='asis'}
d <- read_study("worked-transform.csv")
names(d) <- c("Quantity", "Estimate", "Posterior SD", "IJ SE")
table_out(d,
 "Posterior expectation of the median predictor at the fixed reference value x=1.")
```

The same approach supports jointly sampled products such as $a b$ and nonlinear
transformations. It targets $E\{g(\theta)\mid y\}$; for a product, this need not
equal $E(a\mid y)E(b\mid y)$. Pairing draws from separately fitted models does
not generally construct the required joint posterior. Uncertainty from estimating
a reference covariate distribution is not added by a fixed-reference transform.

## Computation, diagnostics and reproducibility

For a brms fit, `IJ_se(fit, cluster_var)` preserves the vector interface and now
skips diagnostics that would be discarded. `ijse(..., diagnostics = "posterior")`
computes quantity-level posterior diagnostics without additional chain-specific
IJ calculations. `diagnostics = "all"` includes those stability checks, while
`"none"` computes pooled inference only. These settings leave the pooled IJ
estimate unchanged. Unequally retained posterior data-frame chains use their
actual chain metadata. R-hat and ESS are omitted when chain lengths cannot
support the implemented equal-length calculation.

Chain agreement is useful evidence about numerical stability, but the reported
chain ranges are not calibrated IJ Monte Carlo SEs. Part B demonstrates why
sampling precision deserves attention even with a straightforward posterior.
Good diagnostics for a parameter alone do not prove precise estimation of all
parameter-likelihood covariances.

For supported univariate brms likelihoods, `block_size = 256` evaluates pointwise
likelihood blocks. A model-independent callback can provide the same memory
control through `ijse_draws()`. Clusters may span blocks. Already materialized
matrices use a separate dense covariance path. Multivariate and currently
unsupported ARMA pointwise likelihoods require ordinary extraction with
`block_size = NULL`. Streaming is a memory option; it is not a general speed
guarantee, and backend prediction objects may still be large.

The source package contains the complete study scripts:

```{r reproduce-commands, eval=FALSE}
# From the source directory, after installing IJSE 0.3.0:
dir.create("study-output", showWarnings = FALSE)
dir.create("model-cache", showWarnings = FALSE)
system2("Rscript", c(
  "tools/run-study.R", ".", "study-output", "4000", "80"))
system2("Rscript", c(
  "tools/run-quantile-study.R", ".", "study-output",
  "model-cache", "200"))
system2("Rscript", c(
  "tools/summarise-study.R", ".", "study-output", "figures"))
```

The normal study requires IJSE; figure generation additionally uses ggplot2.
The brms study requires brms, cmdstanr and a working CmdStan installation. The
package vignette reads compact precomputed CSV summaries and redraws the figures;
it does not rerun 400 MCMC fits during package checks. Replicate datasets are
regenerated from recorded scenario-specific seeds. Each brms attempt has a
recorded seed, diagnostic result, error and warning message. Session records and
the raw CSVs distinguish the study calculations from the cached integration
fixtures used for unit tests.

## Interpretation and remaining limits

This study supports a specific practical conclusion: posterior SDs can differ
substantially from the sampling variability of posterior expectations, and
properly chosen IJ sampling units can improve calibration in the examined
working-likelihood examples. It also demonstrates three distinct limitations.
Few or uneven clusters can leave normal IJ intervals undercalibrated; prior bias
can invalidate coverage despite accurate variance estimation; and finite
posterior simulation adds numerical uncertainty of its own.

The examples evaluate means under a known-variance normal working likelihood
and median-regression coefficients under an estimated-scale AL likelihood.
They do not establish performance for arbitrary quantiles, high-dimensional
latent effects, informative cluster sampling, multiway dependence, or all
transformed quantities. General high-dimensional IJ consistency is not
automatic [@giordano2024]. Further study of cluster leverage, suitable reference
distributions and covariance-specific Monte Carlo error would address those
questions; this release adds no universal small-cluster correction.

## References