In some of mypapers I’ve used 80% predictive intervals instead of the standard 95% predictive intervals. I didn’t say why in the papers and I was thinking through it again the other day so I thought I’d write it down. The focus here is on prediction intervals (the range of predicted values supported by the model) rather than confidence/credible intervals (the range of parameter values supported by the model). Some of the arguments might transfer, but probably not all of them.
I think most people know that 95% is an arbitrary number plucked out of thin air by an unsavoury racist one hundred years ago. So there’s no rule that says we must use 95%, but how then do we decide which value to use? I’m fairly convinced that not all prediction intervals are created equal. A 40% prediction interval (assuming it’s not presented alongside other intervals) is an odd metric that represents a high density interval that the true value probably isn’t in (<50% chance).
Given that, here are some things to consider when chosing an interval. I’ll use the binary decision of 80% versus 95% just to illustrate things. And I will ignore the fact that of course we often want to use the full distribution rather than summarise a distribution as a single interval; there are so many cases where communicating a full distribution is not feasible. If you can use 2 or 3 intervals, or density plots of the full distribution, do that; if you must summarise a distribution as a single interval, perhaps consider these points.
How do policy makers interpret intervals?
This point is probably the most subjective, but possibly the most important. Despite our best efforts, I think most people, even trained scientists, think of a 95% prediction interval as “the true value is almost certainly in this interval”. Obviously, this isn’t true. 1 in 20 95% prediction intervals will not cover the true value. Furthermore, we are talking about prediction intervals, and we will almost certainly be making hundreds of predictions. Therefore, we must expect many of our prediction intervals not to cover their respective true values.
I think perhaps this problem doesn’t exist for an 80% interval. I don’t think people look at an 80% prediction interval and think “the true value is almost certainly in this interval”. It’s more like a “best guess” interval. Perhaps this is more useful. But to reiterate, this is a very subjective point and I have no hard evidence for this point. The counter argument is that perhaps people can’t switch between intervals easily, and that therefore we should somehow settle on a single interval; given it’s history, 95% would be the clear winner in this case.
Wide and confident or thin and unsure?
This point is similar to that above but less subjective. It is a question of what is really useful in a prediction interval. Do we want a very wide interval that we’re 95% sure the true value is within, or do we want a smaller interval that we are 80% sure contains the true value.
Prediction intervals can quickly become massive and controlling this size can sometimes usefully guide our decisions. For example, imagine you are diagnosed with a terminal disease and told that your 95% predictive survival interval is between 1 and 15 years (I’m thinking for myself as a 30 something. I guess the details will change depending on your age). What can you actually do with this information? On the lower end you have enough time to sort your affairs, visit some family and take a great last trip. On the upper end you have enough time to do anything really; start a new career, see your kids grow up, die from a variety of other causes.
Imagine instead you were given an 80% interval of 3 to 6 years. This “probably true” interval is quite useful. You have some time, you don’t have prioritise only 3 or 4 things to do before you die. But starting a new career may well be a waste of time. Of course, all the survival times that were in the 95% interval are still possible, but this “probably true” interval focusses on a reasonable range of very likely values.
Relatedly, prediction intervals typically get wider faster as you increase the probability interval (look at the shape of a normal distribution for intuition). In a normal distribution, a 95% interval is more than 53% bigger than an 80% interval. Is the extra 15% confidence that the true value is in the interval, worth the 50% increase in width? I’d say often it isn’t.
Confident and wrong or unsure and right?
There are many reasons why we would expect an 80% prediction interval to have better calibration than a 95% prediction interval. Model mispecification and approximations or MCMC methods that break down in the tails are two examples. I would prefer a well calibrated 80% interval to a poorly calibrated 95% interval.
With respect to this point, we can consider a few related situations. Did we just choose one interval from the outset and only check the calibration of that one interval? If so, I think it is reasonable to consider which interval is likely to be better calibrated. As long as we don’t claim that we have evidence that the model is well calibrated across all intervals, we have tested one aspect of the model and found it to be adequate. Perhaps instead, we are testing the calibration of the model with respect to both 80% and 95% prediction intervals. How is it reasonable to behave if we find the model is well calibrated for 80% intervals and badly calibrated for 95% intervals? Again, I think it is totally fine to recommend that users of the mode use the 80% interval. This is similar to saying “linear regression works well as long as you don’t extrapolate far outside the range of the covariates”. We are guiding the user as to when the model does and does not work and again this is totally fine.
Confident but badly estimated coverage or unsure but well estimated coverage?
My last point is that estimating the calibration of a model is easier when the interval is smaller. In the same way as having few cases vs controls makes our effective sample size small, having large prediction intervals gives us fewer data points where the prediction intervals do not cover the true value. For example, if our dataset contains 200 datapoints, a well calibrated model will have around 10 datapoints where the 95% prediction interval does not cover the true value. In contrast, the 80% prediction interval would give us 40 failures, a much more reasonable sample size. So if we are working with modest sample sizes, would we prefer a 95% prediction interval, where our estimates of coverage are very noisey, or an 80% prediction interval with much tighter estimates of coverage. In the above example with n = 200, our 95% confidence interval of our coverage, if we observed exactly 5% of intervals to not cover their true values, would be 2.4% - 9%. I would consider 2.4% or 9% to imply fairly poor calibration, so in this case we are really unsure whether our model is well calibrated or not. In contrast, if we observed exactly 20% of intervals to not cover their true values, our confidence intervals for the coverage of the 80% interval would be 14% - 26%. So we’re pretty sure the coverage for the 80% prediction interval is ok.
Conclusions
So in conclusion, perhaps we should think about what intervals we use a bit more. Different considerations will apply in different situations, and almost certainly there are different considerations for prediction intervals and confidence/credible intervals. As in the intro, we should avoid summarise full distributions when we can, but often we can’t.
A practical guide to squishing yearly ‘disag data’ objects
A practical guide to pooling annual data for Bayesian disaggregation regression in R
This post is a guest post by my predoctoral fellow, Samana.
Bayesian spatial disaggregation gives you a high resolution risk surface from low-resolution case counts ….. but only one year at a time. This post shows how to pool several years into a single model fit.
Hi all, my name is Samana and I’m a predoctoral fellow and I’ve been working on the problem of a multi-year disaggregation model!
In this blog I will be presenting a recipe of how to do a multi-year disaggregation model fit.
Currently, using the disaggregation package we can only really fit one model per time slice of data we might have, but this blog aims to change that and let you fit a disag model with multiple years of data.
Part A demonstrates the mechanics without a time trend; Part B adds a time spline to the six-year model. This blog is part A.
This post shows how to combine annual disag_data objects and fit one pooled spatial disaggregation model. It demonstrates shared effects across years and it is not a full spatio-temporal model with separate spatial field for each year.
The resolution mismatch problem:
Simply put, there is a resolution mismatch between satellite covariates (which are very high resolution) and public health data, like case counts of a particular disease in a specific county (this is low resolution).
Public health data is often reported as annual case counts for administrative areas such as counties and districts whereas environmental information that may explain risk, such as temperature, land cover, vegetation is often available on a fine spatial grid.
The problem is that we can’t just combine these together and then fit a model for risk predictions to get a high resolution risk surface. This is because of the ecological fallacy.
These figures above show the resolution mismatch; on the left we have the low resolution public health case counts and on the right we have the high resolution environmental satellite covariates (Nandi et al. 2023, https://doi.org/10.18637/jss.v106.i11).
What is disaggregation and why you’d want to do it:
So these resolutions are linked by modelling a latent disease risk at the pixel level and THEN aggregating the pixel level predictions to the polygon level so they’re consistent with the observed disease counts. all of this is done by Bayesian spatial disaggregation (implemented in the R package disaggregation )
The package fits a Bayesian disaggregation regression model (similar to a GLM) with an INLA style spatial random field via TMB, using the function prepare_data() to line up polygons, covariates and population, and disag_model() to fit the model.
The current single year workflow is like this:
Read in the inputs: a sf object containing the area boundaries and case counts, a SpatRaster of fine-scale covariates, and a population raster when modelling counts with population as exposure.
Prepare these inputs into a coherent data object with prepare_data().
Fit the model with disag_model(), choosing the likelihood and link and whether to include a spatial field and an independent (iid) effect.
Inspect the fit with summary() and plot(), then use predict() to produce fine-scale prediction and uncertainty maps.
For multiple years, the missing step is between 2 and 3: prepare each year separately, then combine the resulting disag_data objects carefully before fitting one pooled model.
Why multi-year disaggregation?
When you have many years of polygon level counts, a natural question is: fit each year separately, or pool several years into one model?
Pooling helps when you want the spatial field and covariate effects to borrow strength across years, while still letting a time trend move through a flexible term (we’ve used a natural spline on centered year in part B of this blog).
I’ve been calling this “squishing” years together, stacking each year’s prepared data into one long object before handing it to disag_model().
Bayesian spatial disaggregation is straight forward for one year of data but multi-year models are harder because the disag_data object contains many components that need to be in the right shape, format and size.
This object contains polygon-level outcomes, pixel-level covariates, aggregation weights, spatial coordinates, start end indices and the spatial mesh. These objects must remain correctly aligned when the years are combined.
This tutorial develops a reproducible workflow for combining yearly disag_data objects informally, “squishing” or “combining” them before fitting a multi-year Bayesian disaggregation model.
Model components
The disaggregation regression models demonstrated here include fixed effects of covariates and a number of random effects, for various different modelling purposes.
In this example, the pooled model estimates shared climate effects. We have different covariate observations for each year, but want to fit a model using all years of data in one go.
The spatial random effect operates at high resolution. If you model each year separately, then each time period has its own spatial field, defined over the mesh and projected to the pixels. It captures smooth residual spatial variation that remains after accounting for measured environmental covariates. However, we often want to pool estimates for this random effect across years.
The iid random effect operates at polygon level, the lower resolution. It gives each polygon-time observation an additional independent residual term, helping account for overdispersion and unmeasured area-level factors.
When yearly data are squished, the model is fitted jointly across all periods. The workflow must preserve the relationships between each polygon, its pixels, its aggregation weights and its appropriate spatial field.
The iid term allows an additional residual effect for each polygon-year observation. Stacking the data does not automatically create a separate spatial field for each year or model how fields evolve over time.
THE DATA:
The example region is Madagascar, deliberately chosen because it’s the same country used in the disaggregation package’s own founding methods paper (Nandi et al. 2020, malaria case counts across Madagascar, [Nandi et al](https://www.jstatsoft.org/article/view/v106i11)) so the mesh settings and covariate naming are the same as that original worked example.
The geometry, climate covariates AND population is real public data; only the case counts are simulated (real disease data can’t be shared here for data-permission reasons)
The dataset pairs real, public district geometry, real climate covariates and real population with simulated case counts:
Real public administrative level-2 (district) boundaries for Madagascar, from the GADM database (https://gadm.org), downloaded via the R geodata package.names.
Year varying covariates (temp_wc, precip_wc): annual mean temperature and annual total precipitation for each simulated year, built from REAL monthly TerraClimate records (Abatzoglou et al. 2018) not long-run climate normals, so these differ from one year to the next. Centred and scaled (mean 0, sd 1, using one mean/sd computed across all years combined so real year-to-year differences are preserved). Source: https://www.climatologylab.org/terraclimate.html.
Simulated case counts drawn from known truth parameters (shipped alongside the data as true_parameters.csv), using the real covariates and population above plus a simulated spatial random field and district-level noise. So this post can check whether the fitted model recovers the true coefficients, a much better teaching device than fitting to arbitrary noise.
None of the numbers in the dataset correspond to any real disease case count.
:::
### Step 4: prepare data for each year
Mesh args here match those used in the `disaggregation` package's own
Madagascar malaria worked example (Nandi et al. 2020) a nice side
benefit of picking the same country is that these don't need retuning
from scratch.
Always retune these to your own study area if you swap in a different
region.
The details for what these parameters mean can be found in the inla
tutorials for example.
::::::: cell
``` {.r .cell-code}
mesh_args_mdg <- list(max.edge = c(0.7, 8), cutoff = 0.05, offset = c(1, 2))
dis_2001 <- prepare_data(
polygon_shapefile = shapes_2001,
covariate_rasters = cov_stack_2001,
aggregation_raster = population_raster,
mesh_args = mesh_args_mdg,
id_var = "district_id",
response_var = "cases",
na_action = TRUE
)
::: {.cell-output .cell-output-stderr} Warning: [rast] CRS do not match :::
::: {.cell-output .cell-output-stderr} Warning: [extract] transforming vector data to the CRS of the raster :::
::: {.cell-output .cell-output-stderr}
Warning: [rast] CRS do not match
Warning: [extract] transforming vector data to the CRS of the raster
:::
``` {.r .cell-code}
dis_2003 <- prepare_data(
polygon_shapefile = shapes_2003,
covariate_rasters = cov_stack_2003,
aggregation_raster = population_raster,
mesh_args = mesh_args_mdg,
id_var = "district_id",
response_var = "cases",
na_action = TRUE
)
::: {.cell-output .cell-output-stderr} Warning: [rast] CRS do not match Warning: [extract] transforming vector data to the CRS of the raster ::: :::::::
Step 5: The squish: stack the 3 years of data into one disag data object
To stack these objects together correctly, we handle different parts of the data differently. polygon_data and covariate_data are just stacked with a year label attached. We offset start_end_index: each prepare_data() call starts its pixel-row numbering again, so the second and third years’ indices must be shifted to point to their new positions in the combined covariate table.
::: cell ``` {.r .cell-code}
polygon_data and covariate_data just stack with a year label attached
poly_0103 <- bind_rows( transform(dis_2001$polygon_data, year = 2001), transform(dis_2002$polygon_data, year = 2002), transform(dis_2003$polygon_data, year = 2003) )
cov_0103 <- bind_rows( transform(dis_2001$covariate_data, year = 2001), transform(dis_2002$covariate_data, year = 2002), transform(dis_2003$covariate_data, year = 2003) )
:::
### step 6: fit the pooled model
:::::: cell
``` {.r .cell-code}
fit_0103 <- disaggregation::disag_model(
data = dis_0103,
iterations = 1000,
field = TRUE,
iid = TRUE,
family = "poisson",
link = "log"
)
::: {.cell-output .cell-output-stderr} Fitting model. This may be slow. :::
Compare that against three separate single-year fits (disag_model() on 2001, 2002, 2003 individually, no squishing) to see whether pooling changes the covariate estimates, tightens their standard errors, or shifts the spatial field. This comparison is the actual point of the exercise the code above just gets you to a model you can compare against.
:::::: cell ``` {.r .cell-code} fit_2001 <- disaggregation::disag_model(dis_2001, iterations = 1000, field = TRUE, iid = TRUE, family = “poisson”, link = “log”)
::: {.cell-output .cell-output-stderr}
Fitting model. This may be slow.
:::
``` {.r .cell-code}
fit_2002 <- disaggregation::disag_model(dis_2002, iterations = 1000, field = TRUE,
iid = TRUE, family = "poisson", link = "log")
::: {.cell-output .cell-output-stderr} Fitting model. This may be slow. :::
``` {.r .cell-code} fit_2003 <- disaggregation::disag_model(dis_2003, iterations = 1000, field = TRUE, iid = TRUE, family = “poisson”, link = “log”)
::: {.cell-output .cell-output-stderr}
Fitting model. This may be slow.
:::
::::::
::::::: cell
``` {.r .cell-code}
summary(fit_2001)
Pooling the three years reduced uncertainty: the pooled model had smaller standard errors for both temperature and precipitation than any separate annual model. The temperature estimate was positive and the precipitation estimate negative, consistent with the values used to simulate the data.
Refrences
Nandi, A., Lucas, T., Arambepola, R., Python, A. (2023). disaggregation: An R Package for Bayesian Spatial Disaggregation Modeling. Journal of Statistical Software, 106(11). https://doi.org/10.18637/jss.v106.i11
Nandi, A., Lucas, T., Arambepola, R., Gething, P., Weiss, D. (2020). disaggregation: An R Package for Bayesian Spatial Disaggregation Modelling (original methods paper, Madagascar malaria worked example). https://arxiv.org/abs/2001.04847
GADM database of global administrative boundaries: https://gadm.org
Abatzoglou, J.T., Dobrowski, S.Z., Parks, S.A., Hegewisch, K.C. (2018). TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015. Scientific Data, 5, 170191. https://doi.org/10.1038/sdata.2017.191
These figures above show the resolution mismatch; on the left we have the low resolution public health case counts and on the right we have the high resolution environmental satellite covariates (Nandi et al. 2023, https://doi.org/10.18637/jss.v106.i11).
