What is the biggest open problem in statistics? I think this might be quite a common question. And I imagine quite a lot of people would answer with something that is quite specific to their research area. Ask someone who works on survival analysis and they might say the biggest problem is something about biased censoring. Ask someone who works on predictive modelling and they might say something about generalisability of prediction models.
For a while now I’ve had my answer in my mind. But I never got round to write anything about it. Partially because I am utterly unsure what the solution is. But the other day I decided that it was best just to get these ideas written down. At least then maybe I can move on until I come up with someway of moving the field forward in this area.
So, my answer to the question “what is the biggest open problem in statistics?” Is that almost every statistical analysis has too many medium sized problems that no-one can in practice can ever fix them all.
I’ll give a couple of examples to clarify a bit what I mean. In malaria mapping (one of the fields I work in) a typical analysis might involve trying to take some data on malaria prevalence, joining it to environmental variables, and fitting a predictive model such that we can predict the prevalence of malaria in other areas without data. Off the top of my head the medium sized problems here are:
the data might not be collected in a representative way
the predictor variables are measured with uncertainty
the true model is extremely complex (lags in space and time etc.)
we need to balance bias-variance
different diagnostic methods have different sensitivities and specificities
the locations of the individual surveys are uncertain
there are missing variables that are spatially correlated
the data aren’t iid
the estimated uncertainties will only account for parameter uncertainty, not model uncertainty
transporting the estimates in space or time is problematic
the predictors are highly correlated with each other
the model is complex and finding an MLE or MAP estimate might be computationally difficult
there may be small computational issues throughout (small number instability for example)
biased missing data
posterior or likelihood surface might be odd shaped
…
…
This was a 3 minute list, not something I have carefully curated.
Now, each of those problems is interesting and most of them have solutions. Uncertainty in predictor variables, there’s a method for that. None iid data due to spatial autocorrelation, there’s a method for that. The different sensitivities and specificities, there’s a method for that.
However, to account for all of these, broadly equally important factors, is close to impossible. For a number of reasons.
First, many of the methods for accounting for these things will be complex. Therefore, coding up the methods from scratch is difficult and likely to be error prone. But using expert built software doesn’t work because each R package (for example) only accounts for a small handful of the issues.
The more general the modelling ecosystem you are using is (and therefore the more factors you can account for) the more time it takes to code up all these issues. glm() in R is simple, and can account for a few of these factors. INLA can model a bunch of different things, but the model code might quickly become 10s of lines of difficult code. STAN is almost completely general, but a complex model might be hundreds of lines of complex code.
Aside from the coding, it is incredibly difficult for one analyst, especially anyone who isn’t extremely experienced, to know that all these issues exist and to know what the methods are and to understand the specifics well enough to use the methods appropriately. Even someone in their third year of PhD or early postdoctoral researcher, with perhaps 8 years of tertiary education, would never know all these things.
Even if you know everything and have the time to code them up, the solution to most or many problems use top use part of the data for something other than the core question. We can handle bias-variance with cross-validation (holding out some data). We can handle measurement error by estimating additional parameters about the uncertainty in the predictor variables. We can include additional model components that model the difference between different diagnostic tests. But these all take data and suck information away from our main question. This might be fine if we handle a few aspects. But when we want to handle 10 or 20 issues, we either need huge amounts of data or will have very limited data left for estimating the parameters we care about.
I don’t know if this next argument is true, so take with a pinch of salt. But I wonder if the data requirements for adding these additional components grows faster than linearly. Handling measurement error takes n1 ‘data points’ to estimate well, while accounting for differing diagnostic tests takes n2 ‘data points’. But if there are correlations here (the areas with high values of predictor variable 1 all have diagnostic test a) then this seems like to tease apart these affects might take more than n1 + n2 data points. I’ll need to think about this more though.
Next, there will always be some modelling decisions that have to be made by the analyst rather than driven by data. We often do sensitivity analysis to see whether these decisions mattered. Sensitivity on 1 variable is fine. Sensitivity on 10 variables either requires you to look at each variable separately (which is prone to miss issues) or do a latin hypercube in 10 dimensions of different values which is impossible.
Next, the analyst needs to understand and interpret the model output. As model complexity increases, this becomes difficult.
And finally, whatever model we fit and conclusions we make, the work needs to be explained, to patients or policy makers or peer reviewers, or someone. This becomes increasingly impossible with a model with 10 different components.
And so, to bring this together, it’s worth thinking about what this means and what to do. The first thing I am relatively confident about, is that as peer reviewers or question askers at conferences, we must accept that our pet issue is not more important than the other 10 issues and if the analyst has had a good go at handling a bunch of issues, we can’t just pile more and more issues on them to fix. It is impossible.
How we actually move beyond this as a problem though, I really do not know. Telling everyone that they can’t do any policy relevant analysis until they know everything about everything. We’ll never get anywhere. I know that developing software from the starting point of ‘let’s be extremely general and fix all problems’ is just very difficult. For some specific problems, like malaria mapping, the solution is to have big teams, 10 data analysts, all working on the same overall problem. This is good, but just doesn’t scale to every problem in applied science.
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).
