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.
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.
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.
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.
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.
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.
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.
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.
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("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.
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.
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.
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.