Adjusted Inference for Pseudo Posterior Distributions with pseudopost

Introducing the pseudopost R package for adjusted inference with pseudo-likelihoods, implementing the Open-Faced Sandwich (OFS), Infinitesimal Jackknife (IJ), and related methods.

# remotes::install_github("marklhc/pseudopost")
library(pseudopost)

These days genAI has really become quite powerful, or persistent, which makes them increasingly capable of doing coding tasks, like converting scripts and ideas into an R package. Just a few months ago, I wouldn’t imagine I would have the time to make a package for the ideas below, but now AI agents can do it in a few days. And it was built by a local model Qwen3.8-27B, which is relatively slow compared to frontier models on the cloud.

So here’s this package, pseudopost, which implements the “Open-Faced Sandwich” (OFS) by Shaby (2014) and the “Infinitesimal Jackknife” by Giordano & Broderick, 2024 for adjusted inference with Markov Chain Monte Carlo (MCMC) draws when using a pseudo-likelihood.

To illustrate, consider the classic Besag pseudo-likelihood for the Ising model. The Ising model, which originated in statistical physics and later widely used in network analysis, describes the joint distribution of a set of binary random variables. Consider a set of \(p\) binary variables \(X = (X_1, X_2, \ldots, X_p)\), where each \(X_i\) can take values in \(\{0, 1\}\) (originally it was {-1, 1}). In the Ising model, the joint distribution of \(X\) is given by: \[ P(x_1, x_2, \ldots, x_p) = \frac{1}{Z(\boldsymbol\theta)} \exp\left(\sum_{i=1}^{p} \theta_{ii} x_i + \sum_{i < j} \theta_{ij} x_i x_j\right), \] where \(\theta_{ij} = \theta_{ji}\) represents the interaction between variables \(i\) and \(j\). The normalizing constant is \[ Z(\boldsymbol\theta) = \sum_{x_1, x_2, \ldots, x_p} \exp\left(\sum_{i=1}^{p} \theta_{ii} x_i + \sum_{i < j} \theta_{ij} x_i x_j\right), \] which is intractable when \(p\) is large. The Besag pseudo-likelihood sidesteps this by replacing the joint likelihood with a product of full conditionals. For variable \(i\), the full conditional distribution is Bernoulli: \[ P(X_i = 1 \mid \boldsymbol{x}_{-i}, \boldsymbol\theta) = \mathrm{logit}^{-1}(\eta_i), \] where \[ \eta_i = \theta_{ii} + \sum_{j \neq i} \theta_{ij} x_j. \]

Here, we consider an Ising model with a single common coupling, where \(\theta_{ij} = \beta\) for all \(i \neq j\), so that one shared coupling parameter \(\beta\) captures the strength of interaction between every pair of variables (rather than a distinct \(\theta_{ij}\) for each pair). We also denote \(\theta_{ii} = \alpha_i\) for \(i = 1, \ldots, p\), the variable-specific intercepts, giving \(p + 1\) parameters in total.

For one individual \(r\) with symptom profile \(x_r = (x_{r1}, \ldots, x_{rp})\), the Besag pseudo-likelihood factorizes over the \(p\) symptoms into a product of full conditionals: \[ \text{PLL}_r(\boldsymbol{\alpha}, \beta) = \prod_{i=1}^{p} \text{Bernoulli}\bigl(x_{ri} \;\big|\; \text{logit}^{-1}(\eta_{ri})\bigr), \] so each individual \(r\) contributes one pseudo-observation, the sum of the \(p\) conditional log-probabilities over their symptoms, \(\ell_r = \sum_{i} \texttt{bernoulli\_logit\_lpmf}(x_{ri} \mid \eta_{ri})\). The full log pseudo-likelihood is then \(\sum_{r=1}^{n} \ell_r\) (equivalently, the full pseudo-likelihood is the product of the individual contributions \(\exp(\ell_r)\)) — one pseudo-observation per individual, not one per symptom.

Because the pseudo-likelihood is not a true likelihood, the naive MCMC pseudo-posterior has asymptotic covariance \(Q^{-1}\) (one “bread” of the sandwich) rather than the correct Godambe covariance \(V = Q^{-1}PQ^{-1}\). The pseudopost package corrects this: the infinitesimal jackknife ("ij") estimates \(V\) directly from the draws without further model evaluation; the closed sandwich ("sandwich") and open-faced sandwich ("ofs") compute \(Q\) (Hessian) and \(P\) (score) explicitly and apply the appropriate draw transform.

A worked example: PTSD symptom network

Consider a network of five binary PTSD symptoms: intrusion, avoidance, negative cognitions, hyperarousal, and emotional numbing. Each variable \(X_i \in \{0, 1\}\) indicates whether symptom \(i\) is present in a given individual. The Ising model here has six parameters: five symptom-specific intercepts \(\alpha_1 \ldots \alpha_5\) (marginal propensities) and one common coupling \(\beta\) (the average tendency for symptoms to co-occur).

We simulate \(n = 250\) individuals from this model with \(\boldsymbol{\alpha} = (-1.5, -0.5, 0, -0.75, 0.25)\) and \(\beta = 0.5\) using a Gibbs sampler. Each individual’s symptom profile contributes one pseudo-observation (the sum of the five conditional log-probabilities), so the Godambe sandwich clusters scores at the individual level.

p <- 5L
n <- 250L
true_alpha <- c(-1.5, -0.5, 0, -0.75, 0.25)
true_beta  <- 0.5

simulate_ising <- function(p, alpha, beta, n, seed = 42L, sweeps = 200L) {
  set.seed(seed)
  x <- matrix(0L, n, p)
  for (r in 1:n) {
    x[r, ] <- rbinom(p, 1, 0.5)
    for (sw in 1:sweeps) {
      for (i in 1:p) {
        eta <- alpha[i] + beta * (sum(x[r, ]) - x[r, i])
        x[r, i] <- rbinom(1, 1, plogis(eta))
      }
    }
  }
  x
}

x_mat <- simulate_ising(p, true_alpha, true_beta, n)
head(x_mat)

     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    0    1    1
[2,]    1    0    0    0    1
[3,]    0    0    0    1    1
[4,]    0    1    0    0    1
[5,]    0    1    1    0    1
[6,]    1    1    1    1    1

dat  <- list(n = n, p = p, x = x_mat,
             magnitude_adj = 1.0, use_priors = 1L)
init <- list(list(alpha = rep(0, 5), beta = 0.0))
params <- c(paste0("alpha[", 1:5, "]"), "beta")

library(cmdstanr)
model <- cmdstan_model(
  system.file("test-models", "ising.stan", package = "pseudopost"),
  force_recompile = TRUE, quiet = TRUE)
fit <- model$sample(data = dat, chains = 2L, iter_warmup = 500L,
                    iter_sampling = 1000L, seed = 42L, refresh = 0,
                    show_messages = FALSE)

pp <- adjust_pseudo_posterior(
  fit, params,
  estimator = c("ij", "sandwich", "ofs"),
  bread = "hessian_mle", meat = "score",
  ci = "quantile",
  model = model, data = dat, init = init,
  verbose = FALSE)
print(pp)

Pseudo-posterior uncertainty adjustment
Methods   : ij.quant, sandwich.mle.score.quant, ofs.mle.score.quant 
CI level  : 0.9 
IJ method: proxy 
Params    : alpha[1], alpha[2], alpha[3], alpha[4], alpha[5], beta 

                   method    param     est     se      lo      hi
                 ij.quant alpha[1] -1.4879 0.2256 -1.8683 -1.1183
                 ij.quant alpha[2] -0.3491 0.1945 -0.6622 -0.0308
                 ij.quant alpha[3]  0.0128 0.1915 -0.3045  0.3266
                 ij.quant alpha[4] -0.7137 0.2087 -1.0578 -0.3753
                 ij.quant alpha[5]  0.2948 0.2047 -0.0342  0.6353
                 ij.quant     beta  0.4557 0.0624  0.3513  0.5574
 sandwich.mle.score.quant alpha[1] -1.4879 0.2495 -1.9045 -1.0770
 sandwich.mle.score.quant alpha[2] -0.3491 0.2135 -0.6909  0.0022
 sandwich.mle.score.quant alpha[3]  0.0128 0.2101 -0.3369  0.3606
 sandwich.mle.score.quant alpha[4] -0.7137 0.2212 -1.0781 -0.3515
 sandwich.mle.score.quant alpha[5]  0.2948 0.2170 -0.0572  0.6576
 sandwich.mle.score.quant     beta  0.4557 0.0695  0.3396  0.5692
      ofs.mle.score.quant alpha[1] -1.4879 0.2357 -1.8897 -1.1012
      ofs.mle.score.quant alpha[2] -0.3491 0.2018 -0.6762 -0.0212
      ofs.mle.score.quant alpha[3]  0.0128 0.1983 -0.3108  0.3368
      ofs.mle.score.quant alpha[4] -0.7137 0.2129 -1.0664 -0.3657
      ofs.mle.score.quant alpha[5]  0.2948 0.2087 -0.0393  0.6440
      ofs.mle.score.quant     beta  0.4557 0.0651  0.3487  0.5652

For a true likelihood the information equality gives \(P = Q\) and every correction is near-identity. For the pseudo-likelihood the three estimators give different SEs. The OFS applies \(\Omega = (Q^{-1}P)^{1/2}\) directly to the centered draws without calibrating to the empirical draw covariance \(C_{\text{emp}}\), so it reproduces the closed sandwich only in the “one-bread” case \(C_{\text{emp}} = Q^{-1}\). Here the bread is the negative Hessian at the MLE, whereas the MCMC draws’ empirical covariance \(C_{\text{emp}}\) is not that inverse; the one-bread OFS accordingly falls between the model-blind IJ and the closed sandwich for every parameter (and is in fact closer to each of them than the two are to one another).

Comparison with the full likelihood

For five symptoms the partition function has only \(2^5 = 32\) terms, so the exact likelihood is computable. We fit the full-likelihood model with the same data, priors, and MCMC settings and compare the resulting posterior standard deviations with the pseudo-likelihood fits.

all_configs <- as.matrix(expand.grid(rep(list(0:1), p)))
dat_full <- list(n = n, p = p, x = x_mat,
                 num_configs = 32L, all_configs = all_configs)

model_full <- cmdstan_model(
  system.file("test-models", "ising_full.stan", package = "pseudopost"),
  force_recompile = TRUE, quiet = TRUE)
fit_full <- model_full$sample(data = dat_full, chains = 2L, iter_warmup = 500L,
                              iter_sampling = 1000L, seed = 42L, refresh = 0,
                              show_messages = FALSE)

draws_pl  <- fit$draws(params, format = "draws_matrix")
se_pl     <- apply(draws_pl, 2, sd)
draws_fl  <- fit_full$draws(params, format = "draws_matrix")
se_fl     <- apply(draws_fl, 2, sd)
se_ij     <- pp$table$se[pp$table$method == "ij.quant"]
se_sand   <- pp$table$se[pp$table$method == "sandwich.mle.score.quant"]
se_ofs    <- pp$table$se[pp$table$method == "ofs.mle.score.quant"]
cmp <- data.frame(
  param      = params,
  pl_naive   = se_pl,
  ij         = se_ij,
  sandwich   = se_sand,
  ofs        = se_ofs,
  full_lik   = se_fl
)
knitr::kable(cmp, digits = 4)

	param	pl_naive	ij	sandwich	ofs	full_lik
alpha[1]	alpha[1]	0.2064	0.2256	0.2495	0.2357	0.2280
alpha[2]	alpha[2]	0.1916	0.1945	0.2135	0.2018	0.2043
alpha[3]	alpha[3]	0.1931	0.1915	0.2101	0.1983	0.2045
alpha[4]	alpha[4]	0.1975	0.2087	0.2212	0.2129	0.2107
alpha[5]	alpha[5]	0.1996	0.2047	0.2170	0.2087	0.2121
beta	beta	0.0556	0.0624	0.0695	0.0651	0.0639

The naive pseudo-posterior SEs fall below the full-likelihood SEs for every parameter — the pseudo-posterior is overconfident, most so for the coupling \(\beta\) (0.0562 against 0.0639). All three corrections inflate the SEs back toward the full-likelihood values: the closed sandwich (Hessian bread) is the most conservative and overshoots, running above the full-likelihood for every parameter, while the IJ and the one-bread OFS straddle the full-likelihood SEs and stay close to them.

The adjusted pseudo-posterior draws are most useful when using the full likelihood is slow and/or infeasible with Stan. For example, Ji et al. (2026) discussed the use of ij for Bayesian quantile regression, and they also have a package IJSE for working with objects from brms.