PaperlayerView exact revision on arXiv

Source-linked AI summary

Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC

Aki Vehtari, Andrew Gelman, Jonah Gabry

arXiv:1507.04544v5stat.COstat.ME

TL;DR

Estimating Bayesian models’ out-of-sample predictive accuracy can require unstable importance sampling or costly refitting. This paper develops PSIS-based LOO and WAIC computations from existing posterior draws, finding that PSIS-LOO is robust in finite cases and performs consistently well across examples.

  • Problem

    Approximate LOO via raw importance sampling can be noisy or unstable because importance weights may have large or infinite variance.

  • Method

    The paper computes pointwise predictive accuracy from posterior log-likelihood simulations and stabilizes LOO importance weights using Pareto-smoothed importance sampling.

  • Results

    PSIS consistently performs well by root mean squared error, while raw importance sampling, truncated importance sampling, and WAIC encounter variance or bias problems.

  • Takeaways & Limitations

    PSIS-LOO provides a fast, stable option for routine Bayesian predictive assessment, with diagnostics indicating when K-fold cross-validation should be used.

  • Takeaways & Limitations

    In hierarchical models with one data point per group, PSIS-LOO and WAIC can dramatically understate prediction accuracy.

Abstract

from arXiv · show

Leave-one-out cross-validation (LOO) and the widely applicable information criterion (WAIC) are methods for estimating pointwise out-of-sample prediction accuracy from a fitted Bayesian model using the log-likelihood evaluated at the posterior simulations of the parameter values. LOO and WAIC have various advantages over simpler estimates of predictive error such as AIC and DIC but are less used in practice because they involve additional computational steps. Here we lay out fast and stable computations for LOO and WAIC that can be performed using existing simulation draws. We introduce an efficient computation of LOO using Pareto-smoothed importance sampling (PSIS), a new procedure for regularizing importance weights. Although WAIC is asymptotically equal to LOO, we demonstrate that PSIS-LOO is more robust in the finite case with weak priors or influential observations. As a byproduct of our calculations, we also obtain approximate standard errors for estimated predictive errors and for comparing of predictive errors between two models. We implement the computations in an R package called 'loo' and demonstrate using models fit with the Bayesian inference package Stan.

1. Introduction

The section motivates Bayesian predictive-accuracy assessment for model comparison and presents PSIS-LOO as a fast, stable approximation that improves importance-sampling LOO, with diagnostics for when more intensive cross-validation is needed.

  • Motivation: Bayesian predictive accuracy supports model comparison, selection, and averaging, and can be estimated with cross-validation or information criteria from within-sample fits.The paper focuses on computations using the log-likelihood evaluated at posterior parameter simulations.
  • Problem: Exact LOO requires refitting across training sets, while ordinary importance sampling can be noisy because importance-weight variance may be large or infinite.PSIS is introduced to regularize the importance weights and obtain a more accurate, reliable approximation.
  • Background: WAIC is fully Bayesian, invariant to parametrization, and applicable to singular models, avoiding problems associated with DIC’s point-estimate basis.The introduction notes that DIC can yield negative effective-parameter estimates and is undefined for singular models.
  • Contribution: PSIS-LOO is more robust than WAIC in finite samples with weak priors or influential observations, despite WAIC’s asymptotic equality to LOO.Diagnostics indicate when approximation errors may be large and K-fold cross-validation should be used instead.

2. Estimating out-of-sample pointwise predictive accuracy using posterior simulations

This section defines pointwise expected log predictive density (elpd) as the target for evaluating out-of-sample predictive accuracy and describes approximating it with cross-validation or WAIC. It develops PSIS-LOO to stabilize importance sampling when posterior and leave-one-out distributions differ, especially with influential observations.

  • Pointwise predictive accuracy: elpd is the expected log pointwise predictive density for a new dataset, evaluated across the n data points one at a time.The true data-generating distributions are unknown, so elpd must be approximated.
  • Pointwise predictive accuracy: Cross-validation or WAIC approximates elpd, using the log score as the default utility for evaluating predictive density.Other utility or cost functions, such as classification error, could also be used.
  • Leave-one-out estimation: Raw importance sampling can estimate leave-one-out predictive densities from full-posterior draws when observations are conditionally independent.The approach becomes unstable because importance ratios may have high or infinite variance when the full posterior is narrower-tailed than the leave-one-out posterior.
  • Pareto smoothed importance sampling: PSIS improves LOO by smoothing importance weights and fitting a generalized Pareto distribution to the 20% largest importance ratios.The fitted Pareto shape parameter k provides sample-based information about the existence of importance-weight moments.
  • Pareto smoothed importance sampling: A robust model can improve importance-sampling performance when marginal and leave-one-out posteriors differ substantially because of non-robust modeling and highly influential observations.Robust modeling may reduce sensitivity to one or several highly influential observations.

3. Implementation in Stan

The authors implement LOO, WAIC, and K-fold cross-validation in R and Stan, but custom Stan programs require pointwise likelihood factors. The rstanarm interface automates these computations for supported regression models.

  • Implementation in Stan: R and Stan implementations provide a quick, convenient way to assess and compare model fits using LOO, WAIC, and K-fold cross-validation.The implementation is designed to use existing model-fitting workflows.
  • Implementation in Stan: Custom Stan programs require users to compute separate pointwise likelihood factors because Stan generally works with joint densities without data-point factorization.Stan does not automatically distinguish prior from likelihood contributions or expose each p(y_i|θ) factor.
  • Implementation in Stan: For supported regression models, rstanarm precomputes pointwise predictive measures and importance ratios, allowing loo to run without additional user programming after model fitting.The high-level interface uses standard R modeling syntax while estimating models with Stan’s algorithms.

4. Examples

The examples show that exact LOO is most robust under weak priors, while PSIS-LOO offers the best approximate balance with diagnostics. They also demonstrate that WAIC can fail at high posterior log-predictive variance, whereas problematic tail-shape estimates identify observations requiring more sampling, direct LOO, or robust model changes.

  • Diagnostics: WAIC starts to fail when a posterior log-predictive-density variance exceeds 0.4, whereas LOO approximations remain effective across larger tail-shape values before their error increases.The difference between approximation methods becomes clearer when the generalized-Pareto tail shape k exceeds 0.7.
  • Hierarchical-model simulations: Exact LOO is the best general choice for robust predictive-performance estimation because its error remains limited under weak priors.Approximate methods can have lower RMSE at medium scales, but their RMSE eventually increases as the population-distribution scale grows.
  • Hierarchical-model simulations: PSIS-LOO offers the best balance among approximations and provides diagnostics for identifying likely failures.Truncated importance-sampling LOO with heavy truncation closely matches WAIC, while shrinking exact LOO reduces variance at the cost of increased bias for larger scales.
  • Influential observations: High tail-shape estimates indicate surprising observations for which the full posterior poorly approximates the leave-one-out posterior.For one example, estimates converged slowly near k = 1; for another, k stayed above 1 and elpd contributions did not converge, motivating direct sampling or a more robust model.
  • Radon example: In the radon example, all tested methods were accurate, with small estimated tail shapes and computation taking less than 5 seconds.Two observations had k slightly above the preferred threshold of 0.7, so their elpd contributions could be computed directly and combined with the remaining PSIS-LOO estimates.

5. Standard errors and model comparison

Standard errors for LOO, WAIC, effective parameter counts, and model-comparison differences can be computed from pointwise contributions, but their interpretation is weakest for small samples and skewed distributions. PSIS-LOO comparisons provide practical uncertainty estimates, illustrated by differences favoring Model A in both radon and arsenic examples.

  • Standard-error estimation: LOO and WAIC standard errors are computed by taking the standard deviation of their n pointwise components and multiplying by √n.The same approach applies to K-fold cross-validation and effective numbers of parameters.
  • Standard-error estimation: Small samples complicate these standard errors because pointwise terms share posterior simulations and may have highly skewed distributions.The authors state that the calculations are most useful when n is large, when non-normality is less problematic.
  • Model comparison: Model-comparison uncertainty is estimated from the difference in elpdloo or elpdwaic, using a paired standard error based on the same n data points.The deviance scale can be obtained by multiplying the difference by −2; the non-parametric Bayesian bootstrap is an alternative.
  • Model comparison: These standard-error formulas may provide more useful uncertainty information than χ2 comparisons of deviance differences, which apply only to nested models and rely on Gaussian linear-model or asymptotic assumptions.The paper notes that further research is needed to evaluate the model-comparison performance of the proposed formulas.
  • Model comparison: 10.2 elpd with a standard error of 5.1 favored Model A over Model B in the radon example using PSIS-LOO.Model B was the corresponding multilevel model without the county-level uranium predictor.
  • Model comparison: 16.4 elpdloo with a standard error of 4.4 distinguished the arsenic models, with uncertainty driven by finite-sample variation rather than Monte Carlo error.This standard error does not decline as the number of Stan simulation draws increases.

6. Discussion

The discussion emphasizes practical implementation of LOO, WAIC, and K-fold cross-validation, while noting that predictive evaluation depends on modeling choices and has important finite-sample failure modes. Despite these limitations, the authors prefer PSIS-LOO and K-fold cross-validation for Bayesian model comparison over point-estimation approximations.

  • Practical contributions: The paper addresses coding the log-likelihood, computing information measures, and stabilizing weights for approximate LOO without refitting.These practicalities are considered within a Bayesian simulation environment.
  • Definitions and limitations: Predictive accuracy depends on defining both the model component called “the likelihood” and its factorization into “data points.”These choices determine potential out-of-sample replications and later expected log predictive density calculations.
  • Definitions and limitations: WAIC’s bias correction is asymptotic and only approximate for finite samples, while cross-validation relies on smaller-subset inference being close to full-data inference.The cross-validation limitation applies both to direct refitting and importance-sampling approximations for LOO.
  • Failure modes: PSIS-LOO and WAIC can dramatically understate prediction accuracy in hierarchical models with only one data point per group.LOO can also fail with weak priors and sparse data, including near-separation in logistic regression after removing one observation.
  • Practical contributions: The authors value quick out-of-sample prediction estimates for summarizing and comparing models, preferring PSIS-LOO and K-fold cross-validation to point-estimation approximations.AIC and DIC are cited as examples of popular approximations serving this purpose.

A. Implementation in Stan and R · A.1. Stan code for computing and storing the pointwise log-likelihood

This section shows how to compute and store pointwise log-likelihood values in Stan, using an arsenic example saved as logistic.stan. The individual log-likelihood terms are defined in the generated quantities block so Stan saves them for later analysis.

  • A.1. Stan code for computing and storing the pointwise log-likelihood: The implementation illustrates Stan code for computing and storing the pointwise log-likelihood.The example is presented in the context of the arsenic example from Section 5.3.
  • A.1. Stan code for computing and storing the pointwise log-likelihood: The example program is saved in the file logistic.stan.It uses the arsenic example from Section 5.3.
  • A.1. Stan code for computing and storing the pointwise log-likelihood: The Stan data block declares N, P, binary outcomes y[N], and the predictor matrix X.The outcomes are constrained to values between 0 and 1, and X has dimensions N by P.
  • A.1. Stan code for computing and storing the pointwise log-likelihood: The parameter block declares the coefficient vector b with length P.The model block then applies a normal(0,1) prior to b.
  • A.1. Stan code for computing and storing the pointwise log-likelihood: The model specifies b ~ normal(0,1) and y ~ bernoulli_logit(X*b).These statements define the prior and logistic likelihood for the example.
  • A.1. Stan code for computing and storing the pointwise log-likelihood: A vector[N] named log_lik is declared in generated quantities, where the individual log-likelihood terms are saved by Stan.The passage notes that this arrangement stores the separate terms for subsequent use.

A.2. The loo R package for LOO and WAIC

The loo R package efficiently computes PSIS-LOO and WAIC from posterior log-likelihood evaluations, returning predictive-accuracy estimates, uncertainty measures, and diagnostics for model assessment and comparison.

  • A.2. The loo R package for LOO and WAIC: The loo package provides loo() and waic() for efficiently computing PSIS-LOO and WAIC from fitted Bayesian models.Both functions use the methods described in the paper.
  • A.2. The loo R package for LOO and WAIC: The functions accept an S × n log-likelihood matrix and use simulation means and variances to compute effective parameter counts and LOO or WAIC.For large datasets, the log-likelihood can instead be computed from the data and posterior draws as needed.
  • A.2. The loo R package for LOO and WAIC: loo() returns elpd_loo, p_loo, looic, pointwise contributions, standard errors, and estimated Pareto k shape parameters, while waic() returns analogous WAIC quantities without Pareto k estimates.looic = −2 elpd_loo, placing the result on the conventional deviance or AIC scale.
  • A.2. The loo R package for LOO and WAIC: 4000 by 3020 log-likelihood matrix produced elpd_loo = -1968.3 (SE 15.6), p_loo = 3.1 (SE 0.1), and looic = 3936.6 (SE 31.2), with all Pareto k estimates below 0.5.The package checks Pareto k estimates by default and warns when estimates fall between 0.5 and 1 or exceed 1.
  • A.2. The loo R package for LOO and WAIC: The compare() function estimates differences in expected leave-one-out prediction errors or WAIC between two models, together with standard errors.The same comparison workflow applies to WAIC objects, whose returned structure is analogous except for Pareto k estimates.

A.3. Using the loo R package with rstanarm models

The section demonstrates fitting radon models with rstanarm and performing PSIS-LOO through the loo package. It also shows model comparison and automatic PSIS-LOO+ refitting for observations with problematic Pareto shape estimates.

  • Workflow: The radon example is fit with rstanarm, then evaluated using PSIS-LOO from the loo package.The workflow uses the rstanarm and loo packages together.
  • Model fitting: The first model includes floor, log_uranium, their interaction, and county-level varying intercepts and slopes.The model is fit with 4 cores, 2,000 iterations, and 4 chains.
  • Model fitting: A second model removes the county-level uranium predictor while retaining floor and county-level varying intercepts and slopes.The second model is created by updating the first model’s formula.
  • Model comparison: Fitted model objects can be passed directly to rstanarm’s loo method, and the resulting loo objects can be compared.The example computes loo(modelA), loo(modelB), and compare(looA, looB).
  • PSIS-LOO+: When estimated Pareto shape parameters ˆk are large, rstanarm can automatically refit models while leaving out problematic observations individually.The resulting direct elpd contributions are combined with PSIS-LOO results for the remaining observations.

A.4. Stan code for K-fold cross-validation

K-fold cross-validation in Stan repeatedly fits models on training partitions and predicts held-out data, with data partitioning handled externally before passing both datasets to Stan. The same structure can implement LOO, although importance sampling is more practical for large datasets.

  • Implementation strategy: K-fold cross-validation repeatedly partitions data, fitting on each training set and predicting its corresponding holdout set.The authors recommend partitioning in R, Python, or another data-processing environment.
  • Stan implementation: Stan code passes training and holdout data separately, then modifies generated quantities to compute holdout predictions.The logistic-regression arsenic example augments the data block with training and holdout observations and predictors.
  • Relation to LOO: LOO can use the same implementation with N_t = N − 1 and N_h = 1, but importance sampling is more practical for large datasets.Importance sampling reuses posterior draws from the model fit to the entire dataset.
Loading 1507.04544v5…