IJSE

Infinitesimal jackknife (IJ) sampling standard errors and covariance matrices for posterior expectations, using one set of posterior draws. Supports independent or clustered data, brms fits, transformed quantities, and streamed likelihood blocks.

Install

For this development release, install the supplied source package:

install.packages("posterior")
install.packages("IJSE_0.3.4.tar.gz", repos = NULL, type = "source")

Install brms separately if you use the fitted-model interface:

install.packages("brms")

The matrix interface requires no Stan compiler and does not load brms.

Existing brms code still works

library(IJSE)
IJ_se(fit)                               # named vector, same finite-sample formula
IJ_se(fit, cluster_var = school_id)       # IDs must match fitted likelihood columns

result <- ijse(fit)
summary(result)                          # estimate, posterior SD, IJ SE, intervals
coef(result)
vcov(result)                            # complete IJ covariance matrix
vcov(result, type = "posterior")        # covariance of posterior draws
confint(result, parm = "b_x", level = 0.95)

ijse() defaults to fixed-effect parameters beginning with b_. Select other parameters explicitly. Posterior expectations and IJ normal intervals have different meanings from posterior SDs and credible intervals.

Use diagnostics = "posterior" for quantity-level R-hat/ESS without additional chain-specific IJ calculations, or "none" for pooled inference only. The default "all" also reports chain IJ stability. IJ_se() skips diagnostics that would be discarded. These settings leave the pooled estimate unchanged.

Transformed quantities

result <- ijse(
  fit, variables = c("b_a", "b_b"),
  transform = function(d) cbind(indirect = d[, "b_a"] * d[, "b_b"])
)

The target is E[a*b | data], computed from the product on each draw. It is not the product of posterior means. For joint quantities, supply draws from the same joint posterior; independent model fits cannot simply be paired by row number. Fixed-reference predictions and contrasts can be passed in the same way.

Model-independent matrix interface

set.seed(12)
y <- c(-1, 0, 0.5, 1, 2)
v <- 1 / (length(y) + 1 / 4)
mu <- rnorm(2000, v * sum(y), sqrt(v))
draws <- cbind(mu = mu)
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE), numeric(length(mu)))
result <- ijse_draws(draws, ll)
summary(result)

Draws are rows, observations are likelihood columns. Match draw rows exactly, including chain order and any subsetting. Row names are checked when provided by both matrices. chain_id enables chain-level diagnostics; it is inferred for posterior draws objects and brms fits.

Control likelihood memory

result <- ijse(fit, block_size = 256)

# A callback from any backend must return a numeric draws-by-length(i) matrix.
result <- ijse_draws(
  draws, log_lik = function(i) evaluate_loglik(i),
  n_obs = number_of_observations, block_size = 256
)

The callback is evaluated once per observation block. Covariance is accumulated without constructing an entire likelihood matrix. Clusters may span blocks: their covariance contributions are aggregated before the sampling variance is computed. save_influence = TRUE retains the parameters-by-units influence matrix when needed for investigation; it is off by default.

When retained, use influence(result) to obtain that matrix. Dividing by the number of sampling units gives case-weight derivatives of the posterior means.

Streaming requires brms pointwise-likelihood support. Single-column matrix outputs from categorical models are supported. Regression ARMA likelihoods in the tested brms release and multivariate streaming require block_size = NULL. Covariance ARMA models are rejected: brms returns LOO conditional densities rather than joint likelihood factors. Use ijse_draws() with one joint density per independent series. Discrete truncated models also require explicit matrix likelihoods because backend lower-bound normalization is not validated.

Measurement-error predictors require save_pars = brms::save_pars(latent = TRUE) when fitting. With unsaved latents, brms may generate new prediction values; these cannot be paired with fitted posterior draws. Saved latents give a conditional response likelihood, not the observed-data joint likelihood. Missing-response likelihoods containing NA are rejected. Explicitly integrated observed-data likelihoods can be supplied to ijse_draws() for either case. Backend prediction objects may still need substantial memory when streaming.

Diagnostics and interpretation

Missing cluster IDs (including used explicit NA factor levels), fewer than two sampling units, non-finite values, invalid block shapes, and mismatched draw counts produce explicit errors. Original likelihood errors are preserved. Cluster IDs refer to the rows actually used for fitting; matching length alone cannot establish their substantive ordering.

result$diagnostics contains chain-specific IJ SEs, relative ranges, and, when chains have equal lengths of at least four, quantity-level R-hat and bulk/tail ESS. These assess stability. The relative range is not an IJ Monte Carlo standard error, and good parameter ESS does not guarantee precise IJ covariance.

IJ describes the sampling variability of posterior expectations. It does not remove estimator bias. Correct factorization, choice of sampling unit, and model regularity remain essential. Cluster aggregation supplies neither multiway clustering nor a small-cluster correction, and arbitrary hierarchical models require separate methodological validation.

Development and verification

testthat::test_local("/path/to/IJSE")
roxygen2::roxygenise("/path/to/IJSE")

The fast suite includes analytic targets, finite-weight derivatives, legacy parity, transformations, covariance properties, and real cached brms fits. Integration fixtures are synthetic, not user data. They do not require Stan compilation during tests. Regenerate them with tools/generate-fixtures.R using brms, cmdstanr and an installed CmdStan toolchain; the generator and fixture version record are included. tools/benchmark.R reproduces the likelihood-memory comparison in fresh R processes; the benchmark excludes brms extraction and MCMC.

References

Ji, F., Lee, J., and Rabe-Hesketh, S. (2026). Valid Standard Errors for Bayesian Quantile Regression With Clustered and Independent Data. Journal of Educational and Behavioral Statistics, 51(5). doi:10.3102/10769986251379738.

Giordano, R., and Broderick, T. (2024). The Bayesian Infinitesimal Jackknife for Variance. arXiv:2305.06466.

Use citation("IJSE") for package and method citation entries.

Vignette study

vignette("ijse-study", package = "IJSE") reports a repeated-sampling study covering variance misspecification, clustered data, a strong-prior boundary, posterior-draw numerical error, and actual brms median regressions. Coverage is reported with Monte Carlo uncertainty, including undercoverage and fitting diagnostics. Complete generators are in tools/ and inst/study/; compact summary data are in inst/extdata/study/. Building the vignette reads those summaries and redraws figures, without rerunning the MCMC simulations.

Executable tutorials

Start with vignette("ijse", package = "IJSE") for a complete calculation without Stan. Continue with vignette("ijse-brms", package = "IJSE") for actual fitted median regressions, clustering and joint contrasts. Use vignette("ijse-diagnostics", package = "IJSE") for chain metadata, saved influence and deliberately reproduced input errors. All evaluated code can run without compiling Stan. The brms tutorial requires brms; its two synthetic fits are installed with separate fitting data and sampling IDs.

The same executable code is installed as scripts:

source(system.file("tutorials", "ijse.R", package = "IJSE"))
source(system.file("tutorials", "ijse-brms.R", package = "IJSE"))
source(system.file("tutorials", "ijse-diagnostics.R", package = "IJSE"))

tools/prepare-tutorials.R extracts scripts from the R Markdown sources. tools/verify-tutorials.R executes installed scripts in fresh R processes. tools/render-tutorials.R creates separate HTML lessons and a collected PDF. tools/generate-tutorial-data.R refits the synthetic examples when CmdStan is available. tools/audit-numerics.R runs an independent algebra/invariance audit; it is distinct from the calibration study and does not refit its 400 models.

Reference manual and identity checks

help("IJSE") provides a package overview. ?ijse_draws, ?ijse, and ?summary.ijse document dimensions, covariance types, influence normalization, return components, diagnostic omissions, ID mappings and executable examples. Names on vector transforms are draw identities and must preserve alignment. The exact original cluster IDs are retained in result$unit_ids, with unique print labels as names. Exact chain IDs are in result$diagnostics$chain_ids when diagnostics are enabled. Do not reconstruct numeric IDs from display text.

tools/render-reference.R generates HTML help and a reference PDF directly from the actual Rd sources. It requires R’s tools, tinytex, a working XeLaTeX setup, and the selected text fonts. No fitting or simulation is run by this renderer.

Deep verification tools

tools/audit-deep.R checks scale-sensitive covariance and influence against explicit centered cross-products, together with joint nonlinear empirical reweighting derivatives. tools/audit-inputs.R checks ID classes, diagnostic modes and posterior object formats. tools/audit-mutations.R copies the source to a separate scratch directory and verifies that a clean control passes while nine prespecified formula/metadata faults fail the relevant tests. Each copy is tested in a fresh R process; the source is not edited. This development audit requires testthat and pkgload and is separate from an ordinary package check.

tools/generate-deep-fixtures.R regenerates two additional synthetic brms integrations: a weighted Bernoulli model with omitted fitting rows and a nonlinear Gaussian model with a predictor-dependent scale. Tests compare their log likelihoods with direct probability densities and exercise actual row IDs, parameter selection, joint transformations and streaming. These fixtures assess implementation and do not constitute a repeated-sampling coverage study.

Labeled callback blocks must agree with each other as well as with available draw labels. Unlabeled rows or consistently incorrect labels still require independent validation against the original fitting records.