This tutorial uses two installed, synthetic brms fits. You can run all evaluated code without compiling Stan or downloading data. Install brms separately if needed. The optional refit at the end requires cmdstanr, CmdStan and a compiler.
Both examples fit a working asymmetric-Laplace median regression, \(y\sim 1+x\), with identity link, estimated scale, and fixed quantile 0.5. Data have generating conditional median \(0.5+0.7x\). The independent example has heteroscedastic normal errors. The clustered example has 30 groups of ten, with within-group correlation 0.6 in the generated errors. The fitted likelihood uses working independence and no random effects. These are worked examples, rather than an additional coverage study.
path <- system.file("extdata", "tutorial",
"median-examples.rds", package = "IJSE")
stopifnot(nzchar(path))
examples <- readRDS(path)
ind <- examples$independent
cl <- examples$clustered
stopifnot(inherits(ind$fit, "brmsfit"),
inherits(cl$fit, "brmsfit"))The bundle keeps fit, original fitting
data, and cluster separately. brms can discard
a variable unused by its model formula, so a missing
fit$data$group must not silently become your
independent-observation default. Use the saved IDs and verify the
fitting rows.
fit <- cl$fit
ids <- cl$cluster
stopifnot(length(ids) == nrow(cl$data),
identical(ids, cl$data$group),
isTRUE(all.equal(fit$data$x, cl$data$x)),
isTRUE(all.equal(fit$data$y, cl$data$y)),
!anyNA(ids), length(unique(ids)) == 30L)For your own univariate regression, select complete fitting rows explicitly before fitting and save their original row IDs and cluster IDs. A length check alone cannot establish substantive row order. The checks above apply to this ordinary regression; likelihoods with special observation orderings require separate care.
independent <- ijse(ind$fit, diagnostics = "posterior")
#> Loading required namespace: rstan
se_table(independent)
#> Estimate Posterior_SD IJ_SE
#> b_Intercept 0.462 0.086 0.134
#> b_x 0.546 0.086 0.134
confint(independent)
#> 2.5 % 97.5 %
#> b_Intercept 0.1987930 0.7259328
#> b_x 0.2823351 0.8087954
stopifnot(!independent$clustered,
independent$n_units == 200L)
stopifnot(isTRUE(all.equal(
IJ_se(ind$fit), independent$se)))By default, ijse() selects parameters whose names start
with b_. coef() returns posterior
expectations. The original IJ_se() returns only the IJ SE
vector; its pooled inference is the same and it skips discarded
diagnostics. Select exact parameter names when your transformation needs
a particular subset or a distributional parameter.
Compare two calculations on the same clustered-data fit. This changes the resampling units, not the posterior draws or point estimates.
observation_ij <- ijse(fit, diagnostics = "none")
cluster_ij <- ijse(fit, cluster_var = ids,
save_influence = TRUE, diagnostics = "all")
se_table(cluster_ij)
#> Estimate Posterior_SD IJ_SE
#> b_Intercept 0.655 0.043 0.091
#> b_x 0.672 0.043 0.047
data.frame(Observation_IJ = observation_ij$se,
Cluster_IJ = cluster_ij$se)
#> Observation_IJ Cluster_IJ
#> b_Intercept 0.04301248 0.09062993
#> b_x 0.04359916 0.04694187
stopifnot(cluster_ij$clustered,
cluster_ij$n_obs == 300L, cluster_ij$n_units == 30L,
isTRUE(all.equal(coef(observation_ij), coef(cluster_ij))))Cluster aggregation sums likelihood sensitivities within each group before forming sampling covariance. It treats groups as independent sampling units. It does not change the working likelihood, fit random effects, or automatically validate arbitrary hierarchical models. Its normal intervals also do not supply a universal correction for few or highly unequal clusters.
For this identity-link median model, the predictor at fixed \(x=1\) is the intercept plus the slope. At \(x=-1\), it is the intercept minus the slope. Compute both on every draw, with their difference, so all covariance terms are retained. These are posterior expectations of fixed-reference predictors; estimating a reference distribution would add another source of sampling uncertainty.
predictions <- ijse(fit, cluster_var = ids,
variables = c("b_Intercept", "b_x"),
transform = function(d) {
cbind(at_minus1 = d[, "b_Intercept"] - d[, "b_x"],
at_plus1 = d[, "b_Intercept"] + d[, "b_x"],
difference = 2 * d[, "b_x"])
}, diagnostics = "posterior")
se_table(predictions)
#> Estimate Posterior_SD IJ_SE
#> at_minus1 -0.017 0.062 0.108
#> at_plus1 1.327 0.060 0.096
#> difference 1.344 0.087 0.094
confint(predictions, parm = "difference", level = 0.90)
#> 5 % 95 %
#> difference 1.18941 1.49826The contrast has coefficient vector \((0,2)\) on the original intercept and slope. Verify its SE using the original full covariance matrix.
a <- c(b_Intercept = 0, b_x = 2)
v <- vcov(cluster_ij)[names(a), names(a), drop = FALSE]
contrast_se <- sqrt(drop(t(a) %*% v %*% a))
stopifnot(isTRUE(all.equal(
unname(predictions$se["difference"]), contrast_se)))Do not add marginal SEs to obtain a contrast SE. Use the covariance matrix or transform the joint draws. A contrast between predictors is an associational quantity here; the simulation and model do not make it a causal effect.
draw_ids uses post-warmup draw indices, ordered by
iteration within chain, then chain. IJSE applies the same indices to
parameter and likelihood extraction. For this example there are 1,000
retained draws in each of two chains. The following selects half of each
chain for a numerical stability comparison.
arr <- posterior::as_draws_array(fit,
variable = c("b_Intercept", "b_x"), regex = FALSE)
n_iter <- posterior::niterations(arr)
n_chain <- posterior::nchains(arr)
keep <- unlist(lapply(seq_len(n_chain), function(j) {
(j - 1L) * n_iter + seq.int(1L, n_iter, by = 2L)
}))
subset_result <- ijse(fit, cluster_var = ids,
draw_ids = keep, diagnostics = "posterior")
stopifnot(subset_result$n_draws == length(keep))
data.frame(All_draws = cluster_ij$se,
Selected_draws = subset_result$se)
#> All_draws Selected_draws
#> b_Intercept 0.09062993 0.09148666
#> b_x 0.04694187 0.04690679The estimates can differ because fewer posterior draws estimate the covariance less precisely. This comparison does not imply thinning is a precision remedy. Inspect convergence and collect sufficient effective posterior information.
streamed <- ijse(fit, cluster_var = ids,
block_size = 40, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(streamed),
vcov(cluster_ij))))This univariate likelihood supports streaming. The tested brms
implementation does not support IJSE streaming for regression ARMA
models, and multivariate streaming is rejected explicitly. Use
block_size = NULL in those cases. Covariance ARMA models
require explicit joint series likelihoods through
ijse_draws(): their native brms output contains LOO
conditional densities. Discrete truncated models also require explicit
likelihood matrices because backend lower-bound normalization is not
validated. Backend prediction objects can still consume memory when
likelihoods are blocked.
For measurement-error predictors, fit with
save_pars = brms::save_pars(latent = TRUE) to retain
matching latent draws. Without them, IJSE rejects prediction-time
resampling. The native likelihood still conditions on those latents; an
observed-data analysis needs a likelihood that integrates them.
Missing-response models can return NA likelihood columns; these require
a separately derived observed-data likelihood, rather than dropping
incomplete cases or substituting zero for NA entries
indiscriminately.
data.frame(Rhat = cluster_ij$diagnostics$rhat,
Bulk_ESS = cluster_ij$diagnostics$ess_bulk)
#> Rhat Bulk_ESS
#> b_Intercept 1.000973 1389.500
#> b_x 1.000191 1334.933
cluster_ij$diagnostics$chain_relative_range
#> b_Intercept b_x
#> 0.00655842 0.12990259These R-hat and ESS values refer to the selected quantities. Check the original fit’s diagnostics too, including varying nuisance parameters and divergent transitions. Chain-specific IJ ranges are stability checks, not calibrated Monte Carlo SEs of IJ covariance. The two tutorial fits have their generation seeds, fit seeds, fitting settings and diagnostic record in the installed data bundle’s companion files.
saveRDS(list(result = cluster_ij, cluster = ids,
data = cl$data, transform = predictions$call,
session = sessionInfo()), "my-ijse-analysis.rds")For reporting, identify the estimand, working likelihood, 30
independent clusters, normal interval reference, posterior settings and
diagnostic limitations. A single example’s interval width does not
demonstrate repeated-sampling coverage. The separate
ijse-study vignette evaluates that question.
The following code is displayed but not run while building the tutorial. It refits this example’s data and preserves the IDs separately. Using another backend, seed or toolchain can change finite-draw results.
model_data <- cl$data
model_ids <- model_data$group
refit <- brms::brm(
brms::bf(y ~ x, quantile = 0.5), data = model_data,
family = brms::asym_laplace(link_quantile = "identity"),
prior = c(
brms::set_prior("normal(0, 2.5)", class = "b"),
brms::set_prior("normal(0, 2.5)", class = "Intercept"),
brms::set_prior("exponential(1)", class = "sigma")),
backend = "cmdstanr", chains = 2, cores = 2,
iter = 2000, warmup = 1000, seed = 41302,
control = list(adapt_delta = 0.99, max_treedepth = 12))
refit_ij <- ijse(refit, cluster_var = model_ids)
stopifnot(refit_ij$clustered, refit_ij$n_units == 30L)To regenerate both installed examples, run
tools/generate-tutorial-data.R from the source release with
an output directory and model-cache directory. Refit execution is
separate from package checks and from the coverage simulation in the
study.