This tutorial runs without brms, Stan or a C++ compiler. You will construct matched posterior draws and log likelihoods, distinguish posterior uncertainty from sampling uncertainty, calculate a nonlinear quantity, and verify a streamed calculation. All observations are simulated.
Install the supplied source release after installing its required dependency:
install.packages("posterior")
install.packages("IJSE_0.3.3.tar.gz",
repos = NULL, type = "source")
library(IJSE)The package supplies three interfaces. ijse_draws()
accepts matched draws and log likelihoods from any backend.
ijse() extracts them from a brms fit. IJ_se()
preserves the original brms interface and returns only a named vector of
IJ standard errors. Start with ijse_draws() to see what the
calculation uses.
Generate 80 independent observations with mean 0.5 and SD 2. Fit a deliberately narrow working likelihood, \(y_i\mid\mu\sim N(\mu,1)\), and prior \(\mu\sim N(0,4)\). The posterior has variance \(v=(N+1/4)^{-1}\) and mean \(m=v\sum_i y_i\). Sampling it directly separates this first lesson from MCMC.
set.seed(431)
y <- rnorm(80, mean = 0.5, sd = 2)
n <- length(y)
v <- 1 / (n + 1 / 4)
m <- v * sum(y)
mu <- rnorm(20000, mean = m, sd = sqrt(v))
draws <- cbind(mu = mu)
rownames(draws) <- paste0("s", seq_along(mu))
ll <- vapply(y, function(yi) dnorm(yi, mu, log = TRUE),
numeric(length(mu)))
rownames(ll) <- rownames(draws)
c(draw_rows = nrow(draws), likelihood_rows = nrow(ll),
observations = ncol(ll))
#> draw_rows likelihood_rows observations
#> 20000 20000 80Rows index posterior draws; likelihood columns index observations. Column \(i\) contains \(\log p(y_i\mid\mu^{(s)})\) for every retained draw \(s\). The two inputs must use exactly the same draw rows. Row names provide a useful check, but do not prove that someone assigned correct labels to the data.
Use the likelihood, excluding the prior. A log posterior or one joint log likelihood column cannot replace the required observation contributions.
result <- ijse_draws(draws, ll, diagnostics = "none")
data.frame(Estimate = unname(coef(result)),
Posterior_SD = unname(result$posterior_sd),
IJ_SE = unname(result$se))
#> Estimate Posterior_SD IJ_SE
#> 1 0.2021575 0.1117145 0.2464216Estimate is the posterior expectation, approximated by
the draw mean. Posterior_SD is the spread of the posterior
distribution. IJ_SE estimates how the posterior expectation
varies across repeated samples of the observations. In this example the
working variance is too small, so the two uncertainty measures have a
substantive reason to differ. A larger IJ SE is a result to interpret,
rather than an automatic sign that a model was repaired.
confint(result, level = 0.95)
#> 2.5 % 97.5 %
#> mu -0.2808199 0.685135
vcov(result)
#> mu
#> mu 0.0607236
vcov(result, type = "posterior")
#> mu
#> mu 0.01248012confint() returns a normal interval centered on the
posterior expectation, using the IJ SE. It is a sampling-uncertainty
interval, rather than an equal-tailed posterior credible interval.
vcov(result) contains IJ sampling covariance;
vcov(result, type = "posterior") contains covariance of the
posterior draws. Both are matrices even with one quantity. The complete
summary(result) also reports the SE ratio and
intervals.
An IJ interval can be poorly calibrated if the estimator is biased, if the sampling units are wrong, or if the variance approximation is inaccurate. Changing a standard error does not remove those problems by itself.
For this normal model, the exact IJ SE under the package’s finite-sample convention is \(v\sqrt{N s_y^2}\).
exact_se <- v * sqrt(n * var(y))
data.frame(Analytic_IJ = exact_se,
Draw_based_IJ = unname(result$se))
#> Analytic_IJ Draw_based_IJ
#> 1 0.2460452 0.2464216The difference here comes from a finite number of posterior draws. There is also an exact identity for this particular model: estimated IJ SE divided by exact IJ SE equals the draws’ sample variance divided by \(v\). The following check tests numerical implementation, without assuming that a random draw estimate equals its population target exactly.
That identity is specific to this example. For MCMC, inspect the selected quantities and the stability of the covariance calculation; a large draw count alone does not establish precise IJ inference.
To estimate \(E(\mu^2\mid y)\), square each draw and then average. Squaring the posterior mean instead gives a different target. Select input parameters before applying a transformation, and return one row per input draw.
nonlinear <- ijse_draws(draws, ll,
transform = function(d) {
cbind(mu = d[, "mu"], mu_squared = d[, "mu"]^2)
}, diagnostics = "none")
coef(nonlinear)
#> mu mu_squared
#> 0.20215755 0.05334717
c(Exact_E_mu_squared = m^2 + v, Squared_E_mu = m^2)
#> Exact_E_mu_squared Squared_E_mu
#> 0.05355981 0.04109876
vcov(nonlinear)
#> mu mu_squared
#> mu 0.06072360 0.02472164
#> mu_squared 0.02472164 0.01006461The off-diagonal covariance matters for joint contrasts. For example, obtain the SE of \(E(\mu^2-2\mu\mid y)\) either directly from transformed draws or by the linear combination of the full IJ covariance matrix.
contrast <- ijse_draws(draws, ll,
transform = function(d) d[, "mu"]^2 - 2 * d[, "mu"],
diagnostics = "none")
a <- c(mu = -2, mu_squared = 1)
contrast_se <- sqrt(drop(t(a) %*% vcov(nonlinear) %*% a))
stopifnot(isTRUE(all.equal(
unname(contrast$se), contrast_se)))The same principle applies to products of jointly sampled parameters. Pairing rows from separate fitted models does not generally create a joint posterior.
A callback accepts observation indices and returns a draws-by-observations matrix for those indices, in the requested order. It must preserve a matrix when there is only one likelihood column.
read_likelihood <- function(ii) {
block <- vapply(y[ii],
function(yi) dnorm(yi, mu, log = TRUE),
numeric(length(mu)))
rownames(block) <- rownames(draws)
block
}
streamed <- ijse_draws(draws, read_likelihood,
n_obs = n, block_size = 7, diagnostics = "none")
stopifnot(isTRUE(all.equal(vcov(streamed), vcov(result))))This tutorial keeps the dense matrix to verify equality. In
production, a callback can generate blocks without ever allocating that
full matrix. Blocking limits likelihood storage; its runtime advantage
depends on the backend. It is not a general speed guarantee. If
subsetting a stored matrix in your callback, use
ll[, ii, drop = FALSE].
record <- list(
estimate = coef(result),
ij_covariance = vcov(result),
posterior_covariance = vcov(result, type = "posterior"),
interval = confint(result),
draws = result$n_draws, observations = result$n_obs,
units = result$n_units, clustered = result$clustered,
version = as.character(packageVersion("IJSE")))
stopifnot(record$units == n, !record$clustered)When reporting an analysis, state the posterior expectation, working
likelihood, sampling units, interval reference, and posterior simulation
settings. Keep IDs, transformation code and row-selection decisions with
the result. Continue with
vignette("ijse-brms", package = "IJSE") for fitted models
and clustering, or
vignette("ijse-diagnostics", package = "IJSE") for
diagnostics and error cases.