Joint leave-group-out cross-validation in Bayesian spatial models

Aki Vehtari; Alex Cooper; Catherine Forbes

REVIEW 2 major objections 6 minor 3 references

Joint leave-group-out cross-validation in Bayesian spatial models

T0 review · 2 major / 6 minor · reviewed 2026-08-16 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Scoring left-out spatial blocks jointly, rather than pointwise, makes cross-validation model selection more accurate when spatial dependence is strong.

desk verdict Useful, honest extension of joint CV scoring to spatial models, but the simulation evidence leans on a Laplace approximation that gets only thin validation. read the letter →

arxiv 2504.15586 v1 pith:N2LAHTVV submitted 2025-04-22 stat.ME

Alex Cooper , Aki Vehtari , Catherine Forbes This is my paper · ORCID

classification stat.ME MSC 62M3062F1562-08

keywords spatialcross-validationjointscoringleave-group-outBayesianmodelselectionsimultaneousautoregressiveGaussianMarkovrandomfieldblockpredictiverules

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper argues that cross-validation model selection for spatial data should score each left-out block as one multivariate prediction rather than as a collection of pointwise predictions. Simulations with Gaussian spatial models show that joint scoring reduces the variability of the cross-validation selection statistic and improves the probability of selecting the correct model under strong spatial dependence, roughly when standardized spatial dependence exceeds $3/4$. The finding matters because spatial cross-validation is a common tool for choosing among models, yet the scoring rule has received less attention than the blocking design.

What carries the argument

The central object is the joint predictive density $p(y_{\text{test}_k} \mid y_{\text{train}_k}, M)$ for each left-out block, compared with the pointwise aggregation $\sum_i \log p(y_{\text{test}_{k,i}} \mid y_{\text{train}_k}, M)$. The paper summarizes model-selection performance through the Z ratio, the mean of the cross-validation selection statistic divided by its standard deviation across independent data draws; joint scoring raises this ratio by reducing variance.

What would settle it

Re-run the pointwise-versus-joint comparison on a strongly dependent spatial process outside the Gaussian SAR family with row-standardized adjacency, such as a CAR or BYM2 model on an irregular geography: if joint scoring does not reduce the cross-validation statistic's variance or improve correct-selection rates, the claimed effect is specific to the SAR framing rather than spatial dependence itself.

Watch

Extended reading notes

Core claim

The paper claims that jointly evaluated scoring rules outperform pointwise aggregation in block-wise spatial cross-validation. For Gaussian SAR and covariance-kernel models on a regular lattice, evaluating the joint predictive density of each left-out block yields higher correct-selection rates and higher ratios of mean to standard deviation of the selection statistic than summing pointwise predictive densities, once standardized spatial dependence is above $3/4$. The improvement comes primarily from lower variability of the selection statistic rather than from a large shift in its mean.

Load-bearing premise

The evidence rests on Gaussian spatial models with a known row-standardized adjacency structure, and assumes the halo-and-block design keeps training and test sets nearly independent so the Gaussian approximation to the selection statistic is valid.

Editorial extensions

If this is right

Practitioners using block cross-validation with a dozen or so folds and joint scoring can expect more reliable model selection when spatial dependence is strong.
Pointwise evaluation can lose accuracy as test blocks grow under strong spatial dependence, while joint evaluation does not show that decline.
When spatial dependence is weak, the choice of scoring method and test-block size matters little for selection accuracy.
The benefit of joint scoring is largest in marginal cases where candidate models differ only subtly in predictive performance.
Real-data examples show joint and pointwise scoring select the same preferred model, but joint scoring yields larger Z ratios when the fitted spatial dependence is strong.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

If the mechanism is not specific to the Gaussian SAR family, joint scoring could also improve spatial model selection for CAR, BYM, and Gaussian-process models whenever the joint predictive distribution is computable.
Because the paper only studies the logarithmic score, other proper multivariate scoring rules such as the energy score or variogram score could show similar or different patterns; this is a direct testable extension.
The threshold of standardized dependence $\rho > 3/4$ is framed in the row-standardized SAR parameterization, so transferring the threshold to other model classes would require a model-independent measure of spatial dependence.
Block cross-validation with joint scoring may serve as a computationally cheaper alternative to leave-one-out CV, since fewer model fits are needed while selection power improves.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

The take-home: this is a solid, modest extension of Cooper et al. (2024) from time series to spatial block CV, and the main recommendation—use joint scoring when spatial dependence is strong—is plausible and supported by several independent experiments. It is not a breakthrough, but it is a genuinely useful workflow result for applied spatial statisticians.

What is new and good: the paper systematically compares pointwise versus joint evaluation of the log score in leave-group-out CV for SAR covariate selection, network structure selection, and covariance kernel selection. The simulation design is clear, with varying test-set sizes and standardized rho values, and the Z-ratio summary is a sensible way to separate location from variability. The case studies with real data (child vaccination, lung cancer) are a nice addition; they use INLA with actual posterior draws rather than the Laplace shortcut, and they agree with the simulations. The authors are also honest about limitations: no iid comparisons, only the log score, and the difficulty of generalizing beyond standardized SAR rho.

The soft spots: the main one is the Laplace approximation. All of the simulation-based CV scores are computed via the Laplace approximation described in Appendix D, but the only validation of that approximation is for two kernel-model marginal posteriors on a single data draw. There is no check for the SAR experiments that deliver the headline rho > 0.75 result, and no check on the held-out predictive densities p(ytest | ytrain) that actually feed the joint scores. Since joint scores use the full predictive covariance while pointwise scores use only marginals, an approximation error in the covariance could affect the two scoring rules asymmetrically—and that asymmetry could in principle create the very effect the paper reports. I am not saying that is what happened; the case studies, which do not use Laplace, point the same direction. But the simulation evidence would be much more convincing with a few MCMC checks for the SAR predictive densities, especially at high rho and larger test blocks. That is a real gap, not a technicality.

Minor points: the covariance kernel experiments (Section 3.2) are noisy and use trimmed means without formal statistical comparison, but the SAR experiments are clean and carry the main argument. No code is shipped, which makes the results harder to reproduce and check.

The paper deserves serious peer review. I would recommend conditional acceptance after the authors add at least a spot-check of Laplace against MCMC for the SAR setting and for the joint predictive densities, and ideally make the simulation code available. If those checks confirm the simulations, this will be a cite-worthy contribution to the spatial CV workflow literature.

Referee Report

2 major / 6 minor

Summary. The paper studies block-wise cross-validation (CV) for Bayesian spatial models with Gaussian covariance structures, focusing on whether jointly evaluating the predictive performance of left-out blocks (via the joint log predictive density) is preferable to aggregating pointwise predictive scores. Through simulation experiments on a regular lattice, covering SAR-based covariate selection, network structure selection, and covariance kernel selection, the paper reports that joint scoring yields higher model selection accuracy and lower relative variability (measured by a Z ratio) than pointwise scoring, particularly when spatial dependence is strong (standardized rho > 0.75). Two real-data case studies (Australian child vaccination rates and Pennsylvania lung cancer incidence) illustrate the approach and produce results broadly consistent with the simulations. The paper acknowledges the limitations of using the Laplace approximation for computational feasibility and discusses the limited generalizability of the standardized-rho interpretation.

Significance. The paper addresses an under-explored aspect of spatial cross-validation: the scoring rule used to evaluate left-out blocks, as opposed to the more studied blocking design. If the empirical claim is reliable, the recommendation to use joint evaluation for spatially clustered CV would be useful for practitioners, as it extends findings from time-series settings (Cooper et al., 2024) to spatial models and to model components beyond covariates, such as network structure and covariance kernels. The simulation design is a strength: common data draws are used across test set sizes, multiple model selection tasks are considered, and the paper is transparent about the limitations of its simulation-based approach. However, the central claim currently rests on simulation results whose computational approximation is only weakly validated and whose statistical uncertainty is not quantified, so the strength of the conclusions exceeds what the evidence presently supports.

major comments (2)

[Section 3, Appendix D (Eqs. D.1–D.2, Figures D.1–D.2)] All simulation-based CV scores in Section 3 are computed with the Laplace approximation described in Appendix D, but the only validation against exact MCMC inference is for marginal posterior distributions of hyperparameters in the two covariance-kernel models (Section 3.2), for a single data draw. The SAR experiments (Sections 3.1.1 and 3.1.2) that support the headline rho > 0.75 claim are not validated, nor are the held-out predictive densities p(ytest_k | ytrain_k) that enter the joint and pointwise scores in Eqs. (2) and (3). This is a load-bearing gap because joint scores use the full predictive covariance, while pointwise scores use only marginal variances; if the Laplace approximation becomes more biased as rho and block size grow, it could differentially affect joint and pointwise scores and create a spurious joint advantage. The paper itself states that the Laplace approximation 'can be a crude posterior approximation' (Section 3, page 8). Please add a validation study comparing Laplace-based CV score contributions (both joint and pointwise) with MCMC- or exact-based contributions for a subset of SAR settings, especially at high rho and large test block sizes, and report the magnitude of any discrepancies.
[Section 3, Figures 2–4] The central empirical claim that joint scoring improves selection accuracy and reduces variability is supported only by visual comparison of point estimates, without confidence intervals or formal statistical tests. For example, in the covariate selection experiment (Section 3.1.1, Figure 2), correctness percentages over 1,000 replications are plotted with no error bars, making it difficult to judge whether the joint-vs-pointwise differences are beyond sampling noise, particularly at moderate rho values where the differences appear small. In the kernel experiments (Section 3.2, Figure 4), the paper reports 98% trimmed means and variances but again no uncertainty quantification. Please provide standard errors or bootstrap confidence intervals for the correctness percentages and Z-ratio estimates, or conduct a paired comparison across the common data draws, so the reader can assess the size and precision of the reported effects.

minor comments (6)

[Section 3.1.1, Eq. (6) and surrounding text] The text below Eq. (6) refers to 'the sample standard error and mean' as defining the Z ratio, but the denominator of Eq. (6) is the sample standard deviation, not the standard error; please correct the terminology.
[Figure 4 caption] The caption's panel references appear to be mismatched: the text refers to 'Panel (1a)' and 'Panel (1b)' for both correctness and Z-ratio panels, and the parenthetical descriptions do not match the figure layout; please revise for clarity.
[Section 3.2 and Section 5] The conclusion states that leave-group-out CV with 'a dozen or few dozen folds' performs well, but the simulations include test set sizes up to ntest = 144, which corresponds to only 4 folds; please qualify the claim or exclude the very large block sizes.
[Abstract, Section 3.1.2, and Section 5] The threshold 'standardized rho > 3/4' for a clear joint advantage appears stronger than what the simulations show, since joint evaluation also appears to help for rho = 0.5 in the network selection experiment (Figure 3) and in some kernel settings (Figure 4); please align the wording with the full range of results.
[Section 4, Eq. (10)] The definition of \hat Z_y uses within-fold variance from a single dataset, whereas the Z ratio in Section 3 is computed across independent data draws; this distinction is not explained, so please clarify the relationship between the two quantities.
[Section 5 and Appendix D] The paper would benefit from a statement on code and data availability, since the simulation study is computationally intensive and reproducibility is otherwise difficult to assess; also, there is a minor typo in Section 5 ('invovling' should be 'involving').

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the simulation results are generated from independent DGP draws, and the Z ratio is a descriptive summary of the simulated distribution rather than a fitted parameter or a consequence of the defining equations.

full rationale

The paper's central claim—that jointly evaluated scoring rules reduce variability and improve model selection for spatial block CV, especially under strong spatial dependence—is supported by explicit simulation experiments in Section 3, not by a derivation from the defining equations. The CV estimators in Eqs. (2) and (3) are defined independently, and the comparison between them is made on independently generated data draws with known true models. The Z ratio in Eq. (6) is a sample mean-to-standard-deviation summary of the resulting distribution of model selection statistics, not a quantity whose superiority is encoded in its definition. The paper is transparent that it extends the time-series findings of Cooper et al. (2024) rather than invoking that work as proof, and the citations to Cooper et al. and to Sivula et al. are used for motivation and for the Gaussian approximation of CV uncertainty, respectively, not as load-bearing substitutes for the simulation evidence. The Laplace approximation used for computation is admittedly approximate, and its validation is limited to the kernel experiments, but an approximation error is a correctness or robustness concern, not circularity: the paper does not define joint scoring so that it must outperform pointwise scoring. No fitted parameter is renamed as a prediction, no uniqueness theorem is imported from the authors' prior work, and no known empirical pattern is repackaged as a new derivation. Accordingly, the derivation chain is self-contained with respect to the stated empirical comparisons, and no circular step is exhibited.

Assumptions & free parameters 3 free parameters · 3 assumptions · 0 invented entities

The paper introduces no new postulated entities. It relies on standard Bayesian spatial modeling assumptions and the design assumptions of block CV for its simulations and case studies.

free parameters (3)

Test set size = s^2
Varied by experimental design, not used for prediction; it is a design variable.
rho (dgp) = 0.0 to 0.99
True autoregressive parameter used to simulate data and stratify results; not fitted to data.
halo order = 1 (order-1)
Chosen by hand as a design decision for all simulations.

assumptions (3)

domain assumption The block CV halo makes training and test sets sufficiently independent for the Gaussian approximation of the CV estimator to be valid.
Section 2 and Section 4: the Gaussian approximation for the Z ratio in Eq. (10) assumes many weakly correlated contributions with finite variance.
domain assumption The true data generating process for the simulations is a member of the Gaussian SAR/covariance kernel families considered.
Section 3: all simulations draw y from Gaussian SAR or Gaussian kernel models, so the results do not cover model misspecification beyond the candidate set.
domain assumption The Laplace approximation is accurate enough for the computed CV scores and model selection statistics.
Section 3 and Appendix D: Laplace approximation is used throughout; its accuracy is only illustrated for two single data draws, not across the simulation.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Joint leave-group-out cross-validation in Bayesian spatial models." pith.science (2026). https://pith.science/paper/N2LAHTVV

@misc{pith2026250415586,
  author       = {Pith},
  title        = {Pith review of: Joint leave-group-out cross-validation in Bayesian spatial models},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/N2LAHTVV}},
  note         = {Machine review of arXiv:2504.15586}
}

read the original abstract

Cross-validation (CV) is a widely-used method of predictive assessment based on repeated model fits to different subsets of the available data. CV is applicable in a wide range of statistical settings. However, in cases where data are not exchangeable, the design of CV schemes should account for suspected correlation structures within the data. CV scheme designs include the selection of left-out blocks and the choice of scoring function for evaluating predictive performance. This paper focuses on the impact of two scoring strategies for block-wise CV applied to spatial models with Gaussian covariance structures. We investigate, through several experiments, whether evaluating the predictive performance of blocks of left-out observations jointly, rather than aggregating individual (pointwise) predictions, improves model selection performance. Extending recent findings for data with serial correlation (such as time-series data), our experiments suggest that joint scoring reduces the variability of CV estimates, leading to more reliable model selection, particularly when spatial dependence is strong and model differences are subtle.

Figures

Figures reproduced from arXiv: 2504.15586 by the authors.

**Figure 2.** Summary of a sequence of model selection experiments to iden [PITH_FULL_IMAGE:figures/full_fig_p010_2.png] view at source ↗

**Figure 3.** Summary of a sequence of model selection experiments to identify [PITH_FULL_IMAGE:figures/full_fig_p012_3.png] view at source ↗

Figures from the paper (4 more)

**Figure 4.** Figure 4: Summary of a sequence of experiments selecting between two [PITH_FULL_IMAGE:figures/full_fig_p014_4.png]

**Figure 5.** Figure 5: A map of Australia, and six metro area regions. The colors [PITH_FULL_IMAGE:figures/full_fig_p016_5.png]

**Figure 6.** Figure 6: Summary of 105 pairwise comparisons among 15 candidate spatial [PITH_FULL_IMAGE:figures/full_fig_p017_6.png]

**Figure 7.** Figure 7: Summary of pairwise comparisons among 6 candidate spatial models [PITH_FULL_IMAGE:figures/full_fig_p020_7.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

3 extracted references · 2 canonical work pages

[20]

ApproximateBayesianinferenceforspatialeconometricsmod- els

issn: 1572-9052. doi: 10.1007/BF00116466. Bivand, Roger S., Virgilio Gomez-Rubio, and Håvard Rue (2014). “ApproximateBayesianinferenceforspatialeconometricsmod- els”. Spatial Statistics9, pp. 146–165. issn: 22116753. doi: 10.1016/j.spasta.2014.01.002. (Visited on 01/15/2024). Blondel, Mathieu et al. (2022). “Efficient and modular implicit differentiation”...

work page arXiv 2014
[1432]

Bayesian model as- sessment and comparison using cross-validation predictive densities

issn: 15731375. doi: 10.1007/s11222-016-9696-4. Vehtari, Aki and Jouko Lampinen (2002). “Bayesian model as- sessment and comparison using cross-validation predictive densities”. Neural Computation14.10, pp. 2439–2468.issn: 08997667. doi: 10.1162/08997660260293292. Vehtari, Aki, Tommi Mononen, et al. (2016). “Bayesian leave-one- out cross-validation approx...

arXiv 2002
[2740]

No jab, no play: Australia and compulsory vaccination

issn: 0885-6125, 1573-0565.doi: 10.1007/s10994-021- 05972-1. (Visited on 07/25/2023). 24 Kelter, Riko (2021). “Bayesian model selection in the M-open setting — approximate posterior inference and subsampling for efficient large-scale leave-one-out cross-validation via the difference estimator”. Journal of Mathematical Psychology 100, p. 102474.issn: 0022-...

work page arXiv 2021

Pith tools

Reviewed August 16, 2026 · model on record in the stance chip above.

[1] [20]

ApproximateBayesianinferenceforspatialeconometricsmod- els

issn: 1572-9052. doi: 10.1007/BF00116466. Bivand, Roger S., Virgilio Gomez-Rubio, and Håvard Rue (2014). “ApproximateBayesianinferenceforspatialeconometricsmod- els”. Spatial Statistics9, pp. 146–165. issn: 22116753. doi: 10.1016/j.spasta.2014.01.002. (Visited on 01/15/2024). Blondel, Mathieu et al. (2022). “Efficient and modular implicit differentiation”...

work page arXiv 2014

[2] [1432]

Bayesian model as- sessment and comparison using cross-validation predictive densities

issn: 15731375. doi: 10.1007/s11222-016-9696-4. Vehtari, Aki and Jouko Lampinen (2002). “Bayesian model as- sessment and comparison using cross-validation predictive densities”. Neural Computation14.10, pp. 2439–2468.issn: 08997667. doi: 10.1162/08997660260293292. Vehtari, Aki, Tommi Mononen, et al. (2016). “Bayesian leave-one- out cross-validation approx...

arXiv 2002

[3] [2740]

No jab, no play: Australia and compulsory vaccination

issn: 0885-6125, 1573-0565.doi: 10.1007/s10994-021- 05972-1. (Visited on 07/25/2023). 24 Kelter, Riko (2021). “Bayesian model selection in the M-open setting — approximate posterior inference and subsampling for efficient large-scale leave-one-out cross-validation via the difference estimator”. Journal of Mathematical Psychology 100, p. 102474.issn: 0022-...

work page arXiv 2021