4 Geostatistical prediction
List of the main functions used in the chapter
| Function | R Package | Used for |
|---|---|---|
create_grid() |
RiskMap | Creating a regular grid of prediction locations over a study area |
predict() |
RiskMap | Predicting scientific targets at individual locations and across areas |
plot() |
RiskMap | Mapping predictive summaries at locations, across areas or both |
assess_extrapolation() |
RiskMap | Locating predictions that extrapolate in space or in the covariates |
create_convex_hull() |
RiskMap | Creating a study area from the convex hull of the sampled locations |
assess_prediction() |
RiskMap | Assessing the predictive performance of geostatistical models using cross-validation |
simulate_surface() |
RiskMap | Simulating surfaces and data from geostatistical models |
assess_simulation() |
RiskMap | Assessing the predictive performance of geostatistical models using simulations |
Geostatistics was originally developed and applied as a predictive tool for improving mining operations and planning. In the seminal paper on model-based geostatistics (MBG) (Diggle et al. 1998), the authors begin by stating:
“Conventional geostatistical methodology solves the problem of predicting the realized value of a linear functional of a Gaussian spatial stochastic process \(S(x)\) based on observations \(Y_i= S(x_i) + Z_i\) at sampling locations \(x_i\), where the \(Z_i\) are mutually independent, zero-mean Gaussian random variables.”
It is no coincidence that, although MBG represents a broad class of statistical models applicable to both estimation and prediction problems, the majority of its applications have focused on prediction. This is evident in many scientific fields where MBG has been used, including global public health, which is the primary application focus in this book. In line with the quote above, for most of the examples considered in this book, the analyses we illustrate aim to use a finite set of locations \(x_i\), typically corresponding to households or villages, to make inferences about an unobserved disease risk surface. Having performed estimation of MBG models in Chapter 3, we can now carry out spatial prediction to answer critical public health questions, such as identifying hot-spots and cold-spots (areas of unusually high or low disease risk, respectively) or determining whether the average prevalence across a region exceeds a predefined policy threshold.
In this chapter, we will explore these issues, beginning with a formal statistical formulation of the spatial prediction problem and a description of the necessary steps for its implementation. We will then demonstrate how spatial prediction is performed when considering predictive targets at individual locations and over areas. Finally, we will conclude the chapter by presenting methods for comparing the predictive performance of different models, with some insights into their relative strengths and weaknesses.
4.1 Spatial prediction using geostatistical models
Consider a generalised linear geostatistical model as defined in Section 1.5. The first step in spatial prediction is to define what we want to predict. We denote a location-specific predictive target by \(T(x)\), where \(x\) can be a sampled or an unsampled location in the study region \(A\). For the river blindness example in Liberia, our target is the underlying nodule prevalence surface, \[ T(x) = p(x) = \operatorname{logit}^{-1}\{d(x)^\top\beta + S(x)\}. \] This target describes the prevalence of palpable nodules in the population, the indicator measured by the Rapid Epidemiological Mapping of Onchocerciasis, rather than the result of a new survey with a particular sample size. One question is where nodule prevalence exceeds 20%, a threshold historically used to identify areas for mass drug administration (Amazigo 2008). We use it here to illustrate prediction of threshold exceedance.
The predictive distribution describes the possible values of the target conditional on the observed data. We write this as \([T(x)\mid y]\). By default, RiskMap treats the estimated regression coefficients and covariance parameters as fixed when making predictions. The resulting uncertainty describes the unobserved surface under the fitted model. It does not account for uncertainty in those parameter estimates. We use this default for the first Liberia predictions, then show how to propagate parameter uncertainty in Section 4.2.3.
With the parameters held fixed, \(T(x)\) is a deterministic transformation of \(S(x)\). We can therefore obtain samples from \([T(x)\mid y]\) by drawing from \([S(x)\mid y]\) and applying the target’s transformation to each draw. The resulting samples can be used to calculate a predictive mean, a prediction interval, or an exceedance probability such as \(P\{T(x)>0.2\mid y\}\). These summaries answer different questions: the mean provides a point prediction, an interval describes uncertainty, and an exceedance probability quantifies how strongly the fitted model and data support exceeding the specified threshold. The two steps are therefore to define the scientific target and then use its predictive distribution to answer the question of interest. For a mathematical derivation of the predictive distribution, including its analytical form for linear geostatistical models, see Section 4.6. Here we use Monte Carlo samples to provide a common approach for linear and generalised linear models.
We distinguish two types of target:
- Spatially continuous targets describe an unobserved surface, \(\mathcal{T}=\{T(x):x\in A\}\), such as prevalence across a country. Maps evaluate this surface at a finite set of prediction locations.
- Areal-level targets summarise that surface over a region. For example, the spatial average is \(\mathcal{T}_A=|A|^{-1}\int_A T(x)\,dx\), where \(|A|\) is the region’s area. This differs from a population-weighted average when population density varies across the region.
To approximate a continuous surface, we first choose prediction locations \(\tilde X=\{\tilde x_1,\ldots,\tilde x_q\}\), often the centres of a regular grid. We then generate predictive draws at those locations. The required dependence between locations depends on the question:
- Marginal prediction provides draws with the correct predictive distribution at each location, without preserving the full dependence between locations. It is sufficient for location-specific means, intervals and exceedance probabilities.
- Joint prediction preserves the predictive dependence between locations. Each draw represents one possible spatial surface and can be used to calculate quantities involving several locations together.
Joint prediction is generally more computationally demanding. For a map of location-specific summaries, marginal prediction is sufficient. For the distribution of a spatial average, total, maximum, or a simultaneous exceedance event, we need joint prediction. The expectation of a spatial average can be calculated from marginal means, but its variance also depends on the covariances between locations. Marginal draws must not be treated as complete spatial surfaces when calculating this uncertainty.
Why dependence matters for a regional prediction
Consider a small region divided into 16 equally sized cells. At each cell, suppose prevalence is either 10% or 30%, with equal probability. In an idealised example of perfect positive dependence, all cells are low together or high together. Alternatively, imagine drawing independently from the same distribution at each cell. Both constructions give the same predictive distribution at every cell and the same expected regional average of 20%, but very different uncertainty about that average (Figure 4.1).
Under perfect dependence, the regional average is either 10% or 30%, so its standard deviation is 10 percentage points. Under independence, let \(H\) be the number of high-prevalence cells. Then \(H\sim\operatorname{Binomial}(16,0.5)\) and the average is \(0.1+0.2H/16\), with a standard deviation of only 2.5 percentage points. To see what this means for a decision, suppose action is needed when the regional average exceeds 25%, a threshold above the expected average of 20%. Under perfect dependence, the probability of exceeding it is 50%. Under independence it is \(P(H\geq13)\approx1.06\%\). Treating dependent cells as independent would therefore make us much too certain that the threshold is not exceeded. Perfect dependence is used here to make the contrast clear. In a fitted spatial model the strength of dependence is estimated from the data.
4.1.1 Example: preparing predictions for Liberia (continuing from Section 3.3.1)
We now apply these ideas to the river blindness example. We first load the packages used in this chapter.
# Load the packages used in this chapter
library(sf)
library(terra)
library(geobounds)
library(elevatr)
library(dplyr)
library(ggplot2)
library(patchwork)
library(RiskMap)We then refit the model of Section 3.3.1, which uses log elevation and a spatial Gaussian process, with no nugget.
# Refit the river blindness model from Chapter 3
data(liberia)
liberia <- st_as_sf(liberia, coords = c("long", "lat"), crs = 4326)
fit_liberia <- glgpm(npos ~ log(elevation) + gp(),
denominator = ntest, data = liberia,
model_crs = 32629,
control_mcmc = set_control_mcmc(seed = 2026),
return_samples = TRUE,
family = "binomial")Next, we prepare the study boundary, a regular prediction grid and elevation values at the grid locations. The newdata argument of predict() takes these locations and covariates together in one spatial object, so their rows remain aligned. We add the covariates as new columns of the grid created by create_grid(), which keeps the record of its cells and study boundary that later calculations and maps use. For elevation we use the SRTM raster prepared in Chapter 2 (Figure 2.2), which is also the source of the elevations recorded for the surveyed villages. Predicting with a covariate from a different source than the one used to fit the model would introduce an avoidable inconsistency.
# Download the national boundary and project it to UTM zone 29N
liberia_adm0 <- gb_get_adm0("Liberia") |>
st_transform(crs = 32629)
# Create a 5 km grid, keeping cells that overlap the boundary
liberia_grid <- create_grid(liberia_adm0, spacing = 5,
distance_units = "km", keep = "overlap")
# Import the SRTM elevation raster and mask non-positive values
liberia_elev <- rast("data/lbr_elevation_srtm.tif")
liberia_elev[liberia_elev <= 0] <- NA
# Extract elevation at the grid locations, in the raster's CRS
extraction_points <- st_transform(liberia_grid, crs(liberia_elev))
liberia_grid$elevation <- extract(liberia_elev, extraction_points,
search_radius = 5000)[[2]]The prediction locations are spaced 5 km apart. A finer grid gives a finer representation of the surface, at greater computational cost. It does not add information to the survey data. With keep = "overlap", the grid also includes cells along the coast and the borders whose centres lie just outside the country. They matter when we add up population later in the chapter, and maps clip them to the boundary. The elevation raster has no values outside Liberia, so for these cells search_radius = 5000 takes the elevation of the nearest land cell within one grid spacing. We first set the few shoreline cells recorded at 0 m or below to missing, because the model uses the logarithm of elevation. In Chapter 2 we cautioned against borrowing a value across a coastline for a survey location. Here the value describes the part of a prediction cell that lies on land, which is the part these cells contribute to the population totals. The elevation raster and extraction points must use the same coordinate reference system.
We use marginal prediction because our immediate goal is to summarise prevalence separately at each grid location. The following call also retains the spatial draws so that we can inspect a predictive distribution. The summaries of prevalence are discussed in the next section. The arguments of predict() for a RiskMap model are documented in ?predict.RiskMap.
liberia_prediction <- predict(fit_liberia,
newdata = liberia_grid,
type = "marginal",
summaries = c("mean", "interval_width",
"cv", "exceedance"),
thresholds = 0.2,
keep_samples = TRUE,
seed = 2026,
messages = FALSE)Printing the result lists what it contains and how to map it.
liberia_predictionRiskMap prediction: 4077 locations
Response: prevalence (binomial model)
Draws: marginal
Quantities: response, link, spatial
Summaries: mean, interval_width, cv, exceedance
Exceedance thresholds: 0.2
Map with plot(x, summary = "mean").The predictions component contains the grid and its predictive summaries. Because we set keep_samples = TRUE, the spatial draws are retained in samples$spatial, with one row per prediction location and one column per draw. Using type = "joint" would preserve dependence between locations, but is unnecessary for the location-specific summaries shown here. Normally, keep_samples = FALSE avoids retaining these matrices after the summaries have been calculated.
Figure 4.2 shows what predict() does with this matrix. The spatial effect is drawn at every location (step 1), each target is calculated in every draw (step 2), and each row is then summarised across its draws (step 3). The last step, aggregation within areas, is used in Section 4.3.
For this non-Gaussian model, prediction reuses the conditional samples stored in the fitted object when they are available. We should check the conditional sampler as described in Section 3.5.2.3. Supplying a set_control_mcmc() object through control_mcmc requests a fresh conditional run. Otherwise a fresh run is needed only if the fitted object has no stored samples. The seed above controls random draws made during prediction.
To inspect uncertainty about the spatial process at the first grid location, we select the first row of the spatial-sample matrix. Its columns are repeated draws for that same location, rather than different locations in a single draw.
# Select the draws for the first grid location
spatial_effect <- liberia_prediction$samples$spatial[1, ]
spatial_draws <- data.frame(spatial_effect)
ggplot(spatial_draws, aes(x = spatial_effect)) +
geom_histogram(bins = 30, fill = "#347B87", colour = "white") +
labs(x = expression(S(tilde(x)[1])), y = "Number of draws") +
theme_minimal()
4.2 Spatially continuous targets
The previous section showed how to obtain predictive draws of the spatial effect at the grid locations. For a target \(T(x)=f\{d(x)^\top\beta+S(x)\}\), we apply \(f\) separately to each complete linear-predictor draw, including the fitted covariate contribution. With a nonlinear transformation, transforming the predictive mean is not the same as averaging the transformed draws. The latter gives the predictive mean of the target. By default, predict() takes \(f\) to be the inverse link function of the fitted model: the inverse logit for a binomial model, which gives prevalence, and the exponential for a Poisson model, which gives a rate. This response-scale target is always calculated, and its summaries are stored in columns whose names begin with response_. Other targets can be requested through the transform argument, a named list of functions. Each function receives the matrix of linear-predictor draws, with one row per location and one column per draw, and must return a matrix of the same size. For example, transform = list(odds = function(link) exp(link)) adds the predictive odds of infection, summarised in columns beginning with odds_ and mapped with plot(..., quantity = "odds"). Common targets are listed in Table 4.1.
| Predictive target \(T(x)\) | Name of the predictive target | Generalised linear model (GLM) family |
|---|---|---|
| \(d(x)^\top \beta + S(x)\) | Linear predictor | Any GLM |
| \(\frac{\exp\{d(x)^\top \beta + S(x)\}}{1+\exp\{d(x)^\top \beta + S(x)\}}\) | Prevalence | Binomial |
| \(\exp\{d(x)^\top \beta + S(x)\}\) | Rate per unit exposure | Poisson |
| \(S(x)\) | Spatial random effects | Any GLM |
| \(d(x)^\top \beta\) | Covariate effects | Any GLM |
These targets describe the underlying spatial surface and omit the nugget \(Z_i\). Whether a nugget belongs in a target depends on its scientific interpretation: measurement error is not part of the true surface, whereas genuine variation at a finer spatial scale may be relevant. For predictions at new locations, predict() excludes the nugget by default. include = "nugget" adds a new nugget realisation when the fitted model contains one. Adding a nugget can change both predictive uncertainty and the mean on a nonlinear response scale. A zero-mean nugget leaves the mean linear predictor unchanged, but generally does not leave the mean prevalence or rate unchanged. It should therefore be included because it belongs to the scientific target, rather than simply to widen an interval. The Liberia model used here has no nugget.
For a binomial model, the inverse link gives a probability, not a simulated positive count from a new survey. For a Poisson model, the inverse link gives a rate per unit exposure. The expected count additionally depends on exposure. The fixed covariate contribution and spatial effect can also be examined separately to understand their contributions to prediction.
4.2.1 Example: mapping river blindness prevalence in Liberia (continuing from Section 4.1.1)
For the fitted Liberia model, the prevalence target is \[
T(x)=\operatorname{logit}^{-1}\{\beta_0+\beta_1\log e(x)+S(x)\},
\tag{4.1}\] where \(e(x)\) is elevation in metres. The binomial inverse link is already part of fit_liberia, so the earlier predict() call has calculated the prevalence summaries without requiring us to supply the logistic transformation again. The requested summaries have the following interpretation:
"mean"averages the prevalence draws at each location."interval_width"is the width of the 95% prediction interval, in percentage points of prevalence."cv"is the standard deviation of the draws divided by their mean, a relative measure of uncertainty."exceedance", withthresholds = 0.2, gives the proportion of draws with prevalence above 20%, an estimate of \(P\{T(x)>0.2\mid y\}\).
The plot() function maps any of these summaries, and its options are documented in ?plot.RiskMap_prediction. Its default colours depend on what is mapped: predicted levels, uncertainty and exceedance probabilities each have their own palette. Several summaries can be mapped in one call. The result is built with the ggplot2 and patchwork packages, so we can add layers to every map with &, here the surveyed villages.
plot(liberia_prediction,
summary = c("mean", "interval_width", "exceedance")) &
geom_sf(data = liberia, shape = 21, size = 0.7, fill = "white",
colour = "black", stroke = 0.3, inherit.aes = FALSE)
The maps in Figure 4.4 separate predicted prevalence from the uncertainty about it. The interval width is read directly in percentage points of prevalence, which is how uncertainty is usually discussed in programme decisions. The coefficient of variation, also calculated above, is a relative measure: at low prevalence a small absolute uncertainty gives a large coefficient of variation, so it can make low-prevalence areas look more uncertain than they are in absolute terms. The exceedance map answers a different question from the mean map: it shows how strongly the fitted model and data support prevalence above 20% at each location, rather than simply whether the predicted mean exceeds 20%.
The prediction combines two components on the log-odds scale: the fixed effects, \(\beta_0 + \beta_1 \log e(x)\), and the spatial effect, \(S(x)\). Mapping them side by side shows how much of the geographical variation each explains. We centre the fixed effects on their average across the grid, because the intercept only shifts the whole surface, and the spatial effect already has mean zero. With a common colour scale, the two maps then compare variation with variation.
# Centre the fixed effects on their average across the grid
components <- liberia_prediction
components$predictions$fixed <- components$predictions$fixed -
mean(components$predictions$fixed)
# Use the same symmetric colour limits for both maps
component_range <- max(abs(c(components$predictions$fixed,
components$predictions$spatial_mean)))
component_limits <- c(-1, 1) * component_range
plot(components, quantity = c("fixed", "spatial"),
limits = component_limits)
In Figure 4.5, elevation produces a broad gradient from the coast inland. The spatial effect adds more local variation, with clusters of higher and lower prevalence than elevation alone would predict. Where the spatial effect is close to zero, either the data agree with the elevation trend or there are too few data nearby to suggest otherwise.
NoteHow many predictive draws are enough?
Summaries such as exceedance probabilities are estimated from a finite number of draws, \(B\). If the draws were independent, the Monte Carlo standard error of an estimated probability \(p\) would be \(\sqrt{p(1-p)/B}\), at most about 1.6 percentage points for the 1,000 draws used here. Draws produced by Markov chain Monte Carlo are correlated, so the effective number of draws is smaller and the error larger. This matters most when a location or area lies close to a classification cut-off. Repeating the prediction with a different seed, or with more draws requested through set_control_mcmc(), shows whether the conclusions are stable.
4.2.2 Where can we trust the map?
Prediction intervals do not by themselves warn us when a map extrapolates beyond the data. Two kinds of extrapolation matter, and they call for different responses.
The first is spatial. Far from the sampled locations, relative to the practical range of the estimated spatial correlation, the data say little about \(S(x)\) and predictions fall back on the covariate trend. For the Liberia model, with an exponential correlation function and scale parameter \(\hat\phi =\) 28.2 km, the practical range is \(\log(20)\hat\phi \approx\) 84 km. We can measure this by the share of the spatial variance at \(x\) that the data leave unexplained, \[ U(x) = 1 - \rho(x)^\top R^{-1} \rho(x), \tag{4.2}\] where \(\rho(x)\) holds the estimated correlations between \(x\) and the sampled locations, and \(R\) is the correlation matrix of the sampled locations. \(U(x)\) is close to 0 next to sampled locations and close to 1 where the data carry no information about \(S(x)\). Unlike the distance to the nearest sampled location, it accounts for the estimated scale of spatial correlation and for the clustering of the sampled locations. It depends only on where the data were collected, so for binomial and Poisson outcomes, which carry less information than an exact observation of \(S(x)\), it gives an optimistic picture.
The second kind concerns the covariates. Where a covariate lies outside the range observed in the data, predictions rely on extending the fitted relationship beyond the conditions from which it was estimated. The check should be made on the scale used by the model. The Liberia model uses the logarithm of elevation, so a few metres below the lowest surveyed village is a large step on the model’s scale. A combination of covariates can also be unusual even when each covariate is within its range, which the leverage of the location in the fixed-effects regression detects. Plug-in prediction intervals do not widen in either case, because they treat the regression coefficients as known.
The assess_extrapolation() function computes these indicators at each prediction location. It classifies each location as supported, spatial extrapolation (more than 80% of the spatial variance unexplained, by default), covariate extrapolation (a covariate outside its observed range, or a leverage larger than at any sampled location), or both.
liberia_extrapolation <- assess_extrapolation(fit_liberia,
newdata = liberia_grid)
table(liberia_extrapolation$extrapolation)
Supported Spatial Covariate
3581 385 108
Spatial and covariate
3 plot(liberia_extrapolation,
indicator = c("spatial_unexplained", "extrapolation")) &
geom_sf(data = liberia, shape = 21, size = 0.7, fill = "white",
colour = "black", stroke = 0.3, inherit.aes = FALSE)
In Figure 4.6, 10% of the grid locations have more than 80% of the spatial variance unexplained, in the gaps between clusters of surveyed villages. Elevation lies outside its observed range at 111 locations, 82 of them below the lowest surveyed village along the coast and the remaining 29 above the highest, in the northern highlands. Predictions in these areas should be reported with this caveat or supported with additional data: further surveys where the data are sparse, and caution about the form of the elevation effect where the covariate is extrapolated.
4.2.3 Accounting for uncertainty in the model parameters
The predictions above treat the estimated regression coefficients and covariance parameters as fixed. In practice, those estimates are also uncertain. For example, different values of the spatial variance \(\sigma^2\) and scale parameter \(\phi\) can imply different amounts of uncertainty about prevalence at an unsampled location. This additional source of uncertainty can matter particularly when the survey contains little information about those parameters.
RiskMap can propagate parameter uncertainty by setting parameter_uncertainty = TRUE in predict(). It draws parameter values from an approximate multivariate normal sampling distribution of the estimates, using their estimated covariance matrix. At each parameter draw, RiskMap generates new samples of the unobserved effects conditional on the data, calculates the target, and pools the target draws to obtain predictive summaries.
The approximation relies on the estimated covariance matrix and the normal approximation to the sampling distribution. These require care when parameters are poorly identified or near a boundary. Propagating uncertainty cannot compensate for an inadequately fitted model. The sampling-distribution perspective connects to the confidence intervals and parametric bootstrap discussed in Section 3.3.3.
We compare both approaches on the same Liberia grid, using the same fitted model, seed and total number of predictive draws. We request 95% prediction intervals as well as their widths. The second call uses 50 parameter draws. These Monte Carlo budgets control numerical precision and should be checked for stability in a substantive analysis.
prediction_summaries <- c("mean", "lower", "upper", "interval_width")
# Predict with the parameters fixed at their estimates
liberia_fixed <- predict(fit_liberia,
newdata = liberia_grid,
type = "marginal",
summaries = prediction_summaries,
level = 0.95,
seed = 2026,
messages = FALSE)
# Predict again, propagating uncertainty from 50 parameter draws
liberia_uncertain <- predict(fit_liberia,
newdata = liberia_grid,
type = "marginal",
summaries = prediction_summaries,
level = 0.95,
parameter_uncertainty = TRUE,
parameter_draws = 50,
seed = 2026,
messages = FALSE)We compare the two predictions location by location, grouping the locations by elevation, the covariate in the model.
# Interval widths in percentage points, one row per grid location
widths <- data.frame(
fixed = 100 * liberia_fixed$predictions$response_interval_width,
propagated = 100 *
liberia_uncertain$predictions$response_interval_width,
elevation = cut(liberia_grid$elevation,
breaks = c(0, 100, 300, Inf),
labels = c("Below 100 m", "100 to 300 m",
"Above 300 m"))
)
ggplot(widths, aes(x = fixed, y = propagated, colour = elevation)) +
geom_point(size = 0.7, alpha = 0.6) +
geom_abline(linetype = "dashed") +
scale_colour_manual(values = c("#00857F", "grey65", "#B5408F")) +
coord_equal() +
labs(x = "Fixed parameters (pp)", y = "Propagated uncertainty (pp)",
colour = "Elevation") +
theme_minimal() +
theme(legend.position = "top")
Propagating parameter uncertainty widens the interval at 67% of the locations, but the median width increases by only 3%. The changes are concentrated in the coastal lowlands (Figure 4.7). Below 100 m, intervals are typically 14% wider and the mean prevalence rises by a median of 2 percentage points, whereas above 300 m the predictions barely change. Uncertainty about the intercept and the elevation coefficient matters most where log elevation is furthest from that of the typical surveyed village, about 125 m. At low elevations it is far from that value, and the lowest locations lie below the observed range (Section 4.2.2). Plug-in prediction treats both coefficients as known, so its intervals miss this source of uncertainty. Propagating it partly corrects this, although it cannot tell us whether the log-linear form of the elevation effect still holds below the observed range.
The mean changes because the inverse link is nonlinear and the predictions now average over different parameter values. Some intervals become narrower. Part of this is Monte Carlo error from a finite number of draws, and part reflects that we are comparing a distribution calculated at one set of estimates with a mixture of distributions over different parameter values.
This choice is separate from marginal versus joint prediction. Marginal prediction remains sufficient for the location-specific comparison here. To quantify uncertainty about a regional target, we still need joint prediction, whether parameters are held fixed or their uncertainty is propagated.
4.2.4 Example: mapping malaria using age and elevation as predictors
We now consider an example on malaria mapping, where we use two different types of covariates: elevation, which is available at every location, and individual age, which is not. We use the malkenya_comm subset of the malkenya dataset (see Section 3.1.1.2) and, to reduce the computational burden, we take its first 1,000 rows.
# Select the community survey, as in Chapter 3
data(malkenya)
malkenya_comm <- malkenya[malkenya$Survey == "community", ]
# Take the first 1,000 individuals and project to UTM zone 36S
malkenya_comm1000 <- malkenya_comm[1:1000, ] |>
st_transform(crs = 32736)Based on the exploratory analysis shown in Section 3.1.1.2, we consider a geostatistical model for the individual binary outcomes, \(Y_{ij}\), the rapid diagnostic test (RDT) results for the \(j\)-th individual in the \(i\)-th household, with \(Y_{ij}=1\) for a positive test. Using a binomial geostatistical model, we model the probability \(p_{ij}\) of a positive RDT through the linear predictor \[ \begin{aligned} \log\left\{\frac{p_{ij}}{1-p_{ij}}\right\} &= \beta_0 + f(a_{ij}) + \beta_4 e(x_i) + S(x_i) + Z_i,\\ f(a) &= \beta_1 a + \beta_2 \max\{a - 15, 0\} + \beta_3 \max\{a - 40, 0\}. \end{aligned} \tag{4.3}\]
Here, \(e(x_{i})\) is the elevation in metres of the \(i\)-th household and \(a_{ij}\) is the age in years of its \(j\)-th individual. The age effect \(f(a)\) is a linear spline with knots at 15 and 40 years, following the exploratory analysis. The \(Z_i\) are independent zero-mean Gaussian variables with variance \(\tau^2\), a nugget shared by the individuals of the \(i\)-th household. It captures differences between nearby households that are not spatially structured, such as bed-net use or housing quality, which would otherwise be absorbed by \(S(x)\) at an artificially short scale. The nugget can be estimated because most households contribute several individuals. With a single binary outcome per location it could not be separated from the binomial variation.
# Linear splines in age with knots at 15 and 40 years, and a nugget
set.seed(2026)
fit_malkenya <- glgpm(RDT ~ Age + pmax(Age - 15, 0) +
pmax(Age - 40, 0) + elevation +
gp(nugget = TRUE),
data = malkenya_comm1000,
family = "binomial",
messages = FALSE)The estimated nugget variance, \(\hat\tau^2 =\) 0.86, is similar to the spatial variance, \(\hat\sigma^2 =\) 0.71. The scale parameter, \(\hat\phi =\) 4.4 km, is estimated imprecisely, with a 95% confidence interval from 0.7 to 27 km. The data show clear variation between households, but say little about how it is structured in space beyond a few kilometres.
When we use this model to predict RDT prevalence, age raises a question that elevation does not: it is not available as a geo-referenced covariate at every location. The answer depends on the research question. To describe risk at a particular age, we fix \(a_{ij}\) at that age for all prediction locations. To describe the general population, we need the prevalence averaged over the population’s age distribution.
Epidemiologists will recognise the second approach as direct standardisation. For each predictive draw and each location, we calculate prevalence for every age in a reference population and average those prevalences. With reference ages \(a_1, \ldots, a_n\), the age-standardised target is \[ T_{\textrm{std}}(x) = \frac{1}{n} \sum_{k=1}^{n} \operatorname{logit}^{-1}\{\beta_0 + f(a_k) + \beta_4 e(x) + S(x)\}, \] where \(f(a)\) is the spline function of age in Equation 4.3. The average is taken on the prevalence scale, so it is not the same as predicting at an average age. Here the reference population is the surveyed sample itself. This is reasonable if the sample represents the target population and the age distribution varies little across this small study area. A census or a gridded population age structure would provide a better reference when available.
We predict over the two sub-counties that contain most of the surveyed households, Kabondo Kasipul and Kasipul, retrieved with the geobounds package as in Chapter 2. Together they contain 327 of the 350 households in these data. The four neighbouring sub-counties contain only a few households each, near their edges, so predictions for them would rely mostly on spatial extrapolation (Section 4.2.2). The model is still fitted to all the data. We create a 500 m grid over the two sub-counties and extract elevation at the grid points.
# Select the two sub-counties and project them to UTM zone 36S
kenya_subcounties <- gb_get_adm2("Kenya") |>
filter(shapeName %in% c("Kabondo Kasipul", "Kasipul")) |>
st_transform(crs = 32736)
# Create a 500 m grid over the sub-counties
kenya_grid <- create_grid(kenya_subcounties, spacing = 0.5,
distance_units = "km", keep = "overlap")
# Elevation with a 1 km margin, extracted at the grid locations
kenya_elev <- get_elev_raster(locations = kenya_grid, z = 9,
clip = "bbox", expand = 1000) |>
rast()
kenya_points <- st_transform(kenya_grid, crs(kenya_elev))
kenya_grid$elevation <- extract(kenya_elev, kenya_points)[[2]]
# Fix age at 15 years for the first prediction
kenya_grid$Age <- 15The first prediction uses age 15 at every location. The second supplies the reference population through standardize, whose values replace the Age column of the grid. Neither includes the nugget, which predict() leaves out by default (Section 4.2), so the maps describe the underlying surface rather than prevalence in one particular household. Adding include = "nugget" would add a new draw of \(Z\) and widen the predictive intervals accordingly.
# Predict prevalence at age 15
kenya_age15 <- predict(fit_malkenya, newdata = kenya_grid,
summaries = "mean", seed = 2026,
messages = FALSE)
# Standardise to the ages of the surveyed individuals
reference_ages <- data.frame(Age = malkenya_comm1000$Age)
kenya_standardised <- predict(fit_malkenya, newdata = kenya_grid,
standardize = reference_ages,
summaries = "mean",
seed = 2026,
messages = FALSE)We map the two predictions with common colour limits, so that colours mean the same in both maps.
# Use common colour limits for the two maps
kenya_means <- c(kenya_age15$predictions$response_mean,
kenya_standardised$predictions$response_mean)
kenya_limits <- range(pretty(kenya_means))
age15_map <- plot(kenya_age15, limits = kenya_limits) +
labs(fill = "Mean prevalence, age 15")
standardised_map <- plot(kenya_standardised, limits = kenya_limits) +
labs(fill = "Mean age-standardised prevalence")
# Draw the two maps side by side
age15_map | standardised_map
The exploratory analysis of Section 3.1.1.2 showed that prevalence around the age of 15 is higher than at other ages. This is reflected in Figure 4.8, where prevalence at age 15 is higher everywhere than the age-standardised prevalence for the general population. The spatial patterns are very similar, because age and location enter the linear predictor additively on the log-odds scale.
4.3 Areal-level targets
When determining areal-level targets, we must first address the following questions:
- What spatially continuous target, \(T(x)\), are we seeking to aggregate?
- What are the spatial units, denoted as \(A_i\) for \(i = 1, \dots, N\), over which the aggregation is required?
- What aggregating function should be applied?
- Should spatial weights be incorporated in the aggregation, and if so, what weights are appropriate?
The first question was discussed in earlier sections, where examples of predictive targets that define \(T(x)\) can be found in Table 4.1. The answers to the remaining questions are context-dependent and closely linked to the primary research objectives. In public health settings, the areas \(A_i\) often correspond to administrative units where decisions are made and resources are allocated.
A common aggregating function for \(T(x)\) is the mean, \[\mathcal{T}_i = \frac{1}{|A_i|}\int_{A_i} T(x) \: dx. \tag{4.4}\]
We approximate the integral with the regular grid used for spatially continuous targets. Each grid point represents a square cell, and a cell on the edge of \(A_i\) may lie only partly inside it. Writing \(a_{ij}\) for the fraction of cell \(j\) that lies inside \(A_i\), \[\mathcal{T}_i \approx \frac{\sum_{j} a_{ij} T(\tilde{x}_j)}{\sum_{j} a_{ij}}. \tag{4.5}\] For grids created by create_grid(), predict() calculates these fractions exactly. Other aggregating functions follow the same idea. For example, the maximum of \(T(x)\) over \(A_i\) describes the worst-affected part of an area, and the median is less sensitive to a few extreme locations. A total of prevalence over an area has no meaning on its own, but weighting by population counts and summing gives the number of people affected, as we show below.
The formulation so far gives every location within an area equal weight. With spatial weights \(w(x) \geq 0\), the weighted mean is \[\mathcal{T}_i = \frac{\int_{A_i} w(x) T(x) \: dx}{\int_{A_i} w(x) \: dx}, \tag{4.6}\] approximated by replacing \(a_{ij}\) in Equation 4.5 with \(w(\tilde{x}_j)\,a_{ij}\). The most common choice in public health is population. A population-weighted mean prevalence describes the prevalence experienced by a typical resident, whereas the unweighted mean describes the typical location. Weighting by population counts and summing instead of averaging gives the expected number of people with the condition, the burden that planning and costing require. Other weights can reflect exposure or proximity to a source of infection, and binary weights restrict a target to part of an area.
These targets require joint prediction, and the order of operations matters (step 4 in Figure 4.2). The target is first calculated in every draw. For an areal target, locations are then aggregated within each joint draw, and only then are the resulting areal values summarised. For example, averaging the lower interval limits of the locations in an area does not give a lower interval limit for the area’s average.
4.3.1 Example: river blindness prevalence and burden in Liberia’s counties (continuing from Section 4.2.1)
We continue the Liberia analysis with the 15 counties as areas \(A_i\) and nodule prevalence, Equation 4.1, as \(T(x)\). A programme deciding where to deliver mass drug administration needs three things from this analysis: the average prevalence in each county, how certain we are that it exceeds 20%, and how many adults are affected.
The surveys examined adults aged 20 and over, so the prevalence we predict is the prevalence of nodules among adults, and population weights should count adults. We use WorldPop estimates of the number of people aged 20 and over in 2011, the year of the survey, at about 1 km. As in Section 2.6.1, we sum the counts into 5 km cells and use the cell centres as prediction locations, so that each location carries the number of adults it represents. This grid replaces the one built from the national boundary in Section 4.1.1, and we extract elevation at its locations in the same way.
# Download WorldPop counts of adults aged 20 and over in 2011, at 1 km
liberia_adults <- get_worldpop("LBR", year = 2011, ages = c(20, Inf),
resolution = "1km")
# Project to UTM 29N at 1 km, keeping totals, and sum to 5 km cells
adults_5km <- project(liberia_adults, "EPSG:32629", res = 1000,
method = "sum") |>
aggregate(fact = 5, fun = "sum", na.rm = TRUE)
names(adults_5km) <- "adults"
# Use the 5 km cells that overlap Liberia as prediction locations
county_grid <- create_grid(liberia_adm0, raster = adults_5km,
keep = "overlap")
county_grid$adults[is.na(county_grid$adults)] <- 0
# Adults in thousands, so that weighted sums are in thousands
county_grid$adults <- county_grid$adults / 1000
# Extract elevation at the new locations, as for the national grid
county_points <- st_transform(county_grid, crs(liberia_elev))
county_grid$elevation <- extract(liberia_elev, county_points,
search_radius = 5000)[[2]]
# Drop the few coastal cells beyond the extent of the elevation raster
county_grid <- county_grid[!is.na(county_grid$elevation), ]The first call below predicts the unweighted county means. The second weights each location by its number of adults, which gives the prevalence among the adults of each county and, by summing, the expected number of adults with nodules. The third sums over the whole country, because the interval for a national total cannot be obtained by adding the county intervals. Supplying boundaries makes predict() use joint draws automatically.
# Download the county boundaries
liberia_adm1 <- gb_get_adm1("Liberia")
# Average prevalence over the area of each county
liberia_county_area <- predict(fit_liberia,
newdata = county_grid,
boundaries = liberia_adm1,
aggregate = "mean",
summaries = "mean",
seed = 2026,
messages = FALSE)
# Adult prevalence and number of adults with nodules in each county
liberia_county_adults <- predict(fit_liberia,
newdata = county_grid,
boundaries = liberia_adm1,
aggregate = c("mean", "sum"),
weights = "adults",
summaries = c("mean", "lower",
"upper", "exceedance"),
thresholds = 0.2,
seed = 2026,
messages = FALSE)
# Number of adults with nodules in the whole country
liberia_national <- predict(fit_liberia,
newdata = county_grid,
boundaries = liberia_adm0,
aggregate = "sum",
weights = "adults",
summaries = c("mean", "lower", "upper"),
seed = 2026,
messages = FALSE)Mapping the county means next to the location-level prediction, with one colour scale, shows what aggregation does.
plot(liberia_county_area, spatial_scale = "both")
Aggregation hides the variation within counties that the left map shows: a county mean close to 20% can combine areas well above and well below that level, which matters when treatment is planned for smaller units.
The areal summaries are returned as a tidy table in boundary_summary, with one row per county, aggregation, quantity and summary. The to_table() function formats them for a report, with one row per county, the mean and its prediction interval in one column, and probabilities as percentages.
to_table(liberia_county_adults, label = "shapeName", sort_by = "mean",
col_names = c("County", "Adult prevalence (%)",
"P(> 20%) (%)",
"Adults with nodules (thousands)"))| County | Adult prevalence (%) | P(> 20%) (%) | Adults with nodules (thousands) |
|---|---|---|---|
| Gbarpolu | 23.9 (19.2 to 29.0) | 94.3 | 10.1 (8.09 to 12.2) |
| Lofa | 23.5 (18.8 to 29.0) | 92.3 | 29.4 (23.5 to 36.3) |
| Nimba | 21.2 (17.0 to 26.0) | 70.2 | 44.1 (35.3 to 54.1) |
| Bong | 21.0 (17.5 to 25.1) | 70.0 | 34.2 (28.4 to 40.7) |
| Grand Gedeh | 19.7 (15.3 to 24.7) | 41.8 | 12.6 (9.81 to 15.8) |
| River Gee | 18.8 (14.0 to 24.2) | 30.2 | 5.62 (4.19 to 7.24) |
| Margibi | 17.9 (13.6 to 23.0) | 17.2 | 21.9 (16.6 to 28.0) |
| Bomi | 14.4 (9.9 to 19.6) | 1.8 | 5.82 (4.01 to 7.93) |
| Grand Cape Mount | 14.0 (9.7 to 19.0) | 0.7 | 8.33 (5.78 to 11.3) |
| Grand Bassa | 12.4 (9.8 to 15.4) | 0.0 | 13.4 (10.6 to 16.7) |
| Sinoe | 12.1 (9.0 to 15.6) | 0.0 | 5.62 (4.18 to 7.26) |
| Montserrado | 10.6 (6.2 to 17.1) | 0.4 | 44.7 (26.3 to 72.3) |
| Rivercess | 9.8 (6.2 to 14.2) | 0.0 | 3.37 (2.15 to 4.90) |
| Maryland | 9.6 (6.7 to 13.2) | 0.0 | 5.28 (3.72 to 7.29) |
| Grand Kru | 8.9 (6.3 to 11.8) | 0.0 | 2.17 (1.54 to 2.89) |
Weighting by adults matters most where they are concentrated in part of a county. In Montserrado, which contains the capital, Monrovia, prevalence averaged over the county’s area is 18.0%, but among its adults it is only 10.6%, because most of them live in the city, where prevalence is low. The burden column ranks counties differently from prevalence, because a populous county with moderate prevalence can contain more affected adults than a sparsely populated county with high prevalence. Nationally, about 247,000 adults have nodules (95% prediction interval 219,000 to 280,000), out of 1,738,000 adults.
For a decision about mass drug administration, the exceedance probability is the most direct summary. A common way to act on it is to classify each county as likely above the threshold, likely below it, or uncertain. We map the probabilities next to the classes given by breaks, so that readers can see how far each county is from a cut-off.
# Exceedance probabilities on a continuous scale and in classes
probabilities <- plot(liberia_county_adults, aggregate = "mean",
summary = "exceedance")
classes <- plot(liberia_county_adults, aggregate = "mean",
summary = "exceedance", breaks = c(0, 0.2, 0.8, 1))
probabilities | classes
Only Gbarpolu and Lofa are likely above 20%, while Bong, Grand Gedeh, Nimba and River Gee are uncertain. The cut-offs of 20% and 80% are not statistical constants. They express how much risk of misclassification a programme accepts. Wrongly classifying a county as below the threshold leaves people untreated, whereas wrongly classifying it as above wastes resources, and the cut-offs can be chosen to reflect the relative costs of the two errors. Counties in the uncertain class are natural candidates for further surveys. Counties are used here to illustrate the method. Treatment decisions for onchocerciasis are usually made for smaller units, such as health districts, and the same code applies to them with their boundaries in place of liberia_adm1.
4.3.2 Example: Anopheles gambiae mosquitoes in Cameroon (continuing from Section 3.3.2)
For the analysis of Anopheles gambiae mosquitoes, the fitted Poisson model gives the expected number of mosquitoes caught by a trap, \[T(x) = \exp\{\beta_0 + \beta_1 e(x) + S(x)\},\] where \(e(x)\) is elevation in metres. The traps cover only part of the seven departments in which they lie, which together are about three times larger than the trapping area. Predicting over the departments would therefore rely largely on extrapolation. We instead use the convex hull of the trap locations as the study area \(A\).
Note
The convex hull of a set of spatial locations is the smallest convex polygon that encloses all of them. It can be pictured as a rubber band stretched around the outermost points. It is a useful study area when there is no natural boundary, or when administrative units extend far beyond the data.
The model is the one fitted in Section 3.3.2.
# Refit the Anopheles gambiae model from Chapter 3
data(anopheles)
an_fit <- glgpm(An.gambiae ~ elevation + gp(),
data = anopheles,
family = "poisson",
control_mcmc = set_control_mcmc(seed = 2026),
messages = FALSE)The traps were placed between 390 and 838 metres, and predicting at other elevations would extrapolate the fitted linear effect of elevation into conditions that may not suit the species (Section 4.2.2). We therefore use the binary weight \[w(x) = \begin{cases} 1 & \text{if elevation at } x \text{ is within the observed range,} \\ 0 & \text{otherwise,} \end{cases}\] and define the target as the weighted mean of Equation 4.6 over the study area, that is, the mean expected trap catch over its suitable part.
A total would be tempting, but adding the expected trap catch over grid cells does not give the number of mosquitoes in the area. A trap catch is not a density per unit area, so the sum would change with the grid spacing. The mean expected catch is a well-defined target that can be compared across areas and over time.
# Use the convex hull of the traps as the study area, with a 2 km grid
anopheles_hull <- create_convex_hull(an_fit$data)
anopheles_grid <- create_grid(anopheles_hull, spacing = 2,
distance_units = "km")
# Elevation with a 1 km margin, extracted at the grid locations
anopheles_elev <- get_elev_raster(locations = anopheles_grid, z = 9,
clip = "bbox", expand = 1000) |>
rast()
anopheles_points <- st_transform(anopheles_grid, crs(anopheles_elev))
anopheles_grid$elevation <- extract(anopheles_elev,
anopheles_points)[[2]]
# Weight 1 within the observed elevation range of the traps, 0 outside
observed_elevation <- range(an_fit$data$elevation)
anopheles_grid$suitable <- as.numeric(
anopheles_grid$elevation >= observed_elevation[1] &
anopheles_grid$elevation <= observed_elevation[2])
# Average the expected catch over the suitable part of the hull
anopheles_prediction <- predict(an_fit,
newdata = anopheles_grid,
boundaries = anopheles_hull,
aggregate = "mean",
weights = "suitable",
summaries = c("mean", "lower",
"upper"),
seed = 2026,
messages = FALSE)# Mark the traps and the locations excluded by the elevation weight
excluded <- anopheles_prediction$predictions |>
filter(suitable == 0)
plot(anopheles_prediction, spatial_scale = "location") +
geom_sf(data = excluded, shape = 4, size = 0.6, colour = "grey25",
inherit.aes = FALSE) +
geom_sf(data = an_fit$data, shape = 21, size = 0.8, fill = "white",
colour = "black", stroke = 0.3, inherit.aes = FALSE) +
labs(fill = "Mean expected catch per trap")
Of the 2158 grid locations, 96 lie outside the elevation range of the traps and receive zero weight, so the target excludes them. Figure 4.11 shows strong local variation in the expected catch, in line with the short scale of spatial correlation: the estimated scale parameter is 4.3 km, so the practical range is about 13 km. Over the suitable part of the study area, the predicted mean catch per trap is 9.1, with a 95% prediction interval from 7.7 to 10.7. This describes trap catches rather than the abundance of mosquitoes, which would also depend on the area each trap samples and on the trap’s efficiency.
4.4 Assessing the predictive performance of geostatistical models with cross-validation
In this section, we address the problem of identifying the geostatistical model that offers the best predictive performance among a set of candidates. However, we first need a clear definition of predictive performance which can be used to select suitable statistical tools that can be to used to evaluate it. Broadly speaking, predictive performance is defined by how well a model’s predictive distribution aligns with observed data. Evaluating predictive performance requires examining two key characteristics of the predictive distribution: sharpness and calibration.
4.4.1 How to split geostatistical data for model performance comparisons
To carry out the assessment and comparison of the predictive performance of geostatistical models using the methods illustrated in the following sections, the first crucial step is to decide which approach to use to split the dataset into a training set and test set. In the context of geostatistical analysis, the training set is the subset of the original dataset that is used to estimate the model parameters and then predict at the locations of the test set, which is used to assess how well the model can predict unseen data.
When splitting geostatistical datasets for model performance evaluation, it is essential to consider the objective of the spatial prediction. Here, we consider two main prediction objectives: 1) predicting in areas disjoint from the study area where there are no data; 2) inferring the spatial surface of disease risk within a given study area. Let us give an example to better explain the difference between these two objectives. Under objective 1), if survey data exist for a specific country, say Kenya, but not for a neighboring country, say Somalia, a model trained on Kenyan data might then be used to predict disease prevalence in Somalia. In objective 2), instead, high-resolution risk maps of a health outcome might be required within a single country. The splitting of geostatistical data for predictive performance assessment should align with the intended application of the model and reflect whether the goal is to predict in entirely disjoint areas with no data (objective 1) or to infer the spatial surface of interest within the same study area (objective 2). Below, we elaborate on suitable methods for each of these two.
In the first scenario (objective 1) ), splitting the data using k-fold cross-validation with spatially coherent folds is more appropriate. These folds can be created using clustering methods such as k-means or hierarchical clustering (see Yin et al. (2026) for review of the main clustering methods). These methods group spatial locations based on their coordinates or covariates, ensuring each fold represents a distinct spatial region. This approach mimics the challenge of predicting in disjoint areas by withholding entire spatial clusters from the training data. For instance, k-means minimizes within-cluster variance to produce compact, non-overlapping groups, while hierarchical clustering organizes data into nested clusters based on a linkage criterion like geographic distance. In the code below we show an example of the use of these methods for data splitting, using the spatial_clustering_cv() function from the spatialsample package (Mahoney et al. 2023).
library(spatialsample)
set.seed(123)
liberia_sf <- st_transform(liberia, crs = 32629)
# K-means
kmeans_liberia <- spatial_clustering_cv(liberia_sf,
v = 10, cluster_function = "kmeans")
kmeans_plot <- autoplot(kmeans_liberia) + ggtitle("K-means")
# Hierarchical clustering
hclust_liberia <- spatial_clustering_cv(liberia_sf,
v = 10, cluster_function = "hclust")
hclust_plot <- autoplot(hclust_liberia) +
ggtitle("Hierarchical clustering")
library(gridExtra)
grid.arrange(kmeans_plot, hclust_plot, ncol = 2)
The results of the splitting are shown in Figure 4.12, where we have generated 10 folds using the k-means and hierarchical clustering methods with the default options. We can observe that the test sets are very similar, with fold 3 from the k-means and fold 7 from hierarchical clustering being identical. If using only the distance between locations to cluster the data into different folds, the choice of clustering method should have much impact on the model assessment as long as we repeat the split and assessment of the methods multiple times as we will show later in the next sections.
Before proceeding further, it is important to highlight that when the objective is predicting an outcome in entirely disjoint areas, the risks and limitations of this assessment should be carefully considered. In such cases, predictions are predominantly driven by covariates because the spatial Gaussian process \(S(x)\) cannot extrapolate spatial patterns across regions with no connecting data. When the disjoint region requiring predictions is far from the sampled data, the Gaussian process contributes primarily to inflating prediction uncertainty and has no tangible impact on the point predictions, reflecting the lack of ground data in the area. Consequently, this validation approach primarily assesses the predictive power of the covariates rather than the geostatistical model as a whole. This limitation may lead to overestimating the model’s generalization capabilities, especially if the model is being used to predict areas that are far away from the study area. Hence, the results of this validation should be interpreted cautiously, emphasizing the critical role of covariates in driving predictions under such conditions.
We now consider the second prediction objective of inferring the spatial surface within the study area. In this case, a random sampling approach can be employed. However, to preserve a good spatial coverage of the study area, a minimum distance can be imposed between selected locations. For example, methods such as spatial thinning or stratified sampling can be adapted to enforce a minimum spatial separation. This can be done in R using the subsample.distance() function from the spatialEco package (Evans and Murphy 2021).
library(spatialEco)
set.seed(123)
# Regulirized spatial sampling: 30km
dist30_liberia <- subsample.distance(liberia_sf,
size = 20, d = 30000)
dist30_plot <- ggplot(dist30_liberia) + geom_sf() +
theme_minimal() +
ggtitle("Minimum distance: 30km")
# Regulirized spatial sampling: 40km
dist40_liberia <- subsample.distance(liberia_sf,
size = 20, d = 40000)
dist40_plot <- ggplot(dist40_liberia) + geom_sf() +
theme_minimal() +
ggtitle("Minimum distance: 40km")
grid.arrange(dist30_plot, dist40_plot, ncol = 2)
In the code above, we randomly selected 20 locations from the original dataset, enforcing minimum distances of 30 km and 40 km between points. The results are presented in Figure 4.13. The choice of the minimum distance for the test sample should ideally consider the scale parameter \(\phi\) of the spatial Gaussian process \(S(X)\). To minimize correlation within the test set, the minimum distance could be set to a value greater than \(\phi\). However, increasing the minimum distance makes it more challenging to obtain a test set of the desired sample size. Reducing correlation within the test set is important for ensuring that the assessment of predictive performance is as robust as possible.
A key challenge in both approaches is the potential lack of independence between the training and test datasets due to spatial correlation. This can lead to an overly optimistic evaluation of model performance, as the training set may already contain information that is correlated with the test set. However, if we assume that all models benefit equally from this correlation during performance evaluation, the methods we illustrate next can still provide a reliable comparison of model performance in relative terms, even if the considered performance metrics may be somewhat biased.
4.4.2 Assessing calibration using the nonrandomized probability integral transform
A model is considered well-calibrated if its predictions adequately reflect the true uncertainty of the data. Assessment of a model’s calibration can be carried out in several ways. We can examine the agreement between the point predictions, say \(\hat{y}_i\), and the observed outcomes \(y_i\), commonly referred to as accuracy. The mean squared error (MSE), defined as \(\text{MSE} = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2\), is an example of a commonly used metric to evaluate the accuracy of a model. To characterize calibration more fully, it is also essential to quantify the spread of the predictive distributions, which also defines the precision of predictions. For example, prediction intervals generated from the predictive distribution for \(y_i\) can help assess the precision of a geostatistical model. In summary, a well-calibrated model is both accurate and precise.
Assessment of the calibration of the model should precede the assessment of sharpness, since, as we shall see in the next section, this relies on the model being well-calibrated. The probability integral transform (PIT) was originally proposed by Dawid (1984) as a way to assess the calibration of a model. The PIT is based on the simple observation that if we consider a variable \(Y\) and apply the transformation \(Y^* = F_{Y}(Y)\), where \(F_{Y}(\cdot)\) is the cumulative density function (or cumulative distribution if \(Y\) is discrete) of \(Y\), it then follows that \(Y^*\) follows a uniform distribution in the unit interval. The fundamental problem and the reason why statistical modelling exists is that we do not \(F_{Y}\) but we would like to propose a model \(M\) that we believe adequately approximates \(F_{Y}\) with \(F_{M}\). The main advantage of using the PIT is that assessment of whether \(M\) is well calibrated reduces to assessing whether the transform set of data \(F_{M}(y_1), \ldots, F_{M}(y_n)\) follows a uniform distribution. To illustrate this, let us a consider a simple simulated example.
We first create a function that can compute the PIT.
# Define a function for the PIT
# F_M is the cumulative density function generated by the adopted model
pit_transform <- function(data, F_M) {
# Apply the model CDF to the data
pit_values <- F_M(data)
return(pit_values)
}We then apply this function to show how the assessment of calibration through the PIT is carried out under two scenarios: 1) \(F_{Y}\) is a Gamma distribution with shape 2 and scale parameter 1; 2) \(F_{Y}\) is Student’s T distribution with 3 degrees of freedom. In each of the two scenarios, we show how the PIT behaves when our model coincides with thetrue model (i.e. \(F_{M} = F_{Y}\)) and when instead our model \(M\) is a Gaussian distribution, with mean and variance estimated from the data.
# Generate data from a skewed (gamma) distribution
n <- 1000
shape <- 2
rate <- 1
data_gamma <- rgamma(n, shape = shape, rate = rate)
# Define correct and incorrect model CDFs for gamma data
# Correct CDF
model_cdf_correct_gamma <- function(x) pgamma(x, shape = shape, rate = rate)
# Incorrect CDF: Assume data is normal (wrong model)
model_cdf_incorrect_gamma <- function(x) pnorm(x, mean = mean(data_gamma), sd = sd(data_gamma))
# Compute PIT values for gamma data under both models
pit_correct_gamma <- pit_transform(data_gamma, model_cdf_correct_gamma)
pit_incorrect_gamma <- pit_transform(data_gamma, model_cdf_incorrect_gamma)
# Generate data from a heavy-tailed (t) distribution
df <- 3
data_t <- rt(n, df = df)
# Define correct and incorrect model CDFs for t data
# Correct CDF
model_cdf_correct_t <- function(x) pt(x, df = df)
# Incorrect CDF: Assume data is normal (wrong model)
model_cdf_incorrect_t <- function(x) pnorm(x, mean = mean(data_t), sd = sd(data_t))
# Compute PIT values for t data under both models
pit_correct_t <- pit_transform(data_t, model_cdf_correct_t)
pit_incorrect_t <- pit_transform(data_t, model_cdf_incorrect_t)
# Combine results into a data frame for plotting
df_plot <- data.frame(
PIT = c(pit_correct_gamma, pit_incorrect_gamma, pit_correct_t, pit_incorrect_t),
Model = rep(c("Correct Gamma Model", "Incorrect Normal Model for Gamma",
"Correct t Model", "Incorrect Normal Model for t"), each = n)
)
# Plot PIT distributions
hist_plot <- ggplot(df_plot, aes(x = PIT)) +
geom_histogram(aes(y = ..density..), bins = 20, fill = "skyblue", color = "black", alpha = 0.7) +
geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
facet_wrap(~ Model, scales = "free") +
labs(title = "Probability Integral Transform (PIT) Distributions",
x = "PIT Values",
y = "Density") +
theme_minimal()
# Display the plot
print(hist_plot)
In Figure 4.14, we assess whether the PIT transformed data follow a uniform distribution if the bars of the histogram are all approximately at the same height as indicated by the dashed line. Our preference to this diagnostic plot, is to use a qq-plot for a uniform distribution as shown below.
# QQ plot to check uniformity of PIT values
qq_plot <- ggplot(df_plot, aes(sample = PIT)) +
stat_qq(distribution = qunif) +
stat_qq_line(distribution = qunif, color = "red") +
facet_wrap(~ Model, scales = "free") +
labs(title = "QQ Plot of PIT Values against Uniform Distribution",
x = "Theoretical Quantiles (Uniform)",
y = "Sample Quantiles (PIT Values)") +
theme_minimal()
# Display both plot
print(qq_plot)
In the results of Figure 4.15, we consider a model well-calibrated if the qq-plot is as close as possible to the identity line. It is clear that in the first case the use of Gaussian distribution is violated by the skewness of the data generated from a Gamma distribution, and in the second scenario by the heavier tail of the Student’s T distribution.
In the simulated example, the PIT might seem unnecessary because we can directly assess if the data follow the assumed distribution (Gamma, Gaussian, or Student’s T) by comparing empirical and theoretical distributions. In these simple models, including linear Gaussian geostatistical models, the marginal distribution of the data is analytically available, so assessing goodness-of-fit is straightforward.
However, in the case of Binomial and Poisson geostatistical models, the marginal distribution of the data \(Y\) – unconditioned on the spatial random effects \(S(x)\) – cannot be obtained analytically. Specifically, geostatistical models for count data assume that \(Y\), conditioned on \(S(x)\), follows either a Binomial or Poisson distribution. But once we integrate out \(S(x)\) to obtain the marginal (i.e., unconditioned) distribution of \(Y\), this distribution is no longer Poisson or Binomial. Instead, it becomes an overdispersed count distribution that lacks a closed-form expression.
In such cases, the PIT becomes valuable for evaluating model adequacy, as it allows for assessing the uniformity of PIT values against the assumed model distribution, even when the marginal distribution is intractable. However, direct application of the PIT raises the issue that thetransformation \(F_{Y}(Y)\) does not follow a uniform distribution, because of the discrete nature of the random variable \(Y\). Hence, for a given discrete-outcome model \(M\), we use a modified version of the PIT, referred to as nonrandomized PIT (henceforth, nPIT), originally proposed by Czado et al. (2009), taking the form \[\text{nPIT}(u,y) = \begin{cases} 0, & u \leq F_{M}(y-1) \\ \frac{u - F_{M}(y-1)}{F_{M}(y-1) - F_{M}(y-1)}, & F_{M}(y-1) \leq u \leq F_{M}(y) \\ 1, & u \geq F_{M}(y) \end{cases}. \tag{4.7}\]
We then take the average nPIT over the observed outcomes \(y_{1}, \ldots, y_{n}\), hence \[\text{AnPIT}(u) = \frac{1}{n} \sum_{i=1}^n \text{nPIT}(u, y_{i}). \] Assessment of calibration is then carried out by checking wether the AnPIT is as close as possible to the identity function, i.e. \(\text{AnPIT}(u) = u\). Hence, by plotting \(\text{AnPIT}(u)\) against \(u\), we can consider any deviations from the identity line as evidence that the model is not well-calibrated.
We now illustrate the use of the AnPIT through a simulated example. We simulate data from a negative Binomial distribution, with mean \(\lambda = 5\) and dispersion parameter \(\alpha = 1/2\), hence \(E[Y] = \lambda\) and \(\text{Var}[Y] = \lambda \times (1 + \alpha \lambda)\).
# Parameters for the Negative Binomial distribution
set.seed(123)
n <- 200 # number of simulations
lambda <- 5 # mean
dispersion <- 0.5 # dispersion parameter
size <- 1 / dispersion # size parameter for negative binomial
# Simulate data from the Negative Binomial distribution
sim_data_nb <- rnbinom(n, size = size, mu = lambda)
# Define the nonrandomized PIT function
compute_npit <- function(y, u, cdf_func) {
# Calculate cumulative probabilities F_M(y-1) and F_M(y) using the provided CDF function
f_y_minus_1 <- cdf_func(y - 1)
f_y <- cdf_func(y)
# Apply the piecewise formula for nPIT
if (u <= f_y_minus_1) {
return(0)
} else if (u <= f_y) {
return((u - f_y_minus_1) / (f_y - f_y_minus_1))
} else {
return(1)
}
}
# Generate a sequence of u values from 0 to 1
u_values <- seq(0, 1, length.out = 100)
# Calculate the average nPIT for each u value for both models
AnPIT_nb <- sapply(u_values, function(u) {
mean(sapply(sim_data_nb, compute_npit, u = u, cdf_func = function(k) pnbinom(k, size = size, mu = lambda)))
})
AnPIT_poisson <- sapply(u_values, function(u) {
mean(sapply(sim_data_nb, compute_npit, u = u, cdf_func = function(k) ppois(k, lambda = lambda)))
})
# Combine the results into a data frame for plotting
AnPIT_data <- data.frame(
u = rep(u_values, times = 2),
AnPIT = c(AnPIT_nb, AnPIT_poisson),
Model = rep(c("Negative Binomial (correct model)", "Poisson (wrong model)"), each = length(u_values))
)
# Plotting the average nPIT as a function of u for both models
ggplot(AnPIT_data, aes(x = u, y = AnPIT, color = Model)) +
geom_line() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") + # Identity line
facet_wrap(~ Model) +
coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
xlab("u") +
ylab("AnPIT") +
theme_minimal() +
theme(legend.position = "none") 
Figure 4.16 shows clear evidence that the Poisson model is not adequately accounting for the overdispersion of the simulated data from the Negative Binomial distribution.
How can we adapt this approach to generalized linear geostatistical models (GLGMs)? The concept is straightforward, though the implementation is less so (but fortunately, we have already done this in RiskMap). Essentially, in all relevant equations above, we need to replace \(M\) with the predictive distribution of the fitted geostatistical model.
To explain further, suppose we have fitted our model using data \(Y^{(1)} = \left( y^{(1)}_{1}, \ldots, y^{(1)}_{n} \right)\) (the training dataset), and we want to use another set of data \(Y^{(2)} = \left( y^{(2)}_{1}, \ldots, y^{(2)}_{m} \right)\) (the test dataset) to assess calibration; the selection of \(Y^{(1)}\) and \(Y^{(2)}\) can be done using the methods described in Section 4.4.1. Our model \(M\) in the equations above corresponds to the distribution of \(\left[ Y^{(2)} \mid Y^{(1)} \right]\) – the distribution of \(Y^{(2)}\) conditioned to \(Y^{(1)}\) derived from our geostatistical model. Since its derivation cannot be done analytically, we use Monte Carlo methods to approximate \(\left[ Y^{(1)} \mid Y^{(2)} \right]\). We discuss this in more detail in Section 4.6. However, the key point here is to understand why we perform this assessment and how to interpret the results. We next show the application of the AnPIT diagnostic to the example on river blindness mapping in Liberia.
4.4.2.1 Example: assessing the calibration of two geostatistical models for river blindness mapping
In our example on river blindness mapping, we now consider two models. A mode, which we call \(M_{0}\), which is an intercept only geostatistical model, with linear predictor \[M_{0} \: : \: \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_0 + S(x_i) ,\] and a geostatistical model, which we denote as \(M_{1}\), which uses a linear spline with a not in 150 meters, and has linear predictor \[M_{1} \: : \: \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_0 + \beta_{1}e(x_i) + \beta_{2}\max\{e(x_i) - 150, 0\} + S(x_i).\] The reason for using this type of linear spline in \(M_{1}\) is explained in Section 3.1.1.
Hence, we first fit the two models.
liberia <- st_as_sf(liberia, coords = c("long", "lat"), crs = 4326)
# Fitting M0
set.seed(123)
M0_fit <- glgpm(npos ~ gp(), denominator = ntest, data = liberia,
model_crs = 32629,
family = "binomial", messages = FALSE)
# Fitting M1
M1_fit <- glgpm(npos ~ elevation + pmax(elevation - 150, 0) + gp(),
denominator = ntest, data = liberia,
model_crs = 32629,
family = "binomial", messages = FALSE)Since the main goal of this analysis is to draw predictive inference on river blindness prevalence within Liberia, it would be more appropriate to use a regulirized subsampling scheme to split the data. To generate the AnPIT diagnostic plot in R, we can use the assess_prediction() function in RiskMap
regularized <- assess_prediction(list(M0 = M0_fit, M1 = M1_fit),
metrics = "PIT",
method = "regularized", min_dist = 20,
size = 9, iter = 10, messages = FALSE)p1 <- plot(regularized, metric = "PIT", combine_pit = TRUE)
p2 <- plot(regularized, metric = "PIT", pit_mode = "all")
grid.arrange(p1, p2, nrow = 2, heights = c(1, 1))
In the above function, we have used a spatially regularized sampling scheme (method = "clusterized"), where the locations are at least 20 kilometers (min_dist = 20) apart from each other to create 10 test sets (iter = 10) each of size 9 (size = 9) which corresponds to 10% of the original dataset. The results reported in Figure 4.17 and Figure 4.18 show both the AnPIT for each of the 10 test sets and the averaged AnPIT across the 10 test sets, for both \(M_0\) and \(M_1\). The AnPIT plots show curves that are not too far from the identity line, which leads us to conclude that both \(M_0\) and \(M_1\) can be considered to be approximately well calibrated models.
Let us how if these results would change if we were to split the data using a clustering algorithm approach. We then use the assess_prediction(), to create 10 folds and set method = "cluster".
cluster <- assess_prediction(list(M0 = M0_fit, M1 = M1_fit),
metrics = "PIT",
method = "cluster", fold = 10,
messages = FALSE)p1 <- plot(cluster, metric = "PIT", combine_pit = TRUE)
p2 <- plot(cluster, metric = "PIT", pit_mode = "all")
grid.arrange(p1, p2, nrow = 2, heights = c(1, 1))
The AnPIT results in Figure 4.19 and Figure 4.20 are very similar to those in Figure 4.17 and Figure 4.18, hence similar conclusions are drawn.
4.4.3 Assessing calibration and sharpness using continuous ranked probability scores
We now consider another important aspect that helps to describe the predictive performance of a geostatistical model, namely sharpness. As the name suggests, sharpness refers to the spread of the predictive distribution. A sharp predictive distribution has a narrow spread, indicating that the model is highly confident in its predictions. Importantly, sharpness does not depend on how close predictions are to the true values (i.e., it is independent of accuracy). A model with high sharpness will produce predictions with low variance, but this sharpness is only meaningful if the model is well-calibrated, ensuring that the high confidence aligns with observed values. For example, a geostatistical model for disease prevalence mapping that predicts a very tight range of possible prevalence values for a location is producing sharp predictions. However, if the actual prevalence often falls outside this range, then the predictions may be sharp but poorly calibrated, thus limiting the value of quantifying sharpness alone.
To assess predictive performance while considering both sharpness and calibration, we introduce the continuous ranked probability score (CRPS) (Gneiting et al. 2007) and its scaled variant, the scaled CRPS (SCRPS) (Bolin and Wallin 2023). The CRPS evaluates the agreement between the predicted cumulative distribution function of the model \(M\), denoted as in the previous section by \(F_M\), and the observed value \(y_i\). It is formally defined as
\[\text{CRPS}(F_M, y_i) = \int_{-\infty}^\infty \left[ F_M(y) - \mathbf{1}(y \geq y_i) \right]^2 dy,\]
penalizing deviations of the predictive distribution from the observed value. The CRPS balances sharpness and calibration, with smaller values indicating better predictive performance.
The scaled CRPS (SCRPS), proposed by Bolin and Wallin (2023), refines the CRPS by addressing its dependence on the scale of the outcome variable, which can vary across observations. The SCRPS is defined as
\[\text{SCRPS}(F_M, y_i) = -\frac{1}{2} \left[ 1 + \frac{\text{CRPS}(F_M, y_i)}{\text{CRPS}(F_M, F_M)} + \log\left\{2 \, \text{CRPS}(F_M, F_M)\right\} \right]. \tag{4.8}\]
where \(\text{CRPS}(F_M, F_M)\) denotes the expected CRPS when the observation is drawn from the predictive distribution itself, i.e., \(Y \sim F_M\).
In settings where the predictive distribution is discrete, such as the Binomial and Poisson geostatistical models, a discrete analogue to CRPS is used. Let \(F_M(k)\) be the cumulative distribution function over discrete support \(\{0, 1, 2, \dots\}\). The discrete CRPS is given by
\[\text{CRPS}(F_M, y_i) = \sum_{k \in \mathcal{S}}^{\infty} \left[ F_M(k) - \mathbf{1}(k \geq y_i) \right]^2,\]
where \(\mathcal{S}\) represents the support of the outcome \(Y_{i}\), corresponding to \(\{0,1,\ldots,m_i\}\) for the Binomial distribution, and the set of non-negative integers as in a Poisson distribution. This version remains proper and interpretable as a measure of calibration and sharpness for integer-valued outcomes. Analogously, a scaled version of the discrete CRPS can be constructed using the same transformation as in the continuous case, by replcing the continuous CRPS with its discrete counterpart in Equation 4.8.
providing a scale-invariant and proper scoring rule for discrete predictive models. These discrete versions are particularly useful in count data models or models for ordinal responses.
The expressions for the CRPS and SCRPS may initially appear complex, but their appropriate application and interpretation lie in understanding their differences and the contexts in which each is more suitable. These differences become more intuitive when considering the more familiar concept of the standard error. The standard error of a predictive distribution represents its uncertainty, but it is not inherently standardized. In many cases, there exists a mean-variance relationship, where the variance increases with the mean. For instance, in Poisson-distributed data, higher mean values correspond to larger variances. This lack of standardization means that directly comparing predictive uncertainties across different scales can be misleading. Similarly, the CRPS is sensitive to the scale of the predictive distribution, as it penalizes errors in absolute terms, often giving disproportionate weight to observations with larger predictive uncertainty.
The SCRPS addresses this issue by standardizing predictive errors relative to the expected spread of the predictive distribution, ensuring local scale invariance. This means that the penalty for an incorrect prediction is adjusted proportionally to the scale of the predictive distribution. As a result, the SCRPS enables meaningful comparisons across observations with varying levels of uncertainty. To draw an analogy, the SCRPS standardizes the CRPS in much the same way that the coefficient of variation (CV) standardizes the standard error. By dividing the standard error by the mean of the predictive distribution, the CV expresses variability relative to the mean. This property is particularly valuable in geostatistical applications, where uncertainty often varies spatially due to factors such as non-stationarity of the predictive distribution or irregularly spaced observations.
However, while the SCRPS is more robust in settings where variability differs significantly across observations, there are contexts where the CRPS may be preferred. Specifically, in cases where greater emphasis should be placed on accurately predicting observations with larger uncertainty, the CRPS is more appropriate because it assigns greater weight to predictions with higher variance. For example, as highlighted in Bolin and Wallin (2023), when forecasting rare but impactful events (e.g., extreme weather conditions), the larger uncertainties associated with these events may naturally demand greater attention. In such scenarios, the CRPS’s sensitivity to scale can be an advantage, as it prioritizes minimizing errors for these high-variance predictions.
When evaluating the predictive performance of a model \(M\), we proceed as follows. Let \(\left(x_{1}^{(i)}, \ldots, x_{m_i}^{(i)}\right)\) represent the set of locations in the \(i\)-th test set, for \(i=1, \ldots, T\). For each location, we compute the score (either CRPS or SCRPS) and denote it as \(r_{M}\left(x_{j}^{(i)}\right)\).
To assess the overall performance of model \(M\), we calculate a weighted average of the scores across all test sets, where the weights are proportional to the sample size of each test set. The final score assigned to model \(M\) is given by: \[\frac{1}{\sum_{i=1}^T m_i} \sum_{i=1}^{T} \sum_{j=1}^{m_i} r_{M}\left(x_{j}^{(i)}\right). \tag{4.9}\]
Models assumed to be well-calibrated and achieving a lower score based on the above expression are considered to provide sharper predictions and, consequently, better predictive performance.
4.4.3.1 Example: mapping river blindness in Liberia (continuing from Section 4.4.2.1)
In Section 4.4.2.1, we have assessed the calibration of models \(M_0\) and \(M_0\) and found these to be reasonably well calibrated models. We now compare the predictive performance of the two, using both the CRPS and SCRPS metrics introduced previously.
The assess_prediction() function also allows the user to compute the CRPS and SCRPS through the argument metrics as shown below.
regularized_scores <- assess_prediction(list(M0 = M0_fit, M1 = M1_fit),
metrics = c("CRPS", "SCRPS"),
method = "regularized",
min_dist = 20, size = 9,
iter = 10, messages = FALSE)
cluster_scores <- assess_prediction(list(M0 = M0_fit, M1 = M1_fit),
metrics = c("CRPS", "SCRPS"),
method = "cluster", fold = 10,
messages = FALSE)Note that that assess_prediction() can be run once in order to obtain both the AnPIT diagnostic, introduced in Section 4.4.2; by default metrics = c("PIT", "CRPS", "SCRPS"), which computes all three diagnostic using the same test sets. The reason we have separated the AnPIT and the scores computations is only for pedagogical purposes.
By doing a summary() of the output in regularized_scores and cluster_scores, we obtain the average score across all test sets as defined in Equation 4.9.
summary(regularized_scores, view_all = FALSE)Summary of Cross-Validation Scores
----------------------------------
Model: M0
Overall average across test sets:
CRPS: 1.6102
SCRPS: -1.5991
PIT_area: 0.0894
Model: M1
Overall average across test sets:
CRPS: 1.6546
SCRPS: -1.6214
PIT_area: 0.0917summary(cluster_scores, view_all = FALSE)Summary of Cross-Validation Scores
----------------------------------
Model: M0
Overall average across test sets:
CRPS: 2.0627
SCRPS: -1.7287
PIT_area: 0.1683
Model: M1
Overall average across test sets:
CRPS: 1.9317
SCRPS: -1.7060
PIT_area: 0.1426In the code above, we have set view_all = FALSE to display only the overall average. If view_all = TRUE, the scores are instead reported for each test set, in addition to the overall average. The results from the summary() suggest that the two models have a similar predictive performance in terms of both CRPS and SCRPS, regardless of the sampling scheme used to create the test sets. To examine this further we can plot a point map of the CRPS and SCRPS using the locations of the test sets.
crps_maps <- plot(cluster_scores, metric = "CRPS")
scrps_maps <- plot(cluster_scores, metric = "SCRPS")
grid.arrange(crps_maps, scrps_maps, nrow = 2)
In the code above, we evaluated the cross-validation results generated using a clustering algorithm to split the original dataset. The map reveals that the scores across most locations are generally similar, with only three locations near the coast of Liberia showing slightly lower CRPS and SCRPS values for \(M_1\). This suggests a consistent predictive performance across the dataset, with some localized variations.
4.5 Simulation-based assessment of predictive performance
In Section 4.4, we illustrated how to evaluate predictive performance in relative terms, specifically by comparing candidate models under consideration. In this section, we demonstrate how simulation studies can be used to assess the predictive performance of a geostatistical model in predicting a spatially continuous target within a region of interest. Since a spatially continuous target cannot be directly observed, simulations provide a way to generate such targets under specified assumptions. This approach allows us to evaluate how effectively the model can recover the underlying spatial surface using suitable predictive performance metrics.
4.5.1 Step 1: Define the aim of the simulation study
As previously stated, one of the main reasons for carrying out a simulation study is the ability to generate spatial surfaces, which are typically unobservable in the real world. The purpose of the simulation is generally to evaluate how well a given geostatistical model can recover that surface and/or draw inferences on some of its properties.
To fully understand the how to carry out a simulation study, it is essential to define its two main components: the true model and the candidate models (which can indeed include multiple options).
The true model represents the underlying stochastic process that generates the data. In a simulation study, the true model is assumed to be known because it is explicitly specified by us. It serves as the benchmark against which other models are assessed.
The candidate models are the geostatistical models that are tested or evaluated during the simulation study. These models are proposed alternatives to approximate or recover the true spatial surface or its properties. Candidate models may differ in terms of their assumptions and model complexity.
It is important to highlight that the true model not only serves as a benchmark but can also be included among the candidate models. However, the true model does not always necessarily outperform other candidate models. In some cases, certain candidate models may provide better approximations of the spatial surface or its properties, particularly when the true model’s complexity is too high relative to the available sample size. Here, model complexity refers to the number of parameters in a model. For instance, a true geostatistical model with high complexity may include numerous covariates or a highly parameterised covariance function. When the sample size is insufficient to reliably estimate these parameters, simpler candidate models can sometimes yield better performances.
Consider an extreme example where the true model is a spatial Gaussian process \(S(x)\) with a scale parameter of \(\phi = 1\). If only a very small fraction of sampled locations are at distances below 1, estimating \(\phi\) accurately becomes highly challenging, especially with a small sample size. In particular the sparsity of close pairs reduces the information available to estimate the spatial dependence structure, making the estimation of \(\phi\) imprecise or even unreliable. In this case, a candidate model that does not attempt to explicitly estimate the spatial dependence (e.g., a model that relies solely on spatial covariates to interpolate the surface) might perform better in recovering the spatial structure. This example suggests that a simpler or alternative model can sometimes outperform the true model under specific sampling and data constraints.
At this point, a natural question arises: how should we choose the true model for a simulation study? The answer depends on the specific objectives of the study and what we aim to demonstrate. We now provide an example to demonstrate this point.
4.5.1.1 Example
In our analysis of the river blindness dataset, we now aim to evaluate the effectiveness of using elevation as a spatial predictor. To achieve this, we conduct a simulation study to compare the performance of two geostatistical models: one that incorporates elevation as a predictor and another that excludes it. This example will be revisited throughout the steps outlined in this section on simulation. We can chose the following true model (henceforth \(M_T\)) for the simulation \[M_{T}: \: \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_{0} + \beta_{1}e(x_i) + \beta_{2}\max\{e(x_i)-150, 0\} + S(x_i). \tag{4.10}\]
The true model, described by the equation above, assumes a nonlinear relationship with elevation based on a linear spline with a single knot at 150 meters. To simulate data from the true model, which constitutes the next step in our simulation study, we must first define the model parameter values. A reasonable approach is to use the maximum likelihood estimates obtained by fitting the binomial geostatistical model specified in Equation 4.10. These estimates can be derived as shown below.
set.seed(123)
true_model <- glgpm(npos ~ elevation + pmax(elevation - 150, 0) +
gp(),
denominator = ntest,
model_crs = 32629,
family = "binomial",
data = liberia)# The parameters of the true model
coef(true_model)$beta
(Intercept) elevation pmax(elevation - 150, 0)
-2.002985655 0.003160402 -0.001706721
$sigma2
[1] 0.3181267
$phi
[1] 41.02334A natural candidate model to evaluate the usefulness of elevation as a covariate is an intercept-only model (referred to here as \(M_C\)). This model is defined as: \[M_{C}: \: \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_{0} + S(x_i).\] where no covariates are included. This approach makes \(M_C\) a suitable baseline model, commonly selected when the objective is to assess the contribution of covariates to spatial prediction.
4.5.2 Step 2: Simulate the spatial surface and the data from the true model
In Section 3.3.3, we have shown how to simulate from a geostatistical model, when the goal is simulate data at the observed locations. When we want to simulate a spatial surface, a similar approach can be used. First of all, we need to generate a regular grid, say \(\tilde{X}\), covering the study region of interest, say \(A\). Based on the notation used so far, we use \(X\) to denote the set of observed locations which are used to simulate the outcome data \(Y\). Finally, we use \(X_{+} = (X, \tilde{X})\) to denote the full set of locations both on the grid and at the data points. We can then proceed through the following steps to simultaneously simulate values of the surface of interest over \(\tilde{X}\) and \(X\).
Compute the covariance matrix \(\Sigma_{+}\) using the combined set of locations \(X_{+}\).
If covariates are used, compute the covariates effects \(d(x)^\top \beta\) for all the locations in \(X_{+}\).
Simulate \(B\)-times from a multivariate Gaussian distribution with mean vector 0 and covariance \(\Sigma_{+}\) as previously defined. We denote the output form this step as \(S_{(j)}(\tilde{x})\), for simulations generate on \(\tilde{X}\), and with \(S_{(j)}(x)\) for simulations on \(X\), respectively, for \(j=1,\ldots,B\).
Compute the linear predictor at each locations of \(X_{+}\), by adding the covariates effects in 3 to the simulated random effects in 4, hence \(g\{\mu_{(j)}(\tilde{x})\} = d(\tilde{x})^\top + S_{(j)}(\tilde{x})\) for \(\tilde{X}\), and \(g\{\mu_{(j)}(x)\} = d(x)^\top + S_{(j)}(x)\) for \(X\).
Only for locations of the data \(X\), perform the steps 5A to 7A as explained in Section 3.3.3, to simulate the data \(Y\).
The function simulate_surface() in RiskMap can be used to carry out all these steps automatically as shown next.
4.5.2.1 Example
Following the steps previously outlined, we first create a grid within the boundaries of Liberia and extract the values of elevation on the grid from a raster.
shp <- gb_get_adm0("Liberia")
shp <- st_transform(shp, crs= 32629)
pred_grid <- create_grid(shp, spacing = 5)
pred_grid$elevation <- terra::extract(liberia_elev, st_coordinates(pred_grid))The function simulate_surface() also requires a sampling functions that generates the locations \(X\) from which data \(Y\) are simulated. In this example, we shall consider a sampling function that simply returns the locations of the original data.
sampling_f_lib <- function() {
coords_sf <- liberia[, "ntest"]
coords_sf <- st_transform(coords_sf, 32629)
coords_sf$units_m <- coords_sf$ntest
return(coords_sf)
}The output of the sampling function must be an sf object that must contain the number of tested people and its corresponding column must be named units_m. Also, the returned sf object must have a projection that correspond to the one used when fitting the model.
We are now ready to run the simulate_surface() function.
set.seed(123)
lib_surf_sim <- simulate_surface(n_sim = 200,
pred_grid = pred_grid,
formula = ~ elevation + pmax(elevation - 150, 0) +
gp(kappa = 0.5),
sampling_f = sampling_f_lib,
family = "binomial",
par0 = coef(true_model))In the above function, we have simulated \(B = 200\) simulations from the true model as specified in formula. The argument par0 is used to pass the values of the parameters of the true model.
We can now inspect the individual simulated surfaces using the plot_sim_surf() function, as demonstrated below, which produces Figure 4.24. This step is useful verifying that the simulations are consistent with expectations and free from any warning signs or anomalies.
plot_sim_surf(lib_surf_sim, sim = 1)
4.5.3 Step 3: Fit the candidate models to the simulated data and predict the target
To summarize, the results obtained from Step 1 and Step 2 are as follows:
- Simulated values for the linear predictor of the true model (\(M_T\)) over a grid covering the study area.
- Simulated values for the linear predictor of \(M_T\) at the data locations.
- Simulated realizations of the data at the data locations.
Next, we use the third output listed above to fit the candidate models (\(M_C\)) and make predictions over the grid. Note that each of the outputs above has been generated \(B\) times, which denotes the number of simulations.
Using the notation introduced in Section 4.2, let \(T(x)\) denote the spatially continuous target. When the predictive target is at the areal level, the necessary adjustments to Step 3 and Step 4 will be detailed in Section 4.5.5.
The outputs generated in these steps include predictive samples of the target over the grid, represented by \(T_{(j)}^{(h)}(\tilde{x})\), where: \(j\) is the index of the \(j\)-th simulated dataset, for \(j = 1, \ldots, B\); \(h\) is the index of the Monte Carlo samples, for \(h = 1, \ldots, N\), generated by the Monte Carlo Markov Chain (MCMC) algorithm for a given dataset \(j\).
An example in R for Step 3 will be provided in the next section, where we shall use the function assess_simulation(), from the RiskMap package, which combines Step 3 and Step 4 into a single step.
4.5.4 Step 4: Summarize the results using a pre-defined objective function
In this final step, we compare the predictions of the target \(T(x)\) against its true values, originally simulated in Step 1. To do this, we need a metric to summarize the results across all \(B\) simulations for each candidate model \(M_C\).
One approach is to use the predictive mean from \(M_C\) in each simulation and evaluate how it differs from \(T(x)\). A natural choice for this comparison is the mean squared error (MSE), defined as:
\[\frac{1}{B} \sum_{j=1}^B \sum_{i=1}^q \left(T(\tilde{x}_i) - \hat{T}_{C}(\tilde{x}_i)\right)^2. \tag{4.11}\]
where \(\tilde{x}_i\) is the \(i\)-th location on the grid \(\tilde{X}\), \(T(\tilde{x}_i)\) represents the true value of the target, and \(\hat{T}_{C}(\tilde{x}_i)\) is the predictive mean from \(M_C\) at location \(\tilde{x}_i\).
Other options for summarizing the results include using the bias or the absolute median deviation. The bias measures the systematic error between the predicted and true values and is defined as:
\[\frac{1}{B} \sum_{j=1}^B \sum_{i=1}^q \left(T(\tilde{x}_i) - \hat{T}_{C}(\tilde{x}_i)\right)\]
The median absolute deviation (MAD) focuses on the central tendency of the errors, providing a robust measure of variability. It is defined as:
\[\text{median}\left(\left|T(\tilde{x}_i) - T^*_{C}(\tilde{x}_i)\right|\right),\] where \(T^*_{C}(\tilde{x}_i)\) is the median across all \(N\) the Monte Carlo samples \(T_{(j)}^{(h)}(\tilde{x})\).
4.5.4.1 Example
Using the output from the simulate_surface() function, named lib_surf_sim, we can use as the input to the function assess_simulation() which performs Step 3 and Step 4, previously outlined.
res_sim_grid <- assess_simulation(lib_surf_sim,
models = list(M_T = ~ elevation +
pmax(elevation - 150, 0) +
gp(),
M_C = ~ gp()),
f_grid_target = function(x) exp(x)/(1+exp(x)),
pred_objective = "mse",
spatial_scale = "grid")In this function, we specify the models to be assessed using the models argument, which is a list of formula objects. The formula M_T represents the true model \(M_T\), while M_C corresponds to the intercept-only model. Note that the coordinate names used in the gp() function within the formulas are converted according to the CRS specified in lib_suf_sim.
The argument f_grid_target defines our target function \(T(x)\) over the grid as a transformation of the linear predictor. Specifically, exp(x) / (1 + exp(x)) expresses disease prevalence as our predictive target over the grid. The pred_objective argument specifies the objective function from Step 4, with pred_objective = "mse" setting it to the mean squared error from Equation 4.11. Finally, spatial_scle = "grid" indicates that the objective function applies to a pixel-level target according to the grid stored in lib_surf_sim.
We can now look at the summary of res_sim_grid.
summary(res_sim_grid)Summary of Simulation ResultsThe results reported here present the average and standard deviation across 200 simulations for the models \(M_C\) and \(M_T\). The differences between the two models are negligible, indicating that the point predictions generated by \(M_C\) and \(M_T\) perform similarly in recovering the true prevalence surface for river blindness in Liberia.
This result should not be surprising, as estimating regression relationships with prevalence data can be challenging, potentially limiting their predictive advantage. However, incorporating covariates enhances model interpretability, allowing us to understand why the model makes certain predictions in specific areas, which an intercept-only model cannot provide. Thus, even if the model with covariates does not achieve the best predictive performance, we may still choose to retain the covariates, in this case elevation, for their interpretative value.
4.5.5 Assessessment of threshold-based classification of spatial units using geostatistical models
In this section, we address the specific problem of evaluating how effectively a geostatistical model classifies areal units into predefined categories of a disease indicator, such as prevalence. This issue is particularly relevant when dealing with areal-level targets (Section 4.3) as opposed to spatially continuous targets (Section 4.2). This prediction problem arises when policy decisions regarding interventions are often based on the classification of the areal-level target into predefined categories at the level of administrative regions; see, for example, Johnson et al. (2021) and Puranik et al. (2024).
To formalize this problem mathematically, let \(\mathcal{T}_i\) represent the true value of the areal-level target (e.g., the average prevalence over a district), and let \(\hat{\mathcal{T}}_{C,i}\) denote the corresponding prediction produced by the candidate model \(M_C\). Note that \(\mathcal{T}_i\) can be easily calculated in Step 1 during simulations from the true model (Section 4.5.1).
Now, consider that the target \(\mathcal{T}_i\) is categorized into predefined classes, which are typically thresholds or ranges of the disease indicator (such as values that define low, medium, and high levels). Denote the true membership class of areal unit \(i\) by \(\mathcal{D}_i\), which is determined based on \(\mathcal{T}_i\), and the predicted membership class by \(\hat{\mathcal{D}}_{C,i}\), determined based on \(\hat{\mathcal{T}}_{C,i}\). The goal is to assess the overall quality of the classification, which involves evaluating how well the predicted classes \(\hat{\mathcal{D}}_{C,i}\) agree with the true classes \(\mathcal{C}_i\) across all areal units and simulations. This assessment provides insight into the model’s ability to accurately classify regions according to the chosen disease indicator and thus guide effective policy decisions.
The required changes to Step 3 (Section 4.5.3) consist of simply carrying the predictions over the areal units as explained in Section 4.3.
For Step 4, which concerns how we summarize the agreement between the classification generated by \(M_C\) and the true classes \(\mathcal{D}_i\), we can proceed as follows. First, we distinguish between indicators that can be computed for a given region and those that summarize performance across all regions.
Assessing classification performance can be approached similarly to evaluating a diagnostic test. Here, true positives, true negatives, false positives, and false negatives are defined at the level of individual spatial units that subdivide the study area. Hence, key indicators include:
Sensitivity: The proportion of areal units correctly classified as belonging to the true membership class when conditioned on that true class:
\[\text{Sensitivity} = \frac{\text{True Positives (TP)}}{\text{True Positives (TP)} + \text{False Negatives (FN)}}\]
Specificity (True Negative Rate): The proportion of areal units not belonging to a given membership class that are correctly classified as not being in that class:
\[\text{Specificity} = \frac{\text{True Negatives (TN)}}{\text{True Negatives (TN)} + \text{False Positives (FP)}}\]
Positive Predictive Value (PPV): The proportion of areal units predicted to belong to a membership class that are correctly classified, conditioned on the predicted membership class:
\[\text{PPV} = \frac{\text{True Positives (TP)}}{\text{True Positives (TP)} + \text{False Positives (FP)}}\]
Negative Predictive Value (NPV): The proportion of areal units predicted not to belong to a membership class that are correctly classified, conditioned on the predicted membership class:
\[\text{NPV} = \frac{\text{True Negatives (TN)}}{\text{True Negatives (TN)} + \text{False Negatives (FN)}}\]
Correct Classification Rate: The overall proportion of areal units whose predicted membership class matches the true membership class:
\[\text{Correct Classification Rate} = \frac{\text{True Positives (TP)} + \text{True Negatives (TN)}}{\text{Total Cases (TP + TN + FP + FN)}}\]
Note that for multiple classes, the indicators above are computed by comparing each class in turn against all other classes combined. This allows us to create a 2 by 2 table, enabling the calculation of the indicators above.
The choice between sensitivity, specificity, PPV, and NPV depends on the specific public health context and the priorities of the decision-making process. Sensitivity and specificity are conditioned on the true membership class, making them useful for assessing the model’s performance in terms of correctly identifying or excluding true positives or negatives. PPV and NPV are conditioned on the predicted membership class, making them more informative about the trustworthiness of predictions, especially when prevalence varies significantly across regions.
We can thus summarize the difference of each indicator as follows.
Sensitivity is most appropriate when the priority is to ensure that true positive cases (e.g., high-prevalence regions) are not missed. Sensitivity is useful in contexts where failing to identify true cases could have severe consequences, such as delaying critical interventions or underestimating resource needs. For example, during an outbreak, ensuring all high-risk areas are flagged is crucial, even if this increases the number of false positives.
Specificity is more suitable when the focus is on minimizing false positives, which is critical in contexts where over-classification could lead to unnecessary interventions or misallocation of resources. Specificity may be particularly important in low-prevalence scenarios, where the cost of acting on false alarms could outweigh the benefits of capturing all true positives.
PPV is appropriate when the emphasis is on the reliability of positive predictions, i.e., how confident we can be that regions classified as high-prevalence are truly high-prevalence. PPV is heavily influenced by the prevalence of the condition being classified, making it especially relevant in scenarios where actionable decisions depend on the certainty of a positive classification.
NPV is most relevant when it is critical to trust negative predictions, i.e., ensuring that regions predicted as low-prevalence are indeed low-prevalence. NPV is particularly important in low-prevalence settings, where a high proportion of predicted negatives are true negatives.
Ultimately, the selection of indicators should result from a dialogue between statisticians, epidemiologists, and policy makers, ensuring alignment with the specific priorities and policy decisions being informed by the model’s classifications.
4.5.5.1 Example
Using the same models, \(M_T\) and \(M_C\), from the river blindness example in Liberia, we now compare their performance in classifying admin level 1 regions into two prevalence categories, \((0, 0.2)\) and $(0.2, 1) $. The assessment is conducted using assess_simulation(), where the prevalence surface is first aggregated at the areal level before classification.
# Admin level 1 boundaries
shp_adm <- gb_get_adm1("Liberia")
shp_adm <- st_transform(shp_adm, crs= 32629)
res_sim_area <- assess_simulation(lib_surf_sim,
models = list(M_T = ~ elevation +
pmax(elevation - 150, 0) +
gp(),
M_C = ~ gp()),
f_grid_target = function(x) exp(x)/(1+exp(x)),
f_area_target = mean,
pred_objective = "classify", shp = shp_adm,
categories = c(0,0.2,1),
spatial_scale = "area")The implementation first loads the admin level 1 shapefiles for Liberia, transforming them to the appropriate coordinate reference system. The arugment f_area_target is set to mean, ensuring that the target prevalence for each region is computed as the average over the area, corresponding to
\[\frac{1}{|A|} \int_A T(x) \, dx\] The classification is performed by specifying pred_objective = "classify", meaning the models are evaluated based on their ability to assign regions to the correct prevalence class. The categories argument defines the classification thresholds, requiring at least three values to distinguish between two classes. The spatial_scale is set to "area", meaning predictions are made at the level of administrative regions which are specified by the shapefile shp_adm, passed to the argument shp.
Finally, we can examine the results of the simulation through the summary function.
summary(res_sim_area)Summary of Simulation ResultsIn this example, the performance of models \(M_T\) and \(M_C\) is comparable across all considered indicators. This similarity arises for the same reasons previously discussed examining the quality of prevalence surface predictions.
4.6 Theory
4.6.1 Expression of the predictive distribution for a linear geostatistical model
In this section, we consider the simple case of a linear geostatistical model with a single observation per location.
Let \(Y = (Y_1, \ldots, Y_n)^\top\) be the vector of data associated with spatial locations \(x_1, \ldots, x_n\). The model assumes that:
\[\begin{equation} Y_i = d(x_i)^\top {\beta} + S(x_i) + Z_i, \quad i = 1, \ldots, n, \end{equation}\]
where \(d(x_i)\) is a vector of covariates at location \(x_i\), \({\beta}\) is a vector of regression coefficients, \(S(x)\) is a stationary Gaussian process with mean zero and Matern correlation function \(\rho(u)\) (where \(u\) is the distance between locations), and \(Z_i \sim \mathcal{N}(0, \tau^2)\) represents independent measurement errors (nugget effect).
We first consider the simple case of a linear target, hence \(T(x) = d(x)^\top {\beta} + S(x)\). The goal is to predict the unobserved process \(T^* = (T_{n+1}, \ldots, T_{n+q})^\top\) at new locations \(x_{n+1}, \ldots, x_{n+q}\). The predictive distribution of \(T^*\), conditional on the observed data \(Y\), follows a multivariate Gaussian distribution:
\[\begin{equation} T^* | Y \sim \mathcal{N}(\mu^*, \Sigma^*), \end{equation}\]
where the conditional mean and covariance are given by:
\[\begin{align} \mu^* &= D^* {\beta} + C \Sigma^{-1} (Y - D {\beta}), \\ {\Sigma}^* &= V - C \Sigma^{-1} C^\top. \end{align}\]
Here, \(D^*\) is the design matrix for the prediction locations, \(C\) is the cross-covariance matrix between \(T^*\) and \(Y\), \(\Sigma\) is the covariance matrix of \(Y\), and \(V\) is the covariance matrix of \(T^*\). Using a multivariate Gaussian distribution with mean \(\mu^*\) and covariance matrix \(\Sigma^*\), we are now able to predict \(T^*\) and obtain summaries of uncertainties.
For non-linear prediction targets, i.e. targets that involve non-linear transformations of the spatial process \(S(x)\), say \(W(T^*)\), Monte Carlo methods can be employed. More specifically, we generate samples \(T_{(j)}^*\), for \(j=1,\ldots,B\) from the predictive distribution \(T^* | Y\) using the multivariate Gaussian distribution with mean \(\boldsymbol{\mu}^*\) and covariance \(\boldsymbol{\Sigma}^*\). For each sample \(T^{*}_{(j)}\), the non-linear transformation \(W(T^{*}_{(j)} = W^{*}_{(j)}\) of interest is then computed. The predictive summary statistics, such as the mean, variance, or quantiles of \(W(T^*)\), are then estimated by averaging over the simulated samples, e.g.
\[\begin{align} E[W(T^*) | Y] &\approx \frac{1}{B} \sum_{J=1}^B W^{*}_{(j)}), \\ \text{Var}[W(T^*) | Y] &\approx \frac{1}{B-1} \sum_{j=1}^B \left( W^{*}_{(j)} - E[W(T^*)] \right)^2. \end{align}\]
For more detailed explanations on this we refer you to Chapter 6 of Diggle and Ribeiro (2007).
4.6.2 A brief overview of scoring rules
We now provide an overview of the theoretical framework of scoring rules, of which the continuous ranked probability score (CRPS) and its scaled variant (SCRPS), previously presented, are an example. Scoring rules are a tool used in the evaluation of probabilistic forecasts. As shown for CRPS and SCRPS, scoring rules are designed to measure how well a predictive distribution aligns with observed outcomes, balancing both calibration and sharpness.
In general, a scoring rule is a function \(G(F, y)\) that assigns a numerical score to a predictive distribution \(F\) and an observed outcome \(y\). Lower scores typically indicate better predictive performance, as they reflect a closer match between the predicted distribution and the observed value.
A central concept in the theory of scoring rules is that of proper scoring rules. A scoring rule is said to be proper if it favours honest predictions or, in other words, if the expected score is minimized when the predictive distribution \(F\) matches the true data-generating distribution \(F^*\). More formally, a scoring rule \(G\) is proper if, for all distributions \(F\) and \(F^*\), we have: \[E_{Y \sim F^*}[G(F^*, Y)] \leq E_{Y \sim F^*}[G(F, Y)],\] with equality if and only if \(F = F^*\) when the rule is strictly proper. This property is crucial because it ensures that forecasters cannot improve their scores by manipulating their predictions to deviate from their true beliefs.
Proper scoring rules encourage the selection of models that yield predictions that are both well-calibrated (reflecting the true distribution of outcomes) and sharp (concentrated around likely values). However, this does not imply that predictions from the (unknown) true model \(F^*\) will always be favoured in practice, when dealing with finite samples. In other words, even if two scoring rules are proper, this does not guarantee that the true data-generating model will always achieve a lower score compared to a misspecified model in finite samples. Sample variability can in fact lead to the selection of a misspecified model by chance, especially when the sample size is small or the data are noisy. Model complexity is another important factor. If the true model \(F^*\) is highly parameterised and the sample size is insufficient to estimate its parameters accurately, it may perform poorly in predictive tasks due to overfitting or parameter estimation error. In contrast, a simpler misspecified model, despite not capturing the full complexity of the data-generating process, might generalize better in such scenarios. This is related to the bias-variance trade-off, where a model with higher bias (due to misspecification) can still achieve lower overall predictive error if it has substantially lower variance.
The CRPS and SCRPS are examples of a strictly proper scoring rule for continuous outcomes. Other scoring rules are the following:
Logarithmic Score (or Log Score): For a predictive density \(f(y)\), the log score is defined as: \[G_{\text{log}}(f, y) = -\log f(y).\] This rule heavily penalizes low-probability predictions for observed outcomes, making it sensitive to the tails of the predictive distribution.
Brier Score: For binary outcomes, the Brier score measures the mean squared difference between the predicted probability \(p\) of an event and the actual outcome \(y \in \{0, 1\}\): \[G_{\text{Brier}}(p, y) = (p - y)^2.\] The Brier score is both proper and easy to interpret, with connections to mean squared error in classical regression.
Energy Score: A generalization of the CRPS to multivariate outcomes, the energy score for a predictive distribution \(F\) and observation \(y\) is given by: \[G_{\text{energy}}(F, y) = E_{X, X' \sim F}[\|X - X'\|] - 2 E_{X \sim F}[\|X - y\|],\] where \(\|\cdot\|\) denotes a norm, typically the Euclidean norm. The energy score is particularly useful for evaluating spatial or multivariate forecasts.
Many proper scoring rules, including the CRPS, can be decomposed into components that separately capture calibration and sharpness. This decomposition helps to understand the different aspects of predictive performance. Specifically, the expected score can often be expressed as: \[E[G(F, Y)] = \underbrace{\text{Uncertainty}}_{\text{inherent in data}} - \underbrace{\text{Resolution}}_{\text{sharpness of predictions}} + \underbrace{\text{Reliability}}_{\text{calibration error}},\] where: - Uncertainty represents the inherent variability in the data, independent of the model. - Resolution quantifies the model’s ability to produce sharp, confident forecasts. - Reliability (or calibration error) measures the degree to which the forecast probabilities align with observed frequencies.
This decomposition highlights an important trade-off: improving sharpness without compromising calibration is the hallmark of a good probabilistic model. A model that produces overly sharp forecasts may suffer from poor calibration, while a model that is perfectly calibrated but overly diffuse lacks resolution.
As we have shown in this chapter, scoring rules are particularly appropriate in cross-validation frameworks, where the goal is to evaluate how well a model predicts future data realizations. When instead conducting simulation studies, we generate data from a known simulated true surface which becomes our benchmark. Hence, in simulation studies, we do not perform cross-validation but, for a given simulation, we compare the model predictions against the true surface. In such cases, scoring rules are not used because they are designed to assess predictive performance based on observed data, where both calibration and sharpness are critical considerations. Hence, in simulation studies, prediction performance metrics, such as bias, mean squared error (MSE), or other distance-based metrics that directly quantify the deviation from the known truth can be used to quantify the predictive performance of a model.
For a more detailed explanation of scoring rules, Gneiting et al. (2007) is a recommended read.
4.7 Summary
In this chapter, we have demonstrated how geostatistical models can be used to draw predictive inference for targets that are spatially continuous (e.g., a prevalence surface within a study region) and aggregated at the areal level (e.g., the average disease prevalence in a district). The fundamental concept underpinning geostatistical prediction is that of the predictive distribution, which we define as the distribution of the target \(T(x)\) conditioned on the data \(y\). This distribution is the key element that enables the derivation and computation of point predictions, predictive intervals and the probability that the target falls within a pre-specified interval or above/below a given threshold.
To assess the predictive performance of a geostatistical model, we have considered two complementary approaches: cross-validation and simulation studies.
In the cross-validation framework, the concepts of calibration and sharpness are introduced to evaluate the quality of spatial predictions. Calibration refers to the agreement between the predicted probability of an event occurring and its actual observed frequency. Sharpness, on the other hand, quantifies the concentration of the predictive distribution around its point estimate. However, sharper predictions are not necessarily better; they must be well-calibrated to be meaningful. Thus, it is appropriate to discuss sharpness only in the context of well calibrated models.
A simulation study, unlike cross-validation, allows us to assess the ability of our model to predict an entire spatial surface within a study area or a predictive target defined over spatial units (e.g., districts) that make up the study area. This is possible because, in a simulation study, we can generate data from a true model and then evaluate the quality of inferences from candidate models on the simulated surface. The way we summarize predictive performance in a simulation study depends on the objective. If the goal is to assess how closely our point predictions align with the true surface, we can use metrics such as the mean squared error (Equation 4.11), averaged over the grid used to approximate the surface. Conversely, if the objective is to evaluate classification performance by categorizing spatial units into predefined health outcome categories, we can use sensitivity, specificity, negative predictive value, and positive predictive value. The choice of metric depends on whether the priority is minimizing errors in detecting true positives or true negatives (sensitivity and specificity), or ensuring that the spatial units classified as positive or negative are truly so (positive predictive value and negative predictive value). A summary of the metrics considered in this chapter is given in Table 4.3.
| Diagnostic/Metric | Performance characteristic | Recommended use in |
|---|---|---|
| Bias | Accuracy | Simulation |
| Continuous ranked probability score (CRPS) | Calibration and sharpness | Cross-validation |
| Coverage probability | Calibration | Cross-validation and simulation |
| Mean square error | Accuracy | Simulation |
| Nonrandomized probability integral transform (PIT) | Calibration | Cross-validation |
| PIT | Calibration | Cross-validation |
| Scaled CRPS | Calibration and sharpness | Cross-validation |
4.8 FAQs
- What should be the spatial resolution of the grid used for prediction?
The choice of spatial resolution for the prediction grid involves a trade-off between computational efficiency and the ability to capture spatial correlation. Higher spatial resolution provides more detailed predictions but increases computational demands. A reasonable approach is to ensure that the distance between adjacent grid pixels is small enough to maintain a high spatial correlation. For instance, selecting a resolution where the correlation between neighboring pixels is at least 0.95 can be a useful guideline. However, this threshold is somewhat arbitrary, similar to the convention of using 0.95 in confidence intervals. Ultimately, the resolution should be chosen based on the spatial scale of variation in the data, while avoiding excessive computational burden.
- I have covariates with different spatial resolutions. Which spatial resolution should I use?
Two main approaches can be considered for handling covariates at different spatial resolutions in geostatistical analysis. The first involves aggregating all covariates to the coarsest resolution, ensuring consistency and avoiding artificial detail. This prevents high-resolution predictions based on low-resolution data but may sacrifice fine-scale spatial variation. The second approach disaggregates covariates to the finest resolution, enabling higher-detail predictions but risking false precision and artifacts if the data do not support such resolution. To assess the impact on inferences, one might apply both approaches and evaluate how they influence the final results.
- My prediction map looks like noise. What could be the error?
The prediction map appearing noisy could be due to several factors. First, one should verify that the spatial grid was created correctly by plotting it alongside the original data to ensure alignment and consistency. Second, the spatial correlation in the data might be very small relative to the resolution of the grid, causing the model to struggle with capturing meaningful patterns. To address this, increasing the spatial resolution of the grid may help, as finer resolution can better capture the underlying spatial structure.
- The main goal of my analysis is to carry out predictions of the main outcome. However, can I also use the model to draw inferences on the association between the outcome and the covariates that I have included in the model?
The model developed for prediction, in theory, should not be primarily used to draw inferences on the regression relationships between the outcome and covariates. Understanding regression relationships and making predictions are distinct goals that often require different modeling approaches. While it is possible to interpret the estimated regression coefficients, caution is necessary, as these associations are not causative. Predictive models are optimized for accuracy in predicting th outcome rather than for identifying or explaining underlying relationships, so any inferences drawn should be treated as exploratory rather than definitive. An recommended reading on this topic is Shmueli (2010).
- In my model I have some covariates that are attached to the individual (e.g. age, gender, occupation, etc.). How do I specify them when doing prediction?
The answer to this question is given in Section 4.2.4.
4.9 Review questions
- What are the steps that needs to taken when performing spatial prediction?
- Define the predictive distribution of a predictive target and explain how this is used to answer a research question.
- Distinguish between spatially continuous and areal-level targets, giving at least two examples for each.
- Explain the difference between marginal and joint predictions, and for which targets each can be used.
- Define predictive performance in terms of calibration and sharpness.
- Explain how calibration can be assessed using graphical tools, by distinguishing between continuous and count outcomes.
- Give examples of metrics that can be used to evaluate calibration and sharpness.
- Explain how simulation studies can be carried to assess the predictive performance of geostatistical models and how this differ from cross-validation.
4.10 Exercises
Using the code from Section 4.4, redraw the plot of Figure 4.16 by choosing different values for the dispersion parameter of the Negative Binomial distribution. What happens when \(\alpha\) is close to 0 and why?
Consider the example in Section 4.4.2.1. Fit a third model, that extends \(M1\) by adding a nugget effect, hence \[M_{2} \: : \: \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_0 + \beta_{1}e(x_i) + \beta_{2}\max\{e(x_i) - 150, 0\} + S(x_i) + Z_i.\] Generate the AnPIT plot for \(M_2\) and compare it to those of \(M_0\) and \(M_1\). Make sure to apply the function
assess_prediction()to all three models so that the same test sets are used. Is \(M_2\) showing a better calibration? Why?In the analysis of the Anopheles gambiae mosquitoes data, consider two models. One that uses elevation, with linear predictor \[M_1\: : \: \log\{\lambda(x_i)\} = \beta_0 + \beta_1 e(x_i) + S(x_i)\] where \(e(x)\) is the elevation in meters at location \(x\) and \(S(x_i)\) is a stationary and isotropic Gaussian process with variance \(\sigma^2\) and exponential correlation function \(\rho(u) = \exp\{-u/\phi\}, \phi>0\). Consider also an intercept only model \[M_0\: : \: \log\{\lambda(x_i)\} = \beta_0 + S(x_i).\] After fitting both models to the data, assess if they are well calibrated models using the average nonrandomized probability integral transform (AnPIT). Assuming both models are well-calibrated, compare the predictive performance using continuous ranked probability score (CRPS) and the scaled CRPS. In doing this, use different sampling schemes to split the data, as illustrated in Section 4.4.1, and assess if and how the adopted approach to split the data affects the results.
Consider the Liberia dataset on river blindness. Carry out a simulation study where the objective is to assess the impact of mispecifying non-linear relationships between elevation and prevalence. Consider Equation 4.10 as the true model and a model that specifies the relationship as linear on the logit scale as the candidate model, hence \[M_{C}: \: \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_{0} + \beta_{1}e(x_i) + S(x_i).\] Carry out 200 iterations in the simulation and define a suitable metric of predictive performance to compare the two models. Provide a clear answer on the impact of assuming a linear relationship with elevation when this is not valid.