Score cross-validation on the posterior predictive distribution

[goldingn]

Drafted by Claude Code, opened from Nick's account.

Implements the plan in #10, and closes #11. Draft: the code path is complete and tested end to end on a real fold, but the full set of folds has not been re-run. Needs the greta version pin described below.

What changes

Cross-validation previously scored the model on its ability to predict the value of an individual held-out bioassay, with the posterior collapsed to a mean predicted fraction and a mean overdispersion before scoring. An individual bioassay is a noisy, overdispersed measurement of the quantity the model actually targets, so that comparison has an unreachable floor and discards parameter uncertainty entirely. This scores the full out-of-sample posterior predictive distribution instead, and reports calibration alongside accuracy.

New files

file	purpose
R/validation_functions.R	beta-binomial primitives, randomised quantile residuals, coverage curves, CRPS, reliability bins, the noise floor, aggregated predictive distributions
R/check_validation_functions.R	checks of the above against known moments and simulated data
R/estimate_bioassay_rho.R	overdispersion estimated from replicate bioassays
R/validation_folds.R	the three cross-validation splits, extracted so all models share them
R/null_models.R	the two null models, extracted, plus predictive distributions for them
R/validation_covariates.R	fold-independent covariate setup, extracted
R/fit_validation_fold.R	the dynamical model definition, once, returning posterior draws
R/run_validation_folds.R	runs every fold of all three experiments and saves the draws
R/validation_metrics.R	scores the saved draws, writes tidy tables
R/fig_predictive_validation.R	figures and the results table

predictive_validation.R, dynamic_predictive_validation.R and their outputs are left in place, to be retired once the new path has been run.

Points worth attention in review

The posterior is no longer collapsed before scoring. run_validation_folds.R saves the draws of the predicted fraction and the overdispersion rather than their means. No extra MCMC, but all folds must be re-run, since only means were kept last time.

The predicted fraction and the overdispersion are now drawn together. The earlier code took them from two separate calculate() calls, leaving them unpaired across posterior samples. That did not matter when only their means were used; a predictive distribution needs them coupled.

Overdispersion re-estimated from every replicated pixel-year. It was previously fitted to the six most heavily sampled pixel-year-insecticide combinations. Refitting over all 3,713 replicated groups (9,420 assays) gives 0.155 (0.149–0.161), against 0.195 (0.147–0.254) from the six-group version on current data, and the value of 0.14 hard-coded in betabinom_dev(). Class-specific estimates span nearly a factor of two, so a single pooled value is not adequate:

	rho	95% CI	assays
all	0.155	0.149–0.161	9,420
Organophosphates	0.216	0.182–0.254	836
Pyrethroids	0.164	0.156–0.171	5,947
Carbamates	0.131	0.113–0.151	1,180
Organochlorines	0.118	0.106–0.131	1,457

The group fractions are integrated out under a beta hyperprior rather than maximised over: with most groups holding two assays, maximising over one fraction per group is subject to the incidental parameters problem and biases rho downwards.

The nulls now carry predictive distributions. Scoring them at a fixed overdispersion while the dynamical model used its own fitted value let a model win on dispersion rather than on prediction. Each null now has a beta posterior on its fraction and an overdispersion fitted to its own training residuals — for the nearest neighbour model on the internal holdout that already selects the number of neighbours, so nothing is tuned on the test fold.

Three copies of the model definition became one. fit_validation_fold.R holds it. Because its greta arrays are function-local, the manual purge_greta_model() between folds is no longer needed. This is the change most worth checking against the original blocks.

The noise floor is in second-moment space, not information space. As set out in #10, the oracle log score varies by 2.4–3.5 nats across the plausible range of p, so it cannot anchor a scale. Mean squared error can, because only p(1-p) is needed and that is recoverable from the observed proportion. Corrected figures, replacing those in the #10 comment (which used rho = 0.14 and a null MSE computed against the overall rather than per-insecticide mean):

noise floor, class-specific rho	MSE 0.0201	RMSE 0.142
insecticide-type-mean null	MSE 0.0831	RMSE 0.288

so assay noise is 24% of the null's scatter and 76% is in principle explainable.

What is verified

check_validation_functions.R passes 26 checks: the primitives against known moments and the binomial limit, agreement with extraDistr to machine precision and with scoringRules to 1e-16, and the metrics against simulated data where the truth is known — flat for a model whose predictive distribution is the data-generating one, deviating in the expected direction when biased, overconfident or underconfident, and the noise floor recovering the irreducible error.

The extracted fold definitions in validation_folds.R produce splits identical to the originals on master, checked by running both and comparing the fold membership: 6 country folds (1133, 1563, 1638, 913, 1486, 1961 held-out assays), interpolation 24318/1045, forecasting 25900/1461.

The scoring and figure scripts were run end to end against synthetic draws in the saved format, so the path from draws to tables to figures works. Those artefacts were deleted rather than committed, since their content is not real.

greta version

greta 0.6.0 cannot sample a model whose likelihood depends on iterate_dynamic_function() output. Model construction succeeds; the error is raised at mcmc(), which is why it is invisible to anything that only builds the model. Bisected by holding greta.dynamics fixed and moving only greta:

greta	greta.dynamics	result
0.5.0.9000 (njtierney/greta@4cc989f)	db7df31	samples
0.5.0 (CRAN)	db7df31	samples
0.5.0.9000	v0.2.2 (e30d9ad)	samples
0.6.0	db7df31	fails
0.6.0	0.2.3 (CRAN)	fails

Which of two TensorFlow errors appears depends on the greta.dynamics version (enters the loop with shape … or Shape must be rank 3 but is rank 2 … while/Tile), but 0.6.0 fails either way. The same single change flips a fifteen-line example from sampling to failing, so this is not specific to this model, and it is not caused by this branch: master's model definition behaves identically.

packages.R now pins the working combination — tensorflow 2.16.0, njtierney/greta@4cc989f, greta-dev/greta.dynamics@db7df31 — replacing the reference to the greta_2 branch, which no longer exists upstream. On those versions both master's model definition and fit_fold() sample a real validation fold end to end.

Worth reporting upstream; a reprex is ready.

To run

source("R/estimate_bioassay_rho.R")     # minutes; writes outputs/bioassay_rho.csv
source("R/run_validation_folds.R")      # 8 HMC fits: 6 countries, interpolation, forecasting
source("R/validation_metrics.R")        # minutes
source("R/fig_predictive_validation.R")

A 20-iteration test fit of one country fold took about two minutes wall clock on 16 cores, much of it graph construction. The real settings are 8 chains and 3000 iterations, so one fold is worth timing before launching all eight. Folds whose draws already exist are skipped, so the run can be interrupted and resumed, and the first fold is worth checking against the current posterior means before committing to the rest.

Open question

The reliability envelope and the pixel-year aggregation both rest on replicated pixel-years, and those are not fully representative: their mean mortality is 0.685 against 0.744 for singly sampled pixel-years, and the replicated share ranges from 0.21 in Ethiopia to 0.56 in Cameroon. Whether the overdispersion itself differs between sentinel and other sites cannot be tested directly, since estimating it requires replicates. figures/bioassay_rho_representativeness.png records the comparison; how much weight the paper should put on the floor given this is a judgement call.