Bayesian Residual Analysis

Definition

Bayesian residual analysis replaces classical point-estimate residuals (Pearson, deviance) with full posterior distributions over residual quantities. Because classical binary regression residuals have unknown sampling distributions, calibrating outlier plots is difficult; Bayesian residuals have continuous-valued posteriors that can be graphed and summarised with interpretable tail probabilities.

Key Ideas

Parametric residual: $r_i = y_i - p_i(\beta)$ where $p_i = F(x_i^T\beta)$ . Even though $y_i$ is discrete, $r_i$ is bounded on $(y_i-1, y_i)$ and its posterior is continuous. An observation is flagged if $\text{pr}(|r_i| > K \mid y)$ is large; Albert-Chib (1995) use $K = 0.75$ .
Latent data residual (probit): $\varepsilon_i = Z_i - x_i^T\beta$ , where $Z_i$ is the Albert-Chib (1993b) data-augmentation latent variable. A priori $\varepsilon_i \sim \mathcal{N}(0,1)$ , providing a natural calibration reference. The posterior (Albert-Chib 1995, eq. 3) is a truncated normal: $\propto \phi(\varepsilon_i)\cdot I(\varepsilon_i > -x_i^T\beta)$ if $y_i = 1$ , and $\propto \phi(\varepsilon_i)\cdot I(\varepsilon_i < -x_i^T\beta)$ if $y_i = 0$ . Outlying criterion: $\text{pr}(|\varepsilon_i| > K \mid y)$ versus the prior tail $2[1-\Phi(K)]$ .
Free computation: because the Albert-Chib (1993b) Gibbs sampler already draws $Z_i^{(g)}$ and $\beta^{(g)}$ at each iteration, the latent residual $\varepsilon_i^{(g)} = Z_i^{(g)} - x_i^T\beta^{(g)}$ is available at no extra cost.
Rao-Blackwellised density: $\hat\pi(\varepsilon_i \mid y) = G^{-1}\sum_g \pi(\varepsilon_i \mid y_i, \beta^{(g)})$ , obtained by averaging the closed-form truncated-normal density over Gibbs draws (Gelfand-Smith 1990). Smoother than empirical quantiles from raw simulation.

How It Works

Probit case

Run the two-block Albert-Chib (1993b) Gibbs sampler. At each iteration $g$ , record: $\varepsilon_i^{(g)} = Z_i^{(g)} - x_i^T\beta^{(g)}$ Summarise the $G$ draws to get $\hat\pi(\varepsilon_i \mid y)$ and compute $\text{pr}(|\varepsilon_i| > K \mid y) = G^{-1}\sum_g \mathbf{1}[|\varepsilon_i^{(g)}| > K]$ (or use the Rao-Blackwellised form for a smoother estimate).

Logistic link

$Z_i$ cannot be simulated directly from a logistic model, so only Rao-Blackwellised density estimates are available; these are computed by averaging over posterior draws of $\beta$ (Dellaportas-Smith 1993).

Longitudinal random effects probit

For the model $\text{pr}(y_{it} = 1 \mid b_i) = \Phi(x_{it}^T\beta + b_i)$ , the latent residual is $\varepsilon_{it} = Z_{it} - x_{it}^T\beta$ , which has prior $\mathcal{N}(0, 1 + \sigma^2)$ integrating over the random effect. The outlying criterion becomes $\text{pr}(|\varepsilon_{it}| > K(1+\sigma^2)^{1/2} \mid y)$ (Albert-Chib 1995, eq. 5).

Why It Matters

Classical Pearson and deviance residuals for binary regression have unknown finite-sample distributions, making it impossible to assign precise probabilities to observed outlier patterns. Bayesian residuals have interpretable posterior distributions with a known a priori reference ( $\mathcal{N}(0,1)$ for the probit latent residual).
The approach is computationally free for probit models estimated via data augmentation: no additional sampling is required beyond what the standard Gibbs sampler already produces.
Extensions to random effects and longitudinal data are clean because the latent variable structure is preserved.

Open Questions

For the logistic link, Rao-Blackwellised densities are less interpretable than the probit latent residual because no $Z_i$ simulation is available. A more principled Bayesian residual for logistic regression remains an open question.
The approach has not been extended to multivariate probit (correlated binary responses), where the latent $Z_i$ vector is higher-dimensional.
Formal multiple-testing correction for simultaneous outlier detection across all observations is not addressed.