A model with R-squared near 0 can still give valid 90% prediction intervals - here's why (and the catch)

I recently calibrated a recovery-rate model that had only two weak features. Its point accuracy was almost nothing — R² basically zero. I expected its uncertainty estimates to be junk too. They weren't: the 90% conformal prediction intervals covered ~89% of held-out outcomes. Valid, just wide . That surprised me enough to nail it down, because it contradicts a belief a lot of us carry around: "my model isn't accurate, so I can't trust its uncertainty." For split conformal prediction, that's backwards. Here's the precise statement, a runnable demo, and the one caveat that actually bites. Coverage is a property of the procedure, not the model Split conformal prediction gives a distribution-free, finite-sample marginal coverage guarantee : P( Y ∈ Ĉ(X) ) ≥ 1 − α and it holds for any point model, as long as the calibration and test data are exchangeable. The model is a black box. You fit it however you like, then on a held-out calibration set you take the (1−α) quantile of the absolute residuals, and that quantile becomes the half-width of your intervals. Nowhere does that construction require the model to be good. A bad model just has large residuals, so the calibration quantile is large, so the intervals are wide — wide enough to still cover at the stated rate. Accuracy doesn't buy you validity ; it buys you efficiency (narrower intervals at the same coverage). The demo (numbers are reproducible, seed fixed) Same dataset and target, three models from strong to useless, target coverage 90%: model R² marginal coverage mean interval width gradient boosting 0.741 0.895 5.39 weak linear (1 noisy feature) 0.061 0.905 10.39 predict-the-mean −0.000 0.907 10.83 All three land at ~90% coverage. The only thing that changes is width: the good model's intervals are half as wide . That's the whole story in one table — validity is constant, efficiency tracks accuracy. import numpy as np from sklearn.linear_model import LinearRegression from sklearn.ensemble import GradientBoostingReg