Infinitesimal Jackknife Inference with IJSE: A Reproducible Vignette Study

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 Giordano and Broderick (2024) and the clustered-data formulation of Ji et al. (2026). The design and reporting separate aims, data-generating mechanisms, estimands, methods, and performance measures, following Morris et al. (2019). 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.

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]
#>     Estimate Posterior_SD        IJ_SE 
#>   0.05929163   0.07112025   0.15050797
vcov(result)
#>            mu
#> mu 0.02265265
vcov(result, type = "posterior")
#>             mu
#> mu 0.005058089

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

Coverage and calibration

Coverage percentages, with Monte Carlo SE in percentage points in parentheses. Cluster IJ uses normal intervals.

Scenario	Posterior	IJ (ind.)	IJ (cl.)
Matched (50)	94.9 (0.35)	94.2 (0.37)
Matched (200)	95.0 (0.35)	94.6 (0.36)
Variance x4 (50)	66.8 (0.74)	94.5 (0.36)
Variance x4 (200)	68.4 (0.74)	95.4 (0.33)
t5 variance x4 (200)	67.3 (0.74)	95.2 (0.34)
8 equal clusters	56.4 (0.78)	53.4 (0.79)	90.8 (0.46)
30 equal clusters	56.4 (0.78)	56.2 (0.78)	95.0 (0.34)
100 equal clusters	56.6 (0.78)	56.7 (0.78)	95.2 (0.34)
30 unequal clusters	41.6 (0.78)	40.6 (0.78)	92.4 (0.42)
Strong prior (20)	17.3 (0.60)	6.0 (0.37)

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 66.8% at \(N=50\) and 68.4% at \(N=200\). The corresponding independent IJ coverages are 94.5% and 95.4%. The heavy-tailed case has independent IJ coverage of 95.2%.

Ignoring within-cluster dependence produces substantial undercoverage in all clustered designs. Cluster aggregation improves calibration, but its normal intervals cover only 90.8% with eight clusters. The t reference gives 95.3% there. In the unequal design, the normal and t-reference coverages are 92.4% and 93.1%, respectively. That remaining shortfall is evidence against treating a simple reference change as a general small-cluster solution.

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.

Bias, empirical sampling SD, and mean estimated uncertainty divided by empirical SD.

Scenario	Bias	Emp. SD	Post/SD	IJ ind/SD	IJ cl/SD
Matched (50)	0.002	0.141	1.004	0.996
Matched (200)	0.000	0.071	0.993	0.991
Variance x4 (50)	0.000	0.282	0.501	0.998
Variance x4 (200)	0.001	0.139	0.508	1.014
t5 variance x4 (200)	-0.002	0.140	0.506	1.007
8 equal clusters	-0.005	0.283	0.395	0.377	0.968
30 equal clusters	0.006	0.144	0.402	0.397	1.009
100 equal clusters	0.001	0.079	0.401	0.400	1.015
30 unequal clusters	-0.004	0.188	0.275	0.269	0.964
Strong prior (20)	-0.445	0.124	1.346	0.991

The strong-prior example isolates the distinction between precision and bias. Its empirical bias is -0.445, while independent IJ SEs are about 0.991 times the empirical SD. Nevertheless, IJ interval coverage is only 6.0%. 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.

The relative RMSE decreases from 11.5% at 200 draws to 1.0% 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.

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.

Part C: actual brms median regression

Results

All 400 datasets produced numerical results and met the stated diagnostic criterion on their first attempt. The retained fits had zero divergent transitions, maximum R-hat 1.0142, and minimum bulk ESS 403.6.

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.

Actual brms fits: empirical SD and mean SE-to-SD ratios for both regression coefficients.

Scenario	Term	Emp. SD	Post/SD	IJ ind/SD	IJ cl/SD
independent	Intercept	0.084	0.876	1.100
independent	Slope	0.085	0.855	1.043
clustered	Intercept	0.166	0.348	0.431	0.963
clustered	Slope	0.074	0.774	0.931	0.974

Actual brms coverage percentages, with Monte Carlo SE in percentage points in parentheses.

Scenario	Term	Posterior	IJ (ind.)	IJ (cl.)
independent	Intercept	90.5 (2.07)	96.5 (1.30)
independent	Slope	91.0 (2.02)	94.0 (1.68)
clustered	Intercept	51.5 (3.53)	58.0 (3.49)	93.5 (1.74)
clustered	Slope	86.5 (2.42)	91.5 (1.97)	91.5 (1.97)

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.

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 90.5% for the intercept and 91.0% for the slope, compared with 96.5% and 94.0% using independent IJ. In the clustered scenario, intercept coverage is 51.5% using posterior intervals, 58.0% using independent IJ, and 93.5% using clustered IJ. Clustered IJ slope coverage is 91.5%, compared with 91.5% 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 35.5 percentage points, with paired Monte Carlo SE 3.39 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.

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

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:

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)

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:

# 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 (Giordano and Broderick 2024). 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.