1What the method is for

Not how well the model fits, but how well it will land on data not yet seen.

What learning theory really wants to make small is the expected risk \(R(f)=\mathbb{E}_{(x,y)\sim P}[\ell(f(x),y)]\). But the distribution \(P\) that the data follow is unknown, so the only thing computable from the data at hand is the empirical risk \(\hat{R}(f)\). Fit on one set of data and measure on the same set, and the fit will always come out good. How to close the gap between these two quantities is the central question of Perspective 4.

Reading the seven perspectives in equations gives the generalization bound \(R(f) \le \hat{R}(f) + \mathfrak{C}(\mathcal{F},n,\delta)\) as one tool for holding that gap down from above: the route of evaluating a complexity term theoretically. Cross-validation takes a different route. It actually hides part of the data and measures on the part it hid. Not a theoretical upper bound, but a direct estimate obtained by taking the sample apart and putting it back.

Leave-one-out cross-validation is the most extreme form of this. Only one observation is hidden each time, and the fit uses everything else. It leaves as much data as possible available for training, which makes it easy to use on small samples. That the Maibara–Sekigahara study works at a scale of \(n=25\) settlements is a direct reason for choosing this method.

2Origin — as sample reuse, for prediction

Built to evaluate prediction, not fit.

Cross-validation was formalized in the 1970s by Seymour Geisser and Mervyn Stone, independently and at almost the same time. Geisser (1975) called his own procedure “predictive sample reuse”, and the name captures the character of the method well. When fresh data cannot be obtained, the sample at hand is used over and over, split between training and evaluation, in order to measure how good the prediction is.

After becoming a standard evaluation method in machine learning through the 1990s, the machinery was assembled again in the Bayesian setting. Vehtari, Gelman & Gabry (2017) set out a way of estimating leave-one-out predictive performance without regenerating the posterior draws, and introduced Pareto-smoothed importance sampling (PSIS) to stabilize the importance weights. The loo() function in brms is this implementation, and it is what the Maibara–Sekigahara study uses.

It plays two roles in the Maibara–Sekigahara study

The first is model comparison: loo() is used to ask which predicts better, the model with a spatial random effect or the model without. The second is influence analysis, looking at how far the result moves when one settlement is removed. The latter is not model selection; it is closer to the sensitivity-analysis tradition of epidemiology and ecology. The same phrase “leave one out” covers two different purposes, and the computation differs as well. Sections 4 and 5 treat them separately.

3Where it comes from mathematically

It comes down to one move: replace an expectation with a sample mean.

What cross-validation rests on is the operation of replacing the expectation one wants with an average one can compute. The expected risk is an integral over the distribution \(P\), and \(P\) is unknown. So each observation is made to play, once, the role of “data not yet seen”, and the losses there are averaged.

DEFINITION OF LEAVE-ONE-OUT CROSS-VALIDATION
\[ \widehat{R}_{\mathrm{LOO}} \;=\; \frac{1}{n}\sum_{i=1}^{n} \ell\Bigl(f^{(-i)}(x_i),\; y_i\Bigr) \]

\(f^{(-i)}\) … the model fitted with observation \(i\) removed
\(\ell\bigl(f^{(-i)}(x_i), y_i\bigr)\) … how well that model landed on the point it was not fitted to
\(n\) … the number of observations; the refitting is needed \(n\) times as well

In a Bayesian model the loss taken is the log predictive density. What is then evaluated is not a point estimate but the posterior predictive distribution, so predictive quality is measured with the uncertainty left in.

THE BAYESIAN FORM (elpd)
\[ \mathrm{elpd}_{\mathrm{LOO}} \;=\; \sum_{i=1}^{n} \log \int p\bigl(y_i \mid \theta\bigr)\, p\bigl(\theta \mid y_{-i}\bigr)\, d\theta \]

\(p(\theta \mid y_{-i})\) … the posterior with observation \(i\) removed
elpd … expected log pointwise predictive density. Larger is better
For comparing models, what is read is the difference \(\Delta\mathrm{elpd}\) and its standard error

Relation to information criteria

AIC, treated under Perspective 5 (information theory), estimates predictive quality by adding a complexity penalty to the fit. Leave-one-out cross-validation reaches the same target by actually removing points and measuring, instead of deriving the shape of the penalty theoretically. Both aim at “prediction on unseen data”; one gets there analytically, the other by operating on the sample. This is the kind of relation meant when Reading the seven perspectives in equations says that the same expression can be read from another perspective.

4The approximation, and what it costs

Implemented literally, it stops being practical for Bayesian models.

Implement the definition in section 3 as it stands and the model must be refitted as many times as there are observations. For a light fit such as least squares this is no problem, but for a Bayesian model whose posterior comes from MCMC it would mean running \(n\) separate sampling runs, which is not practical.

So the posterior draws obtained once from the full data are reused. The posterior with observation \(i\) removed can be written as a reweighting of the full-data posterior, and the weight is proportional to the reciprocal of that observation's likelihood.

IMPORTANCE WEIGHTS — CANCELLING THE REMOVED POINT BY WEIGHT
\[ w_i^{(s)} \;\propto\; \frac{1}{p\bigl(y_i \mid \theta^{(s)}\bigr)} \qquad \theta^{(s)} \sim p(\theta \mid y) \]

\(\theta^{(s)}\) … the \(s\)-th draw from the full-data posterior
The better a \(\theta\) explains observation \(i\), the smaller its weight. This approximates “the posterior that has not seen \(i\)”

The problem is that these weights can take very large values. When an observation is extreme for the model and only a small part of the posterior can explain it, the distribution of weights has a heavy right tail and a handful of draws dominate the sum. The estimate becomes unstable. PSIS fits a generalized Pareto distribution to the upper weights and smooths them, holding that tail down.

Since it is an approximation, whether it is working has to be checked

The estimated shape parameter \(\hat{k}\) of the generalized Pareto distribution serves directly as the diagnostic. The larger \(\hat{k}\), the heavier the tail and the worse the convergence of the importance sampling. The conventional guide has been that above \(\hat{k} > 0.7\) the approximation cannot be trusted for that observation (a sample-size-dependent threshold has since been proposed).

What deserves attention here is that the observation the diagnostic flags is precisely the observation with the largest influence on the model. A large \(\hat{k}\) means nothing other than that removing that single point moves the posterior a great deal. The breakdown of the approximation and the fragility of the result are two sides of the same thing.

The diagnostics this study actually produced

In the model without a spatial random effect, all 25 observations had \(\hat{k} \le 0.7\). In the model with the ICAR spatial random effect, 24 observations (96.0%) were fine and one observation fell above \(\hat{k} > 0.7\). The run log still carries the warning recommending moment matching. Since this model was not adopted, no recomputation was carried out. It stands as a case where the approximation is not guaranteed to hold everywhere.

Modelelpdp_loolooicPareto \(\hat{k}\)
Hierarchical beta (no spatial term)33.8 ± 2.64.4−67.725/25 fine
+ ICAR spatial random effect33.1 ± 2.44.8−66.324/25 fine, 1 above 0.7
Difference (with − without)−0.7 ± 0.3———

Computed from the log-likelihood matrix of 8,000 posterior draws × 25 observations. Adding the spatial term lowers elpd slightly. Since the difference stays at about twice its standard error it cannot be called a strong preference, but it is used as grounds for the claim that the data give no reason to include a spatial term.

5How it is used in this research — what unit is left out

Match the unit left out to the units of dependence.

The main analysis in this study is an aggregated model with the settlement as its unit. For each of 25 settlements, the response is the proportion of vocabulary used by the grandparent generation and no longer used by the junior-high-school generation (the loss rate), with population growth rate as the explanatory variable. Since the structure is one settlement per row, leave-one-out and “leave one settlement out” coincide in this model.

Onohara calls this LOSO (Leave-One-Settlement-Out) and refits the aggregated beta regression completely, 25 times over. This is the computation as defined, not the importance-sampling approximation of section 4. Fitting a beta regression is light enough that no approximation is needed.

  1. Estimate on the full data With all 25 settlements, estimate the standardized coefficient on population growth rate. \(b = +0.246\), \(p < .0001\).
  2. Remove one settlement at a time Refit the same model on the 24 settlements excluding settlement \(i\), and record the coefficient and the \(p\) value. Repeat 25 times.
  3. Watch how the coefficient moves The coefficient moved within \(b \in [0.115,\ 0.267]\). All signs positive; the direction is consistent.
  4. Count how often significance survives Of the 25 runs, 24 kept \(p<.05\). Only one broke.
  5. Check robustness to the specification as well A version placing school district as a fixed rather than random effect gives \(b = +0.228\), \(p<.0001\), unchanged.

Significance depended on the single least reliable settlement

The run that broke is the one removing Sugisawa. This settlement has only one respondent in the grandparent generation, the highest population growth rate of the 25 at 1.315, and a high loss rate of 0.5. Remove it and the coefficient becomes \(b = 0.115\), \(p = 0.2501\); significance is lost.

So while the positive association itself is consistent across 24 settlements, statistical significance is held up by the settlement with the lowest estimation precision, the one with a single respondent. The main result is driven by a point that is both high-leverage and low-confidence. LOSO showed this as a concrete number.

Why the settlement is the unit — the contrast with the count-level model

In the aggregated model one settlement is one observation, so there is no room to choose a unit. The question arises for the count-level model the same study runs alongside it, which models 5,100 cells of vocabulary item × settlement × generation directly. Naively removing a single cell there lets information leak in from the other cells of the same settlement, because cells within a settlement are dependent on one another through a shared random effect. Hide one cell and the model still knows that settlement almost as well as before.

To estimate predictive performance correctly in this situation, the unit removed has to match the units of dependence. Remove a whole settlement and the model genuinely knows nothing about it. If the question is “how well can a new settlement be predicted”, then the unit removed must be the settlement. The choice of unit in leave-one-out cross-validation is not a matter of computational convenience but a statement of which prediction task is being evaluated.

One of three directions of robustness checking

LOSO is not used on its own; it is one corner of a design that checks the hazards of analysing aggregated values from three directions. The other two are a count-level model that leaves the raw counts unaggregated, and a reanalysis that merges neighbouring settlements to increase the sample. The three results do not point the same way.

AnalysisUnitPopulation growth → lossReading
Aggregated beta regression (main)25 settlements+0.23Credible interval excludes 0
LOSO (influence analysis)24 settlements × 25 runs0.115 – 0.267Direction consistent, significance 24/25
Count-level model5,100 cells+0.17Interval straddles 0
Merged settlements (5+ grandparents)12 units−0.07Association disappears

As the unit shifts from aggregate to raw data to merged units, the claim weakens. That gradient is itself a result, and it is reported in the paper as grounds for a provisional reading.

6Limits — what the method cannot say

A tool for measuring predictive quality, not a guarantee of correctness.

A coefficient that does not move does not mean the model is right

Even if the coefficient is stable under LOSO, all that shows is that the same answer comes out whichever single point is removed. If the form of the model is wrong, a wrong answer will come out stably. What cross-validation measures is predictive performance, not causal structure and not the adequacy of the model.

A small difference in elpd is hard to interpret

The difference between the models with and without the spatial term in this study is \(-0.7 \pm 0.3\): a little over twice the standard error, not a decisive difference. Not being able to say “better” outright when the difference is small is in the character of this method, and indeed the study keeps to the weaker conclusion that the spatial term does not improve fit.

Change the unit removed and the question being evaluated changes too

Removing one cell and removing one settlement will of course give different numbers for predictive performance. It is not that one is correct; the question is different. Whether the task is “predict another vocabulary item in the same settlement” or “predict a new settlement” changes which unit should be evaluated. Numbers cannot be compared without making that choice explicit.

On small samples the diagnostic itself becomes unstable

The aggregated model in this study has \(n=25\). At that scale there is a limit to the precision of both the Pareto \(\hat{k}\) estimate and the standard error of elpd. Good diagnostics show that the approximation has not broken down; they guarantee nothing beyond that.

7Sources and data

The analysis rests on data from a three-generation dialect vocabulary survey covering 25 settlements from Maibara in Shiga Prefecture to Sekigahara in Gifu Prefecture. Bayesian estimation of the aggregated model uses brms and Stan; the LOSO refitting uses betareg. The loo() output for the model comparison and all 25 rows of the LOSO result are stored as run logs and CSV.

References for the method

  • Geisser, S. (1975) The predictive sample reuse method with applications. Journal of the American Statistical Association, 70(350), 320–328. doi:10.1080/01621459.1975.10479865
  • Vehtari, A., Gelman, A. & Gabry, J. (2017) Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27, 1413–1432. doi:10.1007/s11222-016-9696-4
  • Vehtari, A., Simpson, D., Gelman, A., Yao, Y. & Gabry, J. (2024) Pareto smoothed importance sampling. Journal of Machine Learning Research, 25. jmlr.org (first appeared as arXiv:1507.02646, 2015)

Which perspective this method belongs to

Leave-one-out cross-validation does not ask whether a probability model is right; it looks only at the loss on data not yet seen. It is a tool for estimating, directly and by taking the sample apart, the gap between the empirical risk that can be computed at hand and the expected risk that is actually wanted. Where a generalization bound holds that gap down from above with a complexity term, this closes it by measurement. It falls under Perspective 4 (learning theory, risk minimization) in Reading the seven perspectives in equations. The same target — prediction on unseen data — seen from the side of information criteria is the AIC of Perspective 5.

This page sets out, for readers in the same field, a method Onohara uses in actual research. Every figure is quoted from the output of the analysis scripts; no estimates or approximations are given.
The whole layout of where the methods sit is in Reading the seven perspectives in equations; which teaching session falls under which perspective is in Methods & Mathematics.