Checking a forecast evaluation

forecasting

model evaluation

model diagnostics

ecology tutorial

Three checks for an ecological forecast evaluation in R: coverage rewards width, sharpness needs calibration, and small test sets cannot separate two forecasts.

Author

Tidy Ecology

Published

2026-07-24

The three earlier posts built the machinery: proper scores, the CRPS and Brier score, and skill against a baseline. This one is the checking step, the same closing move as every other cluster on this blog. A forecast evaluation can produce a clean headline number and still mislead, in three specific ways. Each has a diagnostic, and each diagnostic is a couple of lines of base R.

Check 1: coverage on its own rewards a wide forecast

A tempting way to grade a forecast is coverage: what fraction of outcomes fell inside the stated interval? It feels intuitive, and it is a trap, because it is improper. Rank forecasters by “more coverage is better” and the winner is whoever quotes the widest interval. A forecast that reports a spread four times too large covers almost everything, and by the coverage yardstick it beats an honest one.

set.seed(909)
mu_t <- 3; sig0 <- 1
y <- rnorm(5000, mu_t, sig0)
cover90 <- function(mu, s, y) { lo <- qnorm(0.05, mu, s); hi <- qnorm(0.95, mu, s)
  mean(y >= lo & y <= hi) }

honest <- c(coverage = cover90(mu_t, sig0, y),     crps = mean(crps_norm(mu_t, sig0, y)))
wide   <- c(coverage = cover90(mu_t, 4 * sig0, y), crps = mean(crps_norm(mu_t, 4 * sig0, y)))
rbind(honest, wide)

       coverage      crps
honest   0.9042 0.5614321
wide     1.0000 1.0321118

The wide forecast covers 1.000 of outcomes against the honest forecast’s 0.904, so by coverage alone it “wins”. But its CRPS is 1.84 times worse. Coverage is a calibration check, not a ranking rule: the target is coverage equal to the nominal level, not coverage as high as possible. A proper score penalises the wasted width that pure coverage rewards.

Two panels of bars for an honest and a wide forecast. The coverage panel shows both near one with the wide one slightly higher; the CRPS panel shows the wide forecast almost twice as high. — Figure 1: Ranking by coverage prefers the wide forecast, but its CRPS is far worse. Coverage checks calibration; it does not rank forecasts.

Check 2: sharpness is a virtue only after calibration

The opposite failure is prizing sharpness: narrow intervals look confident and precise. But a sharp forecast that is wrong is overconfident, and the sign of that is a probability integral transform (PIT) that is not flat. The PIT is the forecast CDF evaluated at each outcome; for a calibrated forecast it is uniform, and for an overconfident one it piles up in the tails, because outcomes keep landing outside the too-narrow intervals.

set.seed(919)
y2 <- rnorm(5000, mu_t, sig0)
width90 <- function(s) qnorm(0.95, 0, s) - qnorm(0.05, 0, s)
pit <- function(mu, s, y) pnorm(y, mu, s)

honest_s <- c(width = width90(sig0),       crps = mean(crps_norm(mu_t, sig0, y2)),
              tail = mean(pit(mu_t, sig0, y2) < 0.05 | pit(mu_t, sig0, y2) > 0.95))
sharp_s  <- c(width = width90(0.5 * sig0), crps = mean(crps_norm(mu_t, 0.5 * sig0, y2)),
              tail = mean(pit(mu_t, 0.5 * sig0, y2) < 0.05 | pit(mu_t, 0.5 * sig0, y2) > 0.95))
rbind(honest_s, sharp_s)

            width      crps   tail
honest_s 3.289707 0.5625069 0.0996
sharp_s  1.644854 0.6087213 0.4122

The overconfident forecast is sharper: its 90 percent interval is 1.64 wide against 3.29. If you chased sharpness you would prefer it. But 41 percent of outcomes fall in its tails instead of the nominal 10, and its CRPS is worse (0.609 against 0.563). The discipline, from Gneiting et al. (2007), is to maximise sharpness subject to calibration: a sharp forecast is better only once the PIT confirms it is honest. A single score cannot tell you which side of that line you are on, which is why the score and the calibration plot are reported together.

Two histograms of the probability integral transform. The calibrated one is flat around the dashed uniform line; the overconfident one is U-shaped with tall bars at both ends. — Figure 2: PIT histograms. Flat means calibrated; the U-shape of the overconfident forecast shows outcomes piling up in the tails.

Check 3: a small test set cannot separate two forecasts

Ecological forecasts are often evaluated on a handful of years or sites. A difference in mean score between two models is itself an estimate, with sampling uncertainty, and on a short test set that uncertainty can swamp the difference. Here model A is the honest forecast and model B is slightly too wide; A is genuinely better. Bootstrapping the per-observation CRPS difference across the test points shows when the evidence is actually there (Diebold and Mariano 1995).

crps_diff_ci <- function(ntest, seed, B = 4000) {
  set.seed(seed)
  yy <- rnorm(ntest, mu_t, sig0)
  d <- crps_norm(mu_t, 1.4 * sig0, yy) - crps_norm(mu_t, sig0, yy)   # B minus A, positive = A better
  bs <- replicate(B, mean(sample(d, replace = TRUE)))
  c(mean = mean(d), lo = quantile(bs, 0.025), hi = quantile(bs, 0.975))
}
t(sapply(c(10, 30, 100, 300), function(nt) crps_diff_ci(nt, seed = 1000 + nt)))

             mean     lo.2.5%   hi.97.5%
[1,] -0.005584860 -0.06403161 0.04713802
[2,]  0.007761495 -0.02474203 0.03743723
[3,]  0.025762271  0.01115841 0.04012961
[4,]  0.023175917  0.01465637 0.03140212

At 10 test points the interval spans zero and the point estimate even comes out the wrong sign; at 30 it still spans zero; only at 100 and beyond does the interval clear zero and confirm that A beats B. A ranking of two forecasts on a dozen years is not evidence of a real difference. The score can favour one model on a small sample purely by chance, so the interval, not the point estimate, is what settles it.

Point estimates with error bars for the mean CRPS difference at test sizes ten, thirty, one hundred and three hundred. The first two bars cross the zero line; the last two sit clearly above it. — Figure 3: Bootstrap intervals for the CRPS difference between two forecasts. Below about a hundred test points the difference is indistinguishable from zero.

Where this leaves the series

Put the four posts together and the workflow is complete. Forecast a distribution, not a point. Grade it with a proper score: the Brier score for presence, the CRPS for abundance. Turn the raw score into skill against a baseline that is genuinely hard, and report the baseline. Then check the evaluation itself: coverage is a calibration test and not a ranking, sharpness only counts once the PIT is flat, and a score difference on a short record needs an interval before you believe it. None of this makes a bad model good, but it stops a good score from hiding a bad forecast.

References

Gneiting T, Balabdaoui F, Raftery AE 2007. Journal of the Royal Statistical Society B 69(2):243-268 (10.1111/j.1467-9868.2007.00587.x)

Broecker J, Smith LA 2007. Weather and Forecasting 22(2):382-388 (10.1175/WAF966.1)

Diebold FX, Mariano RS 1995. Journal of Business and Economic Statistics 13(3):253-263 (10.1080/07350015.1995.10524599)

Dietze MC, Fox A, Beck-Johnson LM, et al. 2018. Proceedings of the National Academy of Sciences 115(7):1424-1432 (10.1073/pnas.1710231115)

--- title: "Checking a forecast evaluation" description: "Three checks for an ecological forecast evaluation in R: coverage rewards width, sharpness needs calibration, and small test sets cannot separate two forecasts." date: "2026-07-24 15:00" categories: [R, forecasting, model evaluation, model diagnostics, ecology tutorial] image: thumbnail.png image-alt: "Two probability integral transform histograms, one flat for a calibrated forecast and one U-shaped for an overconfident forecast." --- The three earlier posts built the machinery: [proper scores](../scoring-ecological-forecasts/), [the CRPS and Brier score](../the-crps-and-brier-score/), and [skill against a baseline](../forecast-skill-and-the-baseline/). This one is the checking step, the same closing move as every other cluster on this blog. A forecast evaluation can produce a clean headline number and still mislead, in three specific ways. Each has a diagnostic, and each diagnostic is a couple of lines of base R. ```{r setup} #| echo: false #| message: false #| warning: false suppressMessages(library(ggplot2)) col_main <- "#275139"; col_acc <- "#B5651D"; col_mut <- "#8C8C8C" base_th <- theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), plot.title = element_text(face = "bold", size = 13), axis.title = element_text(size = 11)) crps_norm <- function(mu, sigma, y) { z <- (y - mu) / sigma sigma * (z * (2 * pnorm(z) - 1) + 2 * dnorm(z) - 1 / sqrt(pi)) } ``` ## Check 1: coverage on its own rewards a wide forecast A tempting way to grade a forecast is coverage: what fraction of outcomes fell inside the stated interval? It feels intuitive, and it is a trap, because it is improper. Rank forecasters by "more coverage is better" and the winner is whoever quotes the widest interval. A forecast that reports a spread four times too large covers almost everything, and by the coverage yardstick it beats an honest one. ```{r coverage} set.seed(909) mu_t <- 3; sig0 <- 1 y <- rnorm(5000, mu_t, sig0) cover90 <- function(mu, s, y) { lo <- qnorm(0.05, mu, s); hi <- qnorm(0.95, mu, s) mean(y >= lo & y <= hi) } honest <- c(coverage = cover90(mu_t, sig0, y), crps = mean(crps_norm(mu_t, sig0, y))) wide <- c(coverage = cover90(mu_t, 4 * sig0, y), crps = mean(crps_norm(mu_t, 4 * sig0, y))) rbind(honest, wide) ``` The wide forecast covers `r sprintf("%.3f", cover90(mu_t, 4 * sig0, y))` of outcomes against the honest forecast's `r sprintf("%.3f", cover90(mu_t, sig0, y))`, so by coverage alone it "wins". But its CRPS is `r sprintf("%.2f", mean(crps_norm(mu_t, 4 * sig0, y)) / mean(crps_norm(mu_t, sig0, y)))` times worse. Coverage is a calibration check, not a ranking rule: the target is coverage equal to the nominal level, not coverage as high as possible. A proper score penalises the wasted width that pure coverage rewards. ```{r fig-coverage} #| echo: false #| fig-cap: "Ranking by coverage prefers the wide forecast, but its CRPS is far worse. Coverage checks calibration; it does not rank forecasts." #| fig-alt: "Two panels of bars for an honest and a wide forecast. The coverage panel shows both near one with the wide one slightly higher; the CRPS panel shows the wide forecast almost twice as high." #| fig-width: 6.6 #| fig-height: 3.6 d4l <- data.frame( forecast = rep(c("honest", "wide (4x sd)"), 2), metric = rep(c("coverage", "crps"), each = 2), val = c(honest["coverage"], wide["coverage"], honest["crps"], wide["crps"])) ggplot(d4l, aes(forecast, val, fill = forecast)) + geom_col(width = 0.6) + facet_wrap(~metric, scales = "free_y") + scale_fill_manual(values = c(col_main, col_acc), guide = "none") + labs(title = "Ranking by coverage picks the worse (wider) forecast", x = NULL, y = NULL) + base_th ``` ## Check 2: sharpness is a virtue only after calibration The opposite failure is prizing sharpness: narrow intervals look confident and precise. But a sharp forecast that is wrong is overconfident, and the sign of that is a probability integral transform (PIT) that is not flat. The PIT is the forecast CDF evaluated at each outcome; for a calibrated forecast it is uniform, and for an overconfident one it piles up in the tails, because outcomes keep landing outside the too-narrow intervals. ```{r sharpness} set.seed(919) y2 <- rnorm(5000, mu_t, sig0) width90 <- function(s) qnorm(0.95, 0, s) - qnorm(0.05, 0, s) pit <- function(mu, s, y) pnorm(y, mu, s) honest_s <- c(width = width90(sig0), crps = mean(crps_norm(mu_t, sig0, y2)), tail = mean(pit(mu_t, sig0, y2) < 0.05 | pit(mu_t, sig0, y2) > 0.95)) sharp_s <- c(width = width90(0.5 * sig0), crps = mean(crps_norm(mu_t, 0.5 * sig0, y2)), tail = mean(pit(mu_t, 0.5 * sig0, y2) < 0.05 | pit(mu_t, 0.5 * sig0, y2) > 0.95)) rbind(honest_s, sharp_s) ``` The overconfident forecast is sharper: its 90 percent interval is `r sprintf("%.2f", width90(0.5 * sig0))` wide against `r sprintf("%.2f", width90(sig0))`. If you chased sharpness you would prefer it. But `r sprintf("%.0f", 100 * sharp_s["tail"])` percent of outcomes fall in its tails instead of the nominal 10, and its CRPS is worse (`r sprintf("%.3f", sharp_s["crps"])` against `r sprintf("%.3f", honest_s["crps"])`). The discipline, from Gneiting et al. (2007), is to maximise sharpness *subject to* calibration: a sharp forecast is better only once the PIT confirms it is honest. A single score cannot tell you which side of that line you are on, which is why the score and the calibration plot are reported together. ```{r fig-pit} #| echo: false #| fig-cap: "PIT histograms. Flat means calibrated; the U-shape of the overconfident forecast shows outcomes piling up in the tails." #| fig-alt: "Two histograms of the probability integral transform. The calibrated one is flat around the dashed uniform line; the overconfident one is U-shaped with tall bars at both ends." #| fig-width: 7.2 #| fig-height: 3.6 n2 <- length(y2) d4b <- data.frame(pit = c(pit(mu_t, sig0, y2), pit(mu_t, 0.5 * sig0, y2)), who = rep(c("calibrated", "overconfident (sharp)"), each = n2)) ggplot(d4b, aes(pit)) + geom_histogram(aes(fill = who), bins = 20, boundary = 0, colour = "white") + geom_hline(yintercept = n2 / 20, linetype = "dashed", colour = col_mut) + facet_wrap(~who) + scale_fill_manual(values = c(col_main, col_acc), guide = "none") + labs(title = "PIT histogram: flat is calibrated, U-shape is overconfident", x = "probability integral transform", y = "count") + base_th ``` ## Check 3: a small test set cannot separate two forecasts Ecological forecasts are often evaluated on a handful of years or sites. A difference in mean score between two models is itself an estimate, with sampling uncertainty, and on a short test set that uncertainty can swamp the difference. Here model A is the honest forecast and model B is slightly too wide; A is genuinely better. Bootstrapping the per-observation CRPS difference across the test points shows when the evidence is actually there (Diebold and Mariano 1995). ```{r diffci} crps_diff_ci <- function(ntest, seed, B = 4000) { set.seed(seed) yy <- rnorm(ntest, mu_t, sig0) d <- crps_norm(mu_t, 1.4 * sig0, yy) - crps_norm(mu_t, sig0, yy) # B minus A, positive = A better bs <- replicate(B, mean(sample(d, replace = TRUE))) c(mean = mean(d), lo = quantile(bs, 0.025), hi = quantile(bs, 0.975)) } t(sapply(c(10, 30, 100, 300), function(nt) crps_diff_ci(nt, seed = 1000 + nt))) ``` At 10 test points the interval spans zero and the point estimate even comes out the wrong sign; at 30 it still spans zero; only at 100 and beyond does the interval clear zero and confirm that A beats B. A ranking of two forecasts on a dozen years is not evidence of a real difference. The score can favour one model on a small sample purely by chance, so the interval, not the point estimate, is what settles it. ```{r fig-diffci} #| echo: false #| fig-cap: "Bootstrap intervals for the CRPS difference between two forecasts. Below about a hundred test points the difference is indistinguishable from zero." #| fig-alt: "Point estimates with error bars for the mean CRPS difference at test sizes ten, thirty, one hundred and three hundred. The first two bars cross the zero line; the last two sit clearly above it." #| fig-width: 6.4 #| fig-height: 3.8 ns <- c(10, 30, 100, 300) M <- t(sapply(ns, function(nt) crps_diff_ci(nt, 1000 + nt))) d4c <- data.frame(n = factor(ns), mean = M[, 1], lo = M[, 2], hi = M[, 3]) ggplot(d4c, aes(n, mean)) + geom_hline(yintercept = 0, colour = col_mut) + geom_errorbar(aes(ymin = lo, ymax = hi), width = 0.2, colour = col_main) + geom_point(size = 2.4, colour = col_main) + labs(title = "Small test sets cannot tell two forecasts apart", x = "number of forecast-observation pairs", y = "mean CRPS difference (B - A)") + base_th ``` ## Where this leaves the series Put the four posts together and the workflow is complete. Forecast a distribution, not a point. Grade it with a proper score: the Brier score for presence, the CRPS for abundance. Turn the raw score into skill against a baseline that is genuinely hard, and report the baseline. Then check the evaluation itself: coverage is a calibration test and not a ranking, sharpness only counts once the PIT is flat, and a score difference on a short record needs an interval before you believe it. None of this makes a bad model good, but it stops a good score from hiding a bad forecast. ## Related tutorials - [Forecast skill and the baseline](../forecast-skill-and-the-baseline/) - [The CRPS and the Brier score](../the-crps-and-brier-score/) - [Scoring ecological forecasts](../scoring-ecological-forecasts/) - [Checking predictive calibration](../checking-predictive-calibration/) ## References Gneiting T, Balabdaoui F, Raftery AE 2007. Journal of the Royal Statistical Society B 69(2):243-268 (10.1111/j.1467-9868.2007.00587.x) Broecker J, Smith LA 2007. Weather and Forecasting 22(2):382-388 (10.1175/WAF966.1) Diebold FX, Mariano RS 1995. Journal of Business and Economic Statistics 13(3):253-263 (10.1080/07350015.1995.10524599) Dietze MC, Fox A, Beck-Johnson LM, et al. 2018. Proceedings of the National Academy of Sciences 115(7):1424-1432 (10.1073/pnas.1710231115)