# Compute prevalence
liberia$prev <- liberia$npos / liberia$ntest
# Scatterplot of prevalence vs. elevation
ggplot(liberia, aes(x = elevation, y = prev)) +
geom_point() +
labs(x = "Elevation (metres)", y = "Prevalence")
3 Model formulation and parameter estimation
List of the main functions used in the chapter
| Function | R Package | Used for |
|---|---|---|
lmer() |
lme4 | Fitting linear mixed models |
glmer() |
lme4 | Fitting generalized linear mixed models |
glgpm() |
RiskMap | Fitting generalized linear geostatistical models |
variogram() |
RiskMap | Computing the empirical variogram and a global envelope for spatial independence |
plot_variogram() |
RiskMap | Plotting an empirical variogram and its global envelope |
summarise_distance() |
RiskMap | Summarising distances between sampling locations |
to_table() |
RiskMap | Formatting model estimates as a table |
set_control_mcmc() |
RiskMap | Configuring Markov chain Monte Carlo settings |
set_control_mcml() |
RiskMap | Configuring repeated Monte Carlo maximum-likelihood updates |
plot_mcmc() |
RiskMap | Checking Markov chain Monte Carlo output |
simulate_glgpm() |
RiskMap | Simulating replicated data from a fitted geostatistical model |
simulated_data() |
RiskMap | Extracting one simulated dataset |
This chapter introduces the process of moving from exploratory analysis to the formal specification of statistical models, through to parameter estimation and interpretation of results. We first outline the main steps involved, highlight the role of covariates, and present simple approaches for exploring and modelling non-linear relationships. We then turn to the assessment of residual spatial correlation, which motivates the use of geostatistical models, and describe how maximum likelihood methods can be applied for their estimation.
3.1 Exploratory analysis
As illustrated in Figure 1.8, exploratory analysis is the first step in any statistical analysis. This stage is crucial for guiding how covariates should be incorporated into the model and for assessing whether the remaining unexplained variation exhibits spatial correlation. In the case of count data, this stage also involves examining the presence of overdispersion, a phenomenon that often arises when important covariates are omitted or when there is unaccounted heterogeneity in the data. In the following sections, we explain what overdispersion is in more detail and why checking for its presence is an essential prerequisite before fitting geostatistical models to count data.
3.1.1 Exploring associations with risk factors using count data
Assessment of the association between the health outcome of interest and non-categorical (i.e. continuous) risk factors can be carried out using graphical tools, such as scatter plots. The graphical inspection of the empirical association between the outcome and the covariates is especially useful to identify non-linear patterns in the relationship which should then be accounted for in the model formulation.
In this section, we focus on the case where the observed outcome is a count, which requires a different treatment from continuously measured outcomes. The latter have been extensively discussed in the statistical literature and are covered in most standard textbooks (see, for example, Chapter 1 of Weisberg (2014)).
3.1.1.1 When the outcome is an aggregated count
Let us first consider the example of the river blindness data in Liberia (Section 1.4.2), and examine the association between prevalence and elevation. We first generate a plot of the prevalence against the measured elevation at each of the sample locations
The plot shown in Figure 3.1 shows that, as elevation increases from 0 to around 150 metres, prevalence rapidly increases to around 0.25 and, for elevation values above 150 metres, the relationship levels off. How can we account for this in a regression model? To answer this question rigorously, the plot in Figure 3.1 cannot be used. This is because, when the modelled outcome is a bounded Binomial count, regression relationships are usually specified on the logit-transformed prevalence (log-odds) scale; see Table 1.3 in Section Section 1.5. To explore regression relationships in the case of prevalence data, it is convenient to use the so-called empirical logit in place of the empirical prevalence. The empirical logit is defined as \[ l_{i} = \log\left\{\frac{y_i + 1/2}{n_i - y_i + 1/2}\right\} \tag{3.1}\]
where \(y_i\) denotes the number of individuals who tested positive for river blindness and \(n_i\) is the total number of people tested at location \(i\). Rather than applying the standard logit transformation to the empirical prevalence, we use the empirical logit transformation, implemented via the elogit() function in the RiskMap package. This transformation introduces a small adjustment to both the numerator and denominator, ensuring that the transformed values remain finite even when the observed prevalence is exactly 0 or 1, for which the standard logit is undefined.
# Empirical logit transformation
liberia$elogit <- elogit(y = liberia$npos, m = liberia$ntest)
# Scatterplot of elogit vs. elevation with different potential models
ggplot(liberia, aes(x = elevation, y = elogit)) +
geom_point() +
# Adding a smoothing spline
geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE) +
# Adding linear regression fit with log-transformed elevation
geom_smooth(method = "lm", formula = y ~ log(x),
col = "darkgreen", lty = "dashed", se = FALSE) +
# Adding linear regression fit with change point at 150 metres
geom_smooth(method = "lm", formula = y ~ x + pmax(x - 150, 0),
col = "orange", lty = "dashed", se = FALSE) +
labs(x = "Elevation (metres)", y = "Empirical logit")
Figure 3.2 shows the scatter plot of the empirical logit against elevation. In this plot, we have also added three curves through the geom_smooth() function from the ggplot2 package. Using this function, we first specify method = "gam" to add a penalised smoothing spline (Wood 2017), represented by the blue solid line. The smoothing spline allows us to better discern the type of relationship and how to best capture it using a standard regression approach. As we can see from Figure 3.2, the smoothing spline corroborates our initial observation of a positive relationship up to about 150 metres, followed by a plateau.
To capture this non-linear relationship, we can use the two following approaches. The first is based on a simple log-transformation of elevation and is represented in Figure 3.2 by the green line. If we express this relationship using a standard Binomial regression model, this would take the form \[ \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_0 + \beta_1 \log\{d(x_i)\} \tag{3.2}\] where \(p(x_i)\) denotes the probability of river blindness infection at sampled location \(x_i\), and \(d(x_i)\) is the elevation at that location.
Alternatively, the non-linear effect of elevation on prevalence could be captured using a linear spline. Put in simple terms, we want to fit a linear regression model that allows for a change in slope above 150 metres. Formally, this is expressed in a Binomial regression model as \[
\log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_0 + \beta_1 d(x_i) + \beta_{2} \max\{d(x_i)-150, 0\}.
\tag{3.3}\] Based on the equation above, the effect of elevation below 150 metres is quantified by the parameter \(\beta_1\). Above 150 metres, instead, the effect of elevation becomes \(\beta_1 + \beta_2\). Note that the function pmax() (and not the standard base function max()) should be used in R when the computation of the maximum between a scalar value and each of the components of a numeric vector is required.
Before proceeding further, it is important to explain the differences between the use of the logarithmic transformation (Equation 3.2) and the linear spline (Equation 3.3). We observe that both curves provide a similar fit to the data, with larger differences observed for larger values in elevation, where the log-transformed elevation models yield larger values for the predicted prevalence. This also suggests that if we were to extrapolate the predictions beyond 600 metres in elevation the implied pattern by the model with the log-transformed elevation would predict an increasingly larger elevation, which is unrealistic, since the fly that transmits the disease cannot breed at those altitudes. The linear spline model instead would generate predictions that would be very similar to those observed between 150 and 600 metres. From this point of view, the linear spline model would thus have more scientific validity than the other model. However, which of the two approaches should be chosen to model the effect of elevation is a question that closely depends on the research question to be addressed.
If the interest of the study was in better understanding the association between elevation and prevalence, the linear spline model does not only provide a more credible explanation but also its regression parameters can be more easily interpreted. In fact, for a unit increase in elevation, the multiplicative change in the odds for river blindness is \(\exp\{\beta_1\}\), if elevation is below 150 metres, and \(\exp\{\beta_1+\beta_2\}\), if elevation is above 150 metres. When instead we use the log-transformed elevation, the interpretation of \(\beta_1\) in Equation 3.2 is more complicated, as it is based on the multiplicative increase in elevation by the same amount given by the base of the algorithm, which is about \(e \approx 2.718\). To avoid this, one could rescale the regression coefficient as, for example, \(\beta_1/\log_{2}(e)\) which would be interpreted as the multiplicative change in the odds for river blindness for a doubling in elevation. However, a doubling in elevation is less meaningful when considering larger values of elevation.
Note
The letter \(e\) stands for the so-called Euler’s number and represents the base of the natural logarithm. In the book, we write \(\log(\cdot)\) to mean the “natural logarithm of \(\cdot\)”.
When the goal of statistical analysis is instead in developing a predictive model for the outcome of interest, the explanatory power and interpretability of the model may be of less concern. For this reason, the model with the log-transformed elevation could be preferred over the model with the linear spline, if it is shown to yield more predictive power. We will come back to this point again in Chapter 4, where we will show how to assess and compare the predictive performance of different geostatistical models.
The other type of aggregated count data that we consider are unbounded counts. The Anopheles mosquitoes dataset (Section 1.4.4) is an example of this, since there is no upper limit to the number of mosquitoes that can be trapped at a location. Let us consider the covariate represented by elevation. In this case, the simplest model that can be used to analyse the data is a Poisson regression, where the linear predictor is defined as the log of the mean number of mosquitoes (Table 1.3). Hence, exploratory plots for the association with covariates should be generated using the log transformed counts of mosquitoes. In this instance, to avoid taking the log of zero, we can add 1 to the reported counts, if required. The variable An.gambiae in the anopheles dataset does not contain any 0, hence we simply apply the log transformation without adding 1.
# Apply log-transformation to Anopheles gambiae counts
anopheles$log_counts <- log(anopheles$An.gambiae)
# Scatterplot of log counts vs. elevation
ggplot(anopheles, aes(x = elevation, y = log_counts)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
labs(x = "Elevation (metres)", y = "Log number of An. gambiae mosquitoes")
The scatter plot in Figure 3.3 shows a weak negative association, with the average number of mosquitoes tending to decrease as elevation increases. In this case, assuming a linear effect of elevation on the log of the expected mosquito count appears reasonable.
3.1.1.2 When the outcome is an individual-level binary indicator
We now consider the malaria data from Kenya (Section 1.4.3) where the main outcome is the result from a rapid diagnostic test (RDT) for malaria from individuals within households. In this case, because the outcome only takes two values, 1 for a positive RDT test result and 0 otherwise, the direct application of the empirical logit from Equation 3.1 would not help us to generate informative scatter plots. Throughout the book, we will consider the data from the community survey only, hence we work with a subset of the data which we shall name malkenya_comm
# Subset to community surveys only
malkenya_comm <- malkenya[malkenya$Survey == "community", ]To show how this issue can be overcome, let us consider the variables age and gender. To generate a plot that can help us understand between the relationship with malaria prevalence and the two risk factors, we proceed as follows.
Using the cut() function, we first split age (in years) into classes through the argument breaks. The classification of age into \([0,5]\), \((5, 10]\) and \((10, 15]\) is common in many malaria epidemiology studies. For older ages, there is no strong requirement for a particular classification, and alternative groupings are equally plausible; the categories used here are intended purely for illustration.
We then compute the empirical logit, using the total number of cases within each combination of age class and gender. For a given age group and gender, which we denote as \(\mathcal{C}\), the empirical logit in Equation 3.1, now takes the form \[
l_{\mathcal{C}} = \log\left\{\frac{\sum_{i \in \mathcal{C}} y_{i} + 0.5}{|\mathcal{C}|- \sum_{i \in \mathcal{C}} y_{i} + 0.5} \right\}
\tag{3.4}\] where \(y_i\) are the individual binary outcomes and \(i\in \mathcal{C}\) is used to indicate that the sum is carried out over all the individuals who belong the class \(\mathcal{C}\), identified by a specific age group and gender. Finally, \(|\mathcal{C}|\) is the number of individuals who fall within \(\mathcal{C}\). In the code above, the empirical logit in Equation 3.4 is computed using the aggregate() function. An inspection of the object age_class_data, a data frame, shows that the empirical logit is found in the column named RDT.
# Compute the average age by age class and gender
age_mean <- aggregate(Age ~ Age_class + Gender,
data = malkenya_comm,
FUN = mean)
age_class_data$age_mean_point <- age_mean$Age
# Compute the number of individuals by age class and gender
age_n <- aggregate(Age ~ Age_class + Gender,
data = malkenya_comm,
FUN = length)
age_class_data$n_obs <- age_n$AgeIn order to generate the scatter-plot, we compute the average age within each age group by gender, and use these as our values for the x-axis. The resulting Age columns are assigned to age_mean_point and n_obs, respectively.
# Scatterplot of empirical logit RDT vs. age
ggplot(age_class_data, aes(x = age_mean_point, y = RDT)) +
geom_point(aes(size = n_obs, colour = Gender)) +
labs(x = "Age (years)", y = "Empirical logit")
The resulting plot in Figure 3.4 shows the empirical logit against age by gender, with the size of each of the points proportional to the number of observations falling within each class. The observed patterns are explained by the fact that young children, especially those under the age of five, are particularly vulnerable to severe malaria infections. This is primarily due to their immature immune systems and lack of acquired immunity. As individuals grow older, they generally develop partial immunity to malaria through repeated exposure to the disease. This acquired immunity can provide some level of protection against severe malaria. At the same time, gender roles and activities can influence exposure to malaria-carrying mosquitoes. For example, men may spend more time outdoors for work or other activities, increasing their exposure to mosquito bites and thus their risk of infection. In addition, there are also biological factors to consider. Hormonal and genetic differences between males and females may also contribute to variations in immune responses to malaria infection. The interaction between age and gender is complex and may vary depending on the specific context and population being studied. A 2020 report from the Bill & Melinda Gates foundation provides a detailed overview of this and other aspects related to gender and malaria (Katz and Bill & Melinda Gates Foundation 2020).
To account for age in a model for malaria prevalence, several approaches are possible, some of which have been developed using biological models (Smith et al. 2007). To model the patterns observed in Figure 3.4, we can follow the same approach used in the previous section to model the relationship between elevation and river blindness prevalence. First, let us consider age without the effect of gender. Let \(p_{j}(x_i)\) denote the probability of a positive RDT for the \(j\)-th individual living in a household at location \(x_i\). Assuming that malaria risk increases during childhood, reaches a peak around 15 years of age, and then changes more gradually in adulthood, we can capture this non-linear relationship using a linear spline with two knots, at 15 and 40 years. This allows the effect of age to vary across three intervals: childhood, early adulthood, and later adulthood. This model can be expressed as \[ \begin{aligned} \log\left\{\frac{p_{j}(x_i)}{1-p_j(x_i)}\right\} = \beta_{0} + \beta_{1}a_{ij}+\beta_{2} \times\max\{a_{ij}-15, 0\} + \beta_{3}\max\{a_{ij}-40, 0\} \end{aligned} \tag{3.5}\] where \(a_{ij}\) is the age, in years, of the \(j\)-th individual at household \(i\). This formulation defines a piecewise linear relationship between age and the log-odds of infection. In particular, the effect of age is \(\beta_{1}\), for \(a_{ij} < 15\), \(\beta_{1}+\beta_{2}\), for \(15 < a_{ij} < 40\), and \(\beta_{1}+\beta_{2}+\beta_{3}\) for \(a_{ij} > 40\).
Figure 3.4 suggests a potential interaction between age and gender, as the difference in RDT prevalence between males and females appears to vary with age, with wider gaps observed among individuals older than 20 years. To examine this pattern using a standard Binomial regression model, the linear predictor for RDT prevalence can be expressed as \[ \begin{aligned} \log\left\{\frac{p_{j}(x_i)}{1-p_j(x_i)}\right\} = \beta_{0} + (\beta_{1} + \beta_{1}^*g_{ij})\times a_{ij}+(\beta_{2} + \beta_{2}^*g_{ij})\times\max\{a_{ij}-15, 0\} + \\ (\beta_{3} + \beta_{3}^*g_{ij}) \times \max\{a_{ij}-40, 0\} \end{aligned} \tag{3.6}\] where \(g_{ij}\) is the indicator for gender, with 1 corresponding to male and 0 to female. The coefficients \(\beta_{1}^*\), \(\beta_{2}^*\) and \(\beta_{3}^*\) thus quantify the differences in risk between the two genders for ages below 15 years, between 15 and 40 years, and above 40 years, respectively. If all of those coefficients were 0, the model in Equation 3.5 would be recovered.
The code above shows how to fit the model specified in Equation 3.6. The terms Age, pmax(Age-15, 0) and pmax(Age-40, 0) respectively correspond to \(\beta_{1}\), \(\beta_{2}\) and \(\beta_{3}\), whilst the Gender:Age, Gender:pmax(Age-15, 0) and Gender:pmax(Age-40, 0) to \(\beta_{1}^*\), \(\beta_{2}^*\) and \(\beta_{3}^*\), respectively.
| Parameter | Term | Estimate | Std. error | z value | p-value |
|---|---|---|---|---|---|
| \(\beta_0\) | Intercept | -1.058 | 0.102 | -10.331 | 0.000 |
| \(\beta_1\) | Age | -0.034 | 0.013 | -2.584 | 0.010 |
| \(\beta_2\) | Age above 15 years | -0.040 | 0.024 | -1.687 | 0.092 |
| \(\beta_3\) | Age above 40 years | 0.092 | 0.025 | 3.695 | 0.000 |
| \(\beta_1^*\) | Male × age | 0.014 | 0.012 | 1.170 | 0.242 |
| \(\beta_2^*\) | Male × age above 15 years | -0.036 | 0.031 | -1.153 | 0.249 |
| \(\beta_3^*\) | Male × age above 40 years | 0.025 | 0.043 | 0.567 | 0.571 |
The estimated coefficients are reported in Table 3.1. The three interaction coefficients are not individually statistically significant at the 5% level. However, removing the interaction on the basis of these individual p-values would be inappropriate, because the interaction is represented by three coefficients. Instead, the model with and without the interaction terms should be compared using a likelihood ratio test, as shown below.
| Model | Residual df | Residual deviance | Difference in deviance | p-value |
|---|---|---|---|---|
| Reduced model (no interaction) | 3348 | 2675.62 | ||
| Full model (age × sex interaction) | 3345 | 2673.82 | 1.81 | 0.61 |
To assess the null hypothesis that \(\beta_{1}^*=\beta_{2}^*=\beta_{3}^*=0\), we compare the model with and without the age-by-sex interaction terms. This is done by fitting the reduced model under the null hypothesis and performing a likelihood ratio test using the anova() function. The results are reported in Table 3.2. The p-value indicates that we do not find evidence against the null hypothesis, hence in our analysis of the data we might favour the simplified model that does not assume an interaction between the two genders. However, it is important to note that the likelihood used here comes from a standard Binomial model that does not account for residual spatial correlation. As will be shown later in the book, ignoring such correlation typically leads to underestimated standard errors for the regression coefficients and, consequently, p-values that are smaller than they should be. The resulting p-values from the likelihood ratio test may therefore be unreliable. Nonetheless, these exploratory analyses remain useful for identifying patterns that may or may not persist once a more flexible model is fitted. In practice, if an interaction appears weak or non-significant at this stage, one would expect even less evidence for it under a model that captures additional variability. Conversely, when certain relationships are scientifically plausible - such as age-specific differences between genders - it may still be appropriate to retain them regardless of their statistical significance.
The approach just illustrated, can also be applied to explore the association with other continuous variables that are a property of the household and not of the individual. Let us, for example, consider the variable elevation from the malkenya dataset.
Following the same approach used for age, we first split elevation into classes. To define these, we use the deciles of the empirical distribution of elevation which we calculate using the quantile() function above. In this way we also ensure that the number of observations falling within each class of elevation is approximately the same.
# Computation of the empirical logit by classes of elevation
elev_class_data <- aggregate(RDT ~ elevation_class,
data = malkenya_comm,
FUN = function(y) elogit(sum(y), length(y)))
# Compute the average elevation within each elevation class
elevation_mean <- aggregate(elevation ~ elevation_class,
data = malkenya_comm,
FUN = mean)
elev_class_data$elevation_mean <- elevation_mean$elevation
# Compute the number of observations within each elevation class
elevation_n <- aggregate(elevation ~ elevation_class,
data = malkenya_comm,
FUN = length)
elev_class_data$n_obs <- elevation_n$elevationWe then compute the empirical logit, the average elevation and total number of observations for each class of elevation. The empirical logit is computed as already defined in Equation 3.4, where now the definition of \(\mathcal{C}\) is given by a specific decile used to split the distribution of elevation.
# Scatterplot of the empirical logit of RDT vs. mean elevation
ggplot(elev_class_data, aes(x = elevation_mean, y = RDT)) +
geom_point(aes(size = n_obs)) +
labs(x = "Elevation (metres)", y = "Empirical logit")
The resulting plot in Figure 3.5 shows an approximately linear relationship with decreasing values of the empirical logit for increasing elevation. This is expected because the cooler environment at higher altitudes is less favourable to the development of the overall mosquito life cycle.
An alternative approach to generate a scatter plot for assessing the association between elevation and the empirical logit would be to aggregate the data at household level, rather than using classes of elevation. However, this approach cannot be applied when only one individual is sampled for each household. In the case of the malkenya data, the great majority of the households only include one individual making this second approach less useful than the one illustrated.
Other more sophisticated approaches for the exploration of the associations between covariates and binary outcomes are available. For example, the use of the empirical logit could be avoided by using non-parametric regression methods for Binomial outcomes (Bowman and Azzalini 1997), also implemented in the sm package in R. Our view is that a careful exploratory analysis based on simpler methods, such as those illustrated above, can be equally effective to inform the model formulation.
3.1.2 Exploring overdispersion in count data
One of the main advantages in the use of covariates is the ability to attribute part of the variation in the outcome to a set of measured variables and, hence, reduce the uncertainty of our inferences. However, it is almost always the case that the available covariates do not fully explain the variation in the outcome. Unmeasured covariates can contribute to the remaining, or residual, variation. In a standard linear regression model, this variation is quantified directly through the residual variance. For count data, the relationship is less direct. Unmeasured heterogeneity can lead to overdispersion, but this can also arise from clustering, dependence between observations, excess zeros, or misspecification of the mean or probability distribution.
Overdispersion occurs when the variability of the data is larger than that implied by the generalized linear model (GLM) fitted to them. For example, if we consider the Binomial distribution, the presence of overdispersion implies that \(V(Y_i) > m_i \mu_{i}(1-\mu_i)\), where we recall that \(m_i\) is the Binomial denominator and \(\mu_i\) is the probability of “success” for each of the \(m_i\) Bernoulli trials; for a Poisson distribution with \(E(Y_i) = \mu_i\), instead, overdispersion implies that \(V(Y_i) > \mu_{i}\).
Assessment of the overdispersion for count data can be carried out in different ways depending on the goal of the statistical analysis. Since the focus of this book is to illustrate how to formulate and apply geostatistical models, the most natural approach to assess overdispersion is through the use of generalized linear mixed models (GLMMs). The class of GLMMs that we consider in this and the next section are obtained by replacing the spatial Gaussian process \(S(x_i)\) introduced in Equation 1.4 with a set of mutually independent random effects, which we denote as \(Z_i\), and thus write \[ g(\mu_i) = d(x_i)^\top \beta + Z_i. \tag{3.7}\] The model above represents one possible source of overdispersion through \(Z_i\). Its variance quantifies unstructured heterogeneity on the scale of the linear predictor. The resulting amount of overdispersion on the outcome scale also depends on the link function and the mean of the outcome, so it cannot be read directly from the variance of \(Z_i\). To illustrate this, we carry out a small simulation. For simplicity, we consider the Binomial mixed model with an intercept only, hence \[ \log\left\{\frac{\mu_i}{1-\mu_i}\right\} = \beta_0 + Z_i \tag{3.8}\] and assume that the \(Z_i\) follow a set of mutually independent Gaussian variables with mean 0 and variance \(\tau^2\). In our simulation we vary \(\beta_0\) over the set \(\{-3, -2, -1, 0, 1, 2, 3\}\) and set \(\tau^2=0.1\) and the binomial denominators to \(n_i = 100\). For a given value of \(\beta_0\), we then proceed through the following iterative steps.
Simulate 10,000 values for \(Z_i\) from a Gaussian distribution with mean 0 and variance \(\tau^2\).
Compute the probabilities \(\mu_i\) based on Equation 3.8.
Simulate 10,000 values from a Binomial model with probability of success \(\mu_i\) and denominator \(n_i\).
Compute the empirical variance of the counts \(y_i\) simulated in the previous step.
Change the value of \(\beta_0\) and repeat the previous steps, for all the values of \(\beta_0\).
The code below shows the implementation of the above steps in R.
# Reproducible simulation
set.seed(2026)
# Number of simulations
n_sim <- 10000
# Variance of the Z_i
tau2 <- 0.1
# Binomial denominator
bin_denom <- 100
# Intercept values
beta0 <- c(-3, -2, -1, 0, 1, 2, 3)
# Vector where we store the computed variance from
# the simulated counts from the Binomial mixed model
var_data <- rep(NA_real_, length(beta0))
mean_prob <- rep(NA_real_, length(beta0))
for (j in seq_along(beta0)) {
# Simulation of the random effects Z_i
Z_i_sim <- rnorm(n_sim, sd = sqrt(tau2))
# Linear predictor of the Binomial mixed model
lp <- beta0[j] + Z_i_sim
# Probabilities of the Binomial distribution conditional on Z_i
prob_sim <- exp(lp) / (1 + exp(lp))
# Simulation of the counts from the Binomial mixed model
y_i_sim <- rbinom(n_sim, size = bin_denom, prob = prob_sim)
# Empirical variance from the simulated counts
var_data[j] <- var(y_i_sim)
# Marginal probability averaged over the random effects
mean_prob[j] <- mean(prob_sim)
}
# Variance from a standard Binomial model with the same marginal mean
var_binomial <- bin_denom * mean_prob * (1 - mean_prob)# Create a data.frame for plotting the results
plot_df <- data.frame(beta0 = beta0,
Binomial.mixed.model = var_data,
Standard.Binomial.model = var_binomial) |>
pivot_longer(cols = c(Binomial.mixed.model, Standard.Binomial.model),
names_to = "Model",
values_to = "Variance") |>
mutate(Model = factor(Model,
levels = c("Binomial.mixed.model",
"Standard.Binomial.model"),
labels = c("Binomial mixed model",
"Standard Binomial model")))
ggplot(plot_df, aes(x = beta0, y = Variance, color = Model)) +
geom_point() +
geom_line() +
labs(x = expression(beta[0]), y = "Variance") +
scale_color_brewer(type = "q", palette = 2) +
theme_minimal() +
theme(legend.title = element_blank(), legend.position = c(0.2, 0.85))
Figure 3.6 shows the results of the simulation. In this figure, the orange line corresponds to the variance of a standard Binomial model with the same marginal mean as the corresponding mixed model. It is computed as \(n_i \bar{\mu}_i (1-\bar{\mu}_i)\), where \(\bar{\mu}_i\) is approximated by averaging the simulated probabilities over the distribution of \(Z_i\). As expected, this plot shows that the variance of the simulated counts from the mixed model in Equation 3.8 is larger than would be expected under the standard Binomial model. In this example, setting \(\tau^2 = 0.1\) therefore induces overdispersion, although the same value of \(\tau^2\) would not imply the same increase in variance for every value of the mean or under a different link function. The magnitude of \(\tau^2\) can instead be interpreted on the scale of the linear predictor. Following from the properties of a univariate Gaussian distribution, a priori the \(Z_i\) will take values between \(-1.96 \sqrt{\tau^2}\) and \(+1.96 \sqrt{\tau^2}\) with approximately 95\(\%\) probability. That implies that \(\exp\{Z_i\}\), which expresses the effect of the random effects on the odds ratios, will be with 95\(\%\) probability between \(\exp\{-1.96 \sqrt{\tau^2}\}\) and \(\exp\{+1.96 \sqrt{\tau^2}\}\). By replacing \(\tau^2\) with the chosen value for the simulation, these limits become about 0.54 and 1.86. Thus, with approximately 95% probability, \(Z_i\) has a multiplicative effect on the odds between 0.54 and 1.86.
To consolidate the concepts introduced so far, we encourage you to complete Exercises 1 and 2 at the end of this chapter, which further explore the use of generalized linear mixed models as a tool for accounting for overdispersion.
3.1.2.1 Maximum likelihood estimation of generalized linear mixed models for count data
We now illustrate how to fit a generalized linear mixed, using the anopheles dataset as an example. We consider two models: an intercept-only model and one that uses elevation as a covariate. Let \(\mu(x_i)\) be the number of mosquitoes captured at a location \(x_i\); then the linear predictor with elevation as a covariate takes the form \[
\log\{\mu_i\} = \beta_{0} + \beta_{1} d(x_i) + Z_i
\tag{3.9}\] where \(d(x_i)\) indicates the elevation in metres at location \(x_i\) and the \(Z_i\) are independent and identically distributed Gaussian variables with mean 0 and variance \(\tau^2\). The model with an intercept only is simply obtained by setting \(\beta_1 = 0\).
We carry out the estimation in R using the glmer() function from the lme4 package (see Bates et al. (2015) for a detailed tutorial). The glmer() function implements the maximum likelihood estimation for generalized linear mixed models. The code below shows how glmer() is used to carry out this step for the model in Equation 3.9 and the one without covariates.
# Create an ID for the location
anopheles$ID_loc <- seq_len(nrow(anopheles))
# Poisson mixed model with elevation
fit_glmer_elev <- glmer(An.gambiae ~ scale(elevation) + (1 | ID_loc),
family = poisson, data = anopheles, nAGQ = 25)
# Poisson mixed model with intercept only
fit_glmer_int <- glmer(An.gambiae ~ (1 | ID_loc),
family = poisson, data = anopheles, nAGQ = 25)To fit the model with glmer(), we first must create a variable in our dataset that allows us to identify the location associated with each count. In this case, since every row corresponds to a different location, we simply use the row number to identify the locations and save this in a new column called ID_loc. The random effects \(Z_i\) are then included in the model by adding (1 | ID_loc) in the formula of the glmer() function.
When introducing the variable elevation, we standardise it so that its mean is 0 and its variance is 1. This helps the model-fitting algorithm converge more reliably and is generally good practice when covariates are measured on very different scales. Importantly, standardising a variable does not change the fit of the model to the data, it only changes how we interpret the regression coefficients. In other words, the model is mathematically equivalent; only the units of measurement differ. For example, suppose we model mosquito counts using elevation measured in metres and obtain a regression coefficient of \(-0.002\). Under our log-link model, this implies that a one-metre increase in elevation multiplies the expected mosquito count by \(\exp(-0.002) \approx0.998\), corresponding to an approximate decrease of \(100\%\times{\exp(−0.002)−1}\approx−0.2\%\). If we instead standardise elevation so that one unit represents one standard deviation (say 150 m), the corresponding coefficient becomes \(-0.002 \times 150 = -0.30\). Both models yield identical fitted values and predictions. The only difference is that in the standardised model, the coefficient refers to a one-standard-deviation change in elevation rather than a one-metre change. Standardization also changes the interpretation of the intercept. In the unstandardised model, the intercept corresponds to the expected response when elevation is 0 m, whereas in the standardised model it corresponds to the expected response at the mean elevation. Since the latter is often more meaningful, standardisation can also improve the interpretability of the intercept.
The argument nAGQ is used to define the precision of the approximation of the maximum likelihood estimation algorithm. By default nAGQ = 1, which corresponds to the Laplace approximation. Values for nAGQ larger than 1 are used to define the number of points of the adaptive Gaussian-Hermite quadrature. The general principle is that the larger nAGQ the better, but at the expense of an increased computing time. Based on the guidelines and help pages of the lme4 package, it is stated that a reasonable value for nAGQ is 25. For more technical details on this aspect, we refer you to Bates et al. (2015).
We can now summarise the fitted models for the mosquito dataset by reporting the quantities that are most relevant for interpretation, rather than printing the full output returned by summary().
| Term | Estimate | Std. error | z value | p-value | Exponentiated estimate | Interpretation |
|---|---|---|---|---|---|---|
| Intercept | 1.530 | 0.094 | 16.342 | 0.000 | 4.62 | Expected mosquito count at mean elevation |
| Elevation (standardised) | -0.198 | 0.089 | -2.212 | 0.027 | 0.82 | Multiplicative change per 1 SD increase |
| Model | AIC | BIC | Log-likelihood | Random-effect variance | Random-effect SD |
|---|---|---|---|---|---|
| Intercept only | 294.599 | 300.106 | -145.300 | 0.761 | 0.872 |
| Elevation | 291.838 | 300.098 | -142.919 | 0.715 | 0.845 |
The estimated fixed effects for the model including elevation are reported in Table 3.3. The estimated regression coefficient associated with elevation, \(\beta_1\), is negative and statistically significant, suggesting that mosquito abundance decreases as elevation increases. Since the model uses a logarithmic link function, exponentiating the coefficient gives the multiplicative effect on the expected mosquito count. In this case, a one-standard-deviation increase in elevation multiplies the expected mosquito count by \(\exp(-0.19794)=0.82\), corresponding to an approximate reduction of \[ 100 \times \{\exp(-0.19794)-1\} \approx -18\%. \]
Because elevation has been standardised, a one-unit increase corresponds to an increase in the original elevation variable equal to its standard deviation, which for the elevation variable is approximately 100 metres.
The exponentiated intercept reported in Table 3.3 also has a useful interpretation. Since elevation has been standardised to have mean 0, the intercept corresponds to the expected mosquito count at the mean elevation, for a location whose random effect is equal to zero.
The comparison between the intercept-only model and the model including elevation is shown in Table 3.4. From this table, we observe that the estimated variance of the random effect, \(\tau^2\), decreases slightly from 0.761 in the intercept-only model to 0.715 in the model including elevation. Recall that \(\tau^2\) represents the variance of the location-specific random effects \(Z_i\). The quantities reported under Random-effect SD are simply the square roots of the corresponding variances.
As expected, including elevation in the model explains a small proportion of the between-location variability captured by the random effects, though a substantial amount of residual variation remains unexplained. The estimated values of \(\tau^2\) therefore suggest the presence of extra-Poisson variation in the mosquito counts that is not accounted for by elevation alone.
In the next section, we will illustrate how to assess the presence of residual correlation for continuous measurements and overdispersed count data.
3.1.3 Exploring residual spatial correlation
In its most basic form, the concept of spatial correlation can be succinctly encapsulated by Tobler (1970) first law of geography, which posits that “everything is interconnected, but objects in close proximity exhibit stronger relationships than those situated farther apart.” After we have identified the key variables to introduce as covariates in the model (Section 3.1.1) and, in the case of count data, assessed the presence of overdispersion (Section 3.1.2), our final exploratory step consists of assessing whether the residuals of the non-spatial model show evidence of spatial correlation. Hence, in geostatistical modelling, the interest is not in the spatial correlation of the data, but rather in understanding whether the variation in the outcome unexplained by the covariates exhibits spatial correlation. We call this residual spatial correlation, to emphasise that spatial correlation is a concept relative to the covariates that we have introduced in the model.
In the context of geostatistical analysis, the tool that is generally used to assess the residual spatial correlation is the so-called empirical variogram. Before looking at the mathematical definition of the empirical variogram, let us consider a generalized linear mixed model as expressed in Equation 3.7. Our goal is then to question the assumption of independently distributed random effects \(Z_i\) by asking whether the \(Z_i\) show evidence of spatial correlation. Let \(Z_i\) and \(Z_j\) be two random effects that are associated with two different locations \(x_i\) and \(x_j\), respectively, and let us take the squared difference between the two \[ V_{ij} = (Z_i - Z_j)^2. \tag{3.10}\] How does the spatial correlation affect the value of \(V_{ij}\)? To answer this question, we can refer to the aforementioned Tobler’s law of geography. When \(x_i\) and \(x_j\) will be closer to each other, then \(Z_i\) and \(Z_j\) will also tend to be more similar to each other, thus making \(V_{ij}\) smaller, on average. On the contrary, when \(x_i\) and \(x_j\) will be further apart, then \(V_{ij}\) will become larger, on average. We can then construct the empirical variogram by considering all possible pairs of locations \(x_i\) and \(x_j\), for which we then compute \(V_{ij}\) and plot this against the distance between \(x_i\) and \(x_j\), which we denote as \(u_{ij}\). If there is spatial correlation in the random effects \(Z_i\), then this will manifest as an average increase in the \(V_{ij}\) as \(u_{ij}\) increases. However, there are still two issues that we have to address before we can generate and plot the empirical variogram.
The first issue is that we do not observe \(Z_i\) as, by definition, this is a latent variable. Hence, we require an estimate for \(Z_i\) which we can then feed into \(V_{ij}\). To emphasise this point, from now on, we shall replace Equation 3.10 with \[
\hat{V}_{ij} = (\hat{Z}_{i} - \hat{Z}_j)^2.
\tag{3.11}\] Several options are available for estimating \(Z_{i}\). Our choice is to use the mean of the predictive distribution of \(Z_i\), that is, the distribution of \(Z_{i}\) conditional on the observed data. This estimator for \(Z_i\) is also readily available from the output of the lmer() and glmer() functions of the lme4 package, as we will illustrate later in our example in this section.
The second issue is that if we simply plot \(\hat{V}_{ij}\) against the distances \(u_{ij}\) (also known as cloud variogram; see Exercise 3 at the bottom of this chapter), due to the high noisiness in the \(\hat{V}_{ij}\), it may be quite difficult to assess the presence of an increasing trend in the \(\hat{V}_{ij}\) and thus detect spatial correlation. Hence, it is general practice to group the distances \(u_{ij}\) into classes, say \(\mathcal{U}\), and then take average of all the \(\hat{V}_{ij}\) that fall within \(\mathcal{U}\). We can now write the formal definition of the empirical variogram as \[ \hat{V}(\mathcal{U}) = \frac{1}{2 |\mathcal{U}|} \sum_{(i, j): u_{ij} \in \mathcal{U}} \hat{V}_{ij} \tag{3.12}\] where \(|\mathcal{U}|\) denotes the number of pairs of locations that fall within the distance class \(\mathcal{U}\). The rationale behind dividing by 2 in \(1/2 |\mathcal{U}|\) from the above equation, will be elucidated in Section 3.2. When creating the empirical variogram plot, each estimate is placed at the mean distance of the pairs in its distance class \(\mathcal{U}\).
Before we can evaluate residual spatial correlation, there remains one crucial concern: relying solely on a visual inspection of the empirical variogram is susceptible to human subjectivity. Furthermore, it is worth noting that even a seemingly upward trend observed in the empirical variogram might be merely a product of random fluctuations, rather than a reliable indication of actual residual spatial correlation. To address these concerns, we compare the observed empirical variogram with variograms generated under spatial independence. This comparison is an exploratory diagnostic, constructed through the following steps.
Randomly reassign the observed values or fitted residuals among the fixed sampling locations. This breaks any spatial association while preserving the values and the spatial configuration of the data.
Compute the empirical variogram \(\hat{V}(\mathcal{U})\) for the permuted dataset.
Repeat steps 1 and 2 many times. The examples below use 999 permutations; together with the observed variogram, this gives 1,000 exchangeable curves.
Rank the complete curves using extreme-rank-length ordering and form a 95\(\%\) global envelope. Unlike separate intervals at each distance, this construction considers departures anywhere along the curve.
Compare the observed curve with the envelope. Departures from the envelope, particularly a coherent pattern across short distances, suggest that spatial independence may be implausible and warrant further modelling.
The extent and pattern of these deviations can vary considerably. A few departures may reflect weak dependence, random variation or sensitivity to the chosen distance classes, so the whole pattern should be considered rather than treating an individual point as a formal test.
Fitting an explicit geostatistical model provides a formal model-based way to estimate residual spatial dependence and its range. The resulting assessment is conditional on the chosen model and its assumptions, so it should not be regarded as definitive. The empirical variogram should be viewed as an exploratory diagnostic rather than conclusive evidence on its own.
We now show an application of all the concepts introduced in this section to the Liberia data on river blindness.
3.1.3.1 Example: assessing spatial correlation for the Liberia data
We consider the Binomial mixed model that uses the log-transformed elevation as a covariate to model river blindness prevalence, hence \[
\log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_{0} + \beta_{1}\log\{d(x_i)\} + Z_i
\tag{3.13}\] where \(d(x_i)\) is the elevation in metres at location \(x_i\) and the \(Z_i\) are i.i.d. Gaussian variables with mean 0 and variance \(\tau^2\). We first fit the model above using the glmer() function.
# Convert the dataset into an sf object
liberia <- st_as_sf(liberia, coords = c("long", "lat"), crs = 4326)
# Create the ID of the location
liberia$ID_loc <- seq_len(nrow(liberia))
# Binomial mixed model with log-elevation
fit_glmer_lib <- glmer(cbind(npos, ntest - npos) ~ log(elevation) +
(1 | ID_loc),
family = binomial, data = liberia, nAGQ = 25)| Parameter | Term | Estimate | Std. error | z value | p-value |
|---|---|---|---|---|---|
| \(\beta_0\) | Intercept | -3.122 | 0.264 | -11.847 | 0 |
| \(\beta_1\) | Log elevation | 0.332 | 0.051 | 6.463 | 0 |
| \(\tau^2\) | Location-level random-effect variance | 0.085 |
The parameter estimates are reported in Table 3.5. The estimate of \(\tau^2\) quantifies the residual variation between locations after accounting for elevation and represents a moderate level of overdispersion.
liberia$Z_hat <- ranef(fit_glmer_lib)$ID_loc[, 1]We use ranef() to extract the estimates of the random effects \(Z_i\) and save them in the dataset. We then use variogram() from RiskMap to compute the empirical variogram for the estimated \(\hat{Z}_i\).
liberia_variog <- variogram(data = liberia,
variable = "Z_hat",
breaks = c(15, 30, 40, 80, 120,
160, 200, 250, 300, 350),
distance_units = "km",
n_permutations = 999,
seed = 2026)Through the argument breaks we can specify the classes of distance, previously denoted by \(\mathcal{U}\); check the help page of variogram() to see how this is defined by default. The values passed to breaks in the code above define the following distance classes \(\mathcal{U}\): \([15, 30]\), \((30, 40]\) and so forth, with the final class being \((300, 350]\). Pairs separated by less than 15 km or more than 350 km are excluded. The argument n_permutations allows the user to specify the number of permutations used to generate the global envelope for spatial independence. The seed argument makes the permutations reproducible.
Both variogram() and summarise_distance() compute distances in the units selected through distance_units. If the data already use a projected CRS, that CRS is retained and its native linear units are converted as needed. If the data use longitude and latitude, the functions automatically transform them to an appropriate UTM CRS and report the transformation. A different projected CRS can be requested through the distance_crs argument.
summarise_distance(data = liberia, distance_units = "km") min max mean median
3.223757 516.399718 200.155078 186.059958 The summarise_distance() function within the RiskMap package can be used for gauging the extent of the area covered by your dataset, aiding in the selection of appropriate values to be passed to the breaks argument.
In the provided output above, we observe that for the Liberia dataset, the minimum and maximum distances are approximately 3 km and 516 km, respectively. While there is not a one-size-fits-all recommendation for setting breaks, two fundamental principles should inform your choice. Firstly, it is advisable to avoid choosing overly large distance intervals, as the uncertainty associated with the empirical variogram tends to increase with distance due to fewer available pairs of observations for estimation. Secondly, especially when spatial correlation is not strong, it is crucial to carefully explore the behaviour of the variogram at smaller distances. Consequently, it is generally advisable to experiment with different breaks configurations and observe how they impact the pattern of the empirical variogram.
plot_variogram(liberia_variog, plot_envelope = TRUE)
The empirical variogram in Figure 3.7 remains within the global envelope. This exploratory diagnostic therefore does not reveal a clear departure from spatial independence. It does not establish that residual spatial correlation is absent: the result depends on the selected distance classes and on how much information the data contain about the spatial pattern. As we will emphasise throughout this book, the empirical variogram should be used to guide further modelling rather than to draw formal statistical conclusions.
3.1.3.2 Exploring residual spatial correlation with linear Gaussian models
When using a linear model to assess spatial correlation, it is important to distinguish two cases: 1) when the data contain only one observation per location; 2) when more than one observation per location is available. We now consider each of these two scenarios separately.
3.1.3.2.1 One observation per location
To illustrate the use of the variogram in this context, we consider the Galicia data on lead concentration in moss samples. Let \(L_i\) denote the observed lead concentration at location \(x_i\). Because lead concentrations are positive and their distribution is right-skewed, we analyse them on the logarithmic scale and define \(Y_i = \log(L_i)\). The logarithmic transformation reduces the influence of unusually large observations, helps to stabilise the variance, and makes the Gaussian modelling assumption more plausible. The simplest possible model for \(Y_i\) is a standard linear model without covariates, which assumes independence among observations:
\[ Y_i = \beta_0 + U_i, \tag{3.14}\]
where the \(U_i\) are independent Gaussian random variables with mean zero and variance \(\omega^2\). Our goal at this stage is to assess whether this assumption of independence is supported by the data, or whether there is evidence of spatial correlation.
However, measurement error from the sampling device may contribute to the observed variability in \(Y_i\). If such error is large relative to the underlying spatial variation, it can obscure the detection of spatial correlation in the residuals. One might therefore consider introducing an additional, location-specific random effect, say \(Z_i\), with mean zero and variance \(\tau^2\), leading to the model
\[ Y_i = \beta_0 + Z_i + U_i, \tag{3.15}\]
where \(Z_i\) represents spatially structured variation due to unmeasured covariates, and \(U_i\) represents unstructured variation, such as measurement error. However, this model is not identifiable because the two random effects, \(Z_i\) and \(U_i\), both contribute additively to the total variance, making it impossible to separate their effects unless one of two conditions is met:
the measurement error variance \(\omega^2\) is known from external information, or
multiple measurements are available at each location, allowing the two sources of variation to be estimated separately.
We will explore the latter scenario in the next section.
For the Galicia data, we do not know the measurement device precision and we only have one observation per location. However, this does not prevent us from using the variogram based on the residuals from Equation 3.14, while keeping in mind the limitations and uncertainty that are inherent to this exploratory tool as remarked at the end of the last paragraph.
# Fitting of the linear model and extraction of the residuals
lm_fit <- lm(log(lead) ~ 1, data = galicia)
# Extract the raw residuals
galicia$residuals <- lm_fit$residuals
# Compute the variogram, using the residuals from the linear model fit,
# and the 95% global envelope for spatial independence
galicia_variog <- variogram(galicia,
variable = "residuals",
distance_units = "km",
breaks = seq(10, 140, length = 15),
n_permutations = 999,
seed = 2026)
# Plotting the results
plot_variogram(galicia_variog, plot_envelope = TRUE)
In the code above, we first fit the linear model in Equation 3.14 and then extract the residuals from this. Note that for this simple model, the residuals of the model are obtained by simply centering the log-transformed outcome to zero by subtracting its mean. Figure 3.8 shows the empirical variogram and the 95% global envelope for spatial independence. Its coherent departure from the envelope at short distances suggests spatial correlation in the log-transformed lead concentrations.
3.1.3.2.2 More than one observation per location
For the case of more than one observation per location, we shall consider the italy_sim dataset. This dataset contains 10 observations at each of 200 locations. The variable ID_loc is a numeric indicator that can be used to identify the location each observation belongs to. For this dataset, we use the population density, named pop_dens, as a log-transformed covariate; we leave you as an exercise to assess whether it is a reasonable modelling choice.
To specify a non-spatial mixed model for the outcome, we then use two subscripts: \(i\) to identify a given location; \(j\) to identify the \(j\)-th observation for a given location \(i\). We denote as \(Z_i\) the location-specific random effect and as \(U_{ij}\) the random variation due to the measurement error inherent to each observation. Hence, we write
\[ Y_{ij} = \beta_0 + \beta_1 \log\{d(x_i)\} + Z_i + U_{ij}, \tag{3.16}\]
where \(d(x_i)\) is the population density at location \(x_i\); as before, we use \(\tau^2\) and \(\omega^2\) to denote the variances of \(Z_i\) and \(U_{ij}\), respectively. To fit this model to the data we use the lmer() function from the lme4 package.
# Fitting a linear mixed model to the italy_sim dataset
# See main text for model specification
lmer_fit <- lmer(y ~ log(pop_dens) + (1 | ID_loc), data = italy_sim)| Parameter | Term | Estimate | Std. error | t value |
|---|---|---|---|---|
| \(\beta_0\) | Intercept | -0.404 | 0.576 | -0.702 |
| \(\beta_1\) | Log population density | 1.338 | 0.076 | 17.506 |
| \(\tau^2\) | Location-level random-effect variance | 2.501 | ||
| \(\omega^2\) | Residual variance | 0.196 |
In the model above, the introduction of the \(Z_i\) random effect is specified in the lmer() function through (1 | ID_loc) in the formula; recall that ID_loc is the numerical indicator that identifies each location. The estimates of the regression coefficients \(\beta_0\) and \(\beta_1\), together with the variance components \(\tau^2\) and \(\omega^2\), are reported in Table 3.6. The location-level random-effect variance \(\tau^2\) is associated with ID_loc, while \(\omega^2\) is the residual variance.
We then extract the estimated location-specific random effects \(Z_i\) and use them to compute the empirical variogram. The code for carrying out these steps is shown below.
# Incorporating the estimated random effects into the data
italy_sim$rand_eff <- ranef(lmer_fit)$ID_loc[italy_sim$ID_loc, 1]
# Compute the variogram, using the random effects from the linear mixed model
# fit, and the 95% global envelope for spatial independence
italy_sim_variog <- variogram(italy_sim,
variable = "rand_eff",
distance_units = "km",
n_permutations = 999,
seed = 2026)
# Plot the variogram
plot_variogram(italy_sim_variog, plot_envelope = TRUE)
In Figure 3.9, the empirical variogram departs from the global envelope over several neighbouring distance classes. This exploratory pattern suggests residual spatial correlation and motivates fitting a geostatistical model.
3.2 The linear geostatistical model
In this section, we focus on spatially referenced outcomes \(Y_i\) that are continuous. We first consider the simpler case of a single measurement \(Y_i\) per location \(x_i\). Recalling the class of generalized linear models introduced in Section 1.5, the linear predictor takes the form
\[ \mu_{i} = d(x_i)^\top \beta + S(x_i). \tag{3.17}\]
Hence, in this case, we interpret \(\beta\) as the effect on \(\mu_i\) for a unit increase in \(d(x_i)\). Let \(U_i\) denote i.i.d. random variables representing the measurement error associated with \(Y_i\), each having mean zero and variance \(\omega^2\). Thanks to the linear properties of Gaussian random variables, we can also express the linear model in a compact expression, as
\[ Y_i = \mu_i + U_i = d(x_i)^\top \beta + S(x_i) + U_i. \tag{3.18}\]
To fully specify a geostatistical model for our dataset, we must address two critical aspects.
Defining the relationship between each covariate and the mean value \(\mu_i\).
Selecting an appropriate correlation function for \(S(x)\).
As illustrated in the previous sections, the initial step of exploratory analysis allows us to handle the first aspect, where we determine the regression relationship between covariates and the mean value \(\mu_i\).
However, based on existing methods of exploratory analysis, it is difficult to understand what is a suitable correlation function at this stage. A commonly recommended starting point is the Matern correlation function (see Equation 1.5), which offers considerable flexibility in capturing a wide range of correlation structures, under the assumption of stationarity and isotropy. As we shall illustrate in the next example, even estimating a Matern correlation function is a task that poses many inferential challenges, especially due to the poor identifiability of its smoothness parameter \(\kappa\).
3.2.1 Evaluating the inclusion of the measurement error term \(U_i\) and the specification of the smoothness parameter \(\kappa\)
In this section we analyse the Galicia data using a linear geostatistical model for the log-transformed lead concentration, recalling that \(Y_i = \log(L_i)\). Since we do not use covariates, we then write the model as
\[ Y_i = \mu + S(x_i) + U_i \tag{3.19}\]
where were \(S(x)\) is a Matern process with variance \(\sigma^2\), scale parameter \(\phi\) and smoothness parameter \(\kappa\); the \(U_i\) correspond to the measurement error term and denote with \(\omega^2\) their variance.
We carry out the parameter estimation of the model using the glgpm() function from the RiskMap package. This function implements maximum likelihood estimation for generalized linear mixed models using a Matern correlation function, while fixing the smoothness parameter \(\kappa\) at a prespecified value. The object passed to the data argument must be an sf object with a coordinate reference system. The galicia dataset supplied by RiskMap already satisfies these requirements.
# Parameter estimation using the geometry and CRS stored in `galicia`
fit_galicia <- glgpm(log(lead) ~ gp(kappa = 1.5),
data = galicia, family = "gaussian",
distance_units = "km", messages = FALSE)The Gaussian process \(S(x)\) is specified through gp() in the formula. The spatial coordinates and CRS are read directly from the sf geometry, while kappa fixes the smoothness parameter at 1.5. By default, nugget = FALSE and the nugget term (denoted in this book by \(Z_i\)) is excluded. Set nugget = TRUE to estimate its variance, or pass a positive numeric value to fix it. Estimating the nugget for a Gaussian model with only one observation per location is not identifiable unless the measurement-error variance is fixed through fix_var_me; Exercise 7 returns to this point. The distance_units argument controls whether spatial distances, and therefore the scale of \(\phi\), are expressed in kilometres or metres.
summary(fit_galicia)Call:
glgpm(formula = log(lead) ~ gp(kappa = 1.5), data = galicia,
family = "gaussian", distance_units = "km", messages = FALSE)
Linear geostatistical model
Link: identity
Inverse link function = x
'Lower limit' and 'Upper limit' refer to 95% confidence intervals
Regression coefficients
Estimate Lower limit Upper limit StdErr z.value p.value
(Intercept) 0.707418 0.552762 0.862075 0.078908 8.9651 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Estimate Lower limit Upper limit
Measurement error var. 0.01546355 0.00019466 1.2284
Spatial Gaussian process
Matern covariance parameters (kappa = 1.5)
Estimate Lower limit Upper limit
Spatial process var. 0.17127 0.10333 0.2839
Spatial corr. scale (km) 9.02085 5.38424 15.1137
Variance of the nugget effect fixed at 0
Log-likelihood: 69.03029
AIC: -130.0606We can then inspect the point and interval estimates for the model through the summary() function of the model as shown in the code chunk above. This output is presented in three sections: in the first section, we have the results for the regression coefficients; in the second section, we have the estimate for the variance of the measurement error component, \(\omega^2\), whose point estimate is about 0.015; in the final section, we have the estimates for the parameters of the spatial covariance function, \(\sigma^2\) and \(\phi\), which are about 0.171 and 9.021 (km), respectively. The message Variance of the nugget effect fixed at 0 indicates that the nugget has not been included in the model.
At this point, you may be wondering why we have used 1.5 for the value of \(\kappa\) and whether there is a statistical approach for finding a suitable value. We show such an approach in the code below.
However, before examining the code, we would like to point out an important aspect that relates to how the value of \(\kappa\) can affect the estimate of the measurement error variance \(\omega^2\). As we have shown, in Section 1.5.1, values of \(\kappa\) that are closer to zero will give a rougher and less regular surface for \(S(x)\). In the case of a single observation, when \(\kappa\) is closer to zero, it will thus become increasingly difficult to estimate \(\omega^2\) since the likelihood function may attribute most of the noisiness to \(S(x)\) and rely less on \(U_i\) to explain the unstructured random variation found in the data. We can therefore expect smaller estimates of \(\omega^2\) when \(\kappa\) is closer to zero and larger estimates when \(\kappa\) increases. The results generated in the code below clearly illustrate this point.
# Number of the values chosen for kappa
n_kappa <- 10
# Set of values for kappa
kappa_values <- seq(0.5, 3.5, length = n_kappa)
# Vector that will store the values of the likelihood function
# evaluated at the maximum likelihood estimate
llik_values <- rep(NA, length = n_kappa)
# Vector that will store the maximum likelihood estimates
# of the variance of the measurement error
sigma2_me_hat <- rep(NA, length = n_kappa)
# List that will contain all the geostatistical models fitted for
# the different values of kappa specified in kappa_values
fit_galicia_list <- vector("list", n_kappa)
for (i in seq_len(n_kappa)) {
formula_i <- as.formula(sprintf("log(lead) ~ gp(kappa = %s)",
kappa_values[i]))
fit_galicia_list[[i]] <- glgpm(formula_i,
data = galicia, family = "gaussian",
distance_units = "km", messages = FALSE)
llik_values[i] <- fit_galicia_list[[i]]$log_lik
sigma2_me_hat[i] <- coef(fit_galicia_list[[i]])$sigma2_me
}
Panel (a) of Figure 3.10 shows that smaller values of \(\kappa\) lead to smaller point estimates of \(\omega^2\), and vice versa. This raises the question of which value should be chosen for \(\kappa\).
To answer this question, a natural approach is to estimate the model for different values of \(\kappa\) and see which one gives the best fit to the data, according to the likelihood function. In panel (b) of Figure 3.10, we show the profile log-likelihood for 10 values of \(\kappa\) between 0.5 and 3.5. The horizontal dashed line gives the approximate 95% cutoff obtained by subtracting half of the 0.95 quantile of a \(\chi^2\) distribution with one degree of freedom from the maximum profile log-likelihood. In R, this cutoff is obtained as
All the values of \(\kappa\) considered lie above the horizontal line, so the profile log-likelihood is essentially flat over this range. The data provide little information to distinguish between these values, which clearly illustrates the weak identifiability of \(\kappa\). We therefore fix \(\kappa=1.5\) as a pragmatic choice that the data do not contradict, rather than one selected by the likelihood. The difficulty of estimating \(\kappa\) is well documented in geostatistics (Zhang (2004)), which supports this pragmatic approach. The problem is further exacerbated for count data, which tend to be less informative about the correlation structure than continuously measured data.
3.2.2 Modelling hierarchical geostatistical data using the re() function
We now consider the analysis of geostatistical data with a hierarchical structure and show how to formulate and fit a geostatistical model that accounts for the effects of the different layers of the data. For this purpose, we use the italy_sim where each of the sampled locations can be grouped according to two administrative subdivisions of Italy, regions (Admin level 2) and provinces (Admin level 3), as shown in Figure 3.11.
library(geobounds)
italy_regions <- gb_get_adm2("Italy")
italy_provinces <- gb_get_adm3("Italy")# Map of the data with the region boundaries
map_regions <- ggplot() +
geom_sf(data = italy_sim, pch = 4, color = "red") +
geom_sf(data = italy_regions, fill = NA) +
theme_void() +
labs(title = "Regions") +
theme(plot.title = element_text(hjust = 1 / 2))
# Map of the data with the province boundaries
map_provinces <- ggplot() +
geom_sf(data = italy_sim, pch = 4, color = "red") +
geom_sf(data = italy_provinces, fill = NA) +
theme_void() +
labs(title = "Provinces") +
theme(plot.title = element_text(hjust = 1 / 2))
italy_sim (red crosses), with the boundaries of the regions (left) and provinces (right) of Italy.
Let \(V_h\), for \(h=1,\ldots,n_V\) and \(Z_k\), for \(k=1,\ldots,n_{Z}\), be a set of mutually independent Gaussian variables with zero means and variances \(\sigma^2_{V}\) and \(\sigma^2_{Z}\), respectively. Finally let \(\mathcal{A}_h\), for \(h=1,\ldots,n_V\) and \(\mathcal{B}_k\), for \(k=1,\ldots,n_{Z}\), denote the areas encompassed by the boundaries of the region and provinces, respectively. Since we have 10 repeated observations at each locations, we shall use \(Y_{ij}\) to denote the \(j\)-th outcome (found in the data under the column y) at the \(i\)-th location \(x_i\). The model from which we generated data \(y_{ij}\) is
\[ \begin{aligned} Y_{ij} = \beta_0 &+ \overbrace{\beta_1 d(x_i) + S(x_i)}^\text{Location-level effect} + \underbrace{\sum_{h=1}^{n_{V}} I(x_i \in \mathcal{A}_h) \times V_h}_\text{Region effect} \\ &+ \overbrace{\sum_{k=1}^{n_{Z}} I(x_i \in \mathcal{B}_k) \times Z_k}^\text{Province effect} + \underbrace{U_{ij}}_\text{Measurement error} \end{aligned} \tag{3.20}\]
for \(j=1,\ldots 10\) and \(i=1,\ldots,200\), where, \(d(x_i)\) is the log-transformed population density, \(I(x \in \mathcal{R})\) is an indicator function that takes value 1 if the location \(x\) falls within the boundaries of the area denoted by \(\mathcal{R}\) and 0 otherwise. The covariance function chosen in the generation of the data was an exponential correlation function. Hence,
\[ \operatorname{Cov}\{S(x), S(x')\} = \sigma^2 \exp\{-\lVert x-x'\rVert/\phi\}. \]
The random effects \(V_h\), for the region, and \(Z_k\), for the province, can be included by using the re() function in the formula passed to glgpm(), as shown below.
# Location-level, region, and province effects
italy_fit <- glgpm(y ~ log(pop_dens) + gp(kappa = 0.5, nugget = FALSE) +
re(region, province),
data = italy_sim, distance_units = "km",
family = "gaussian")In the code above, the factor variables region and province found in italy_sim are passed to re() and estimated as unstructured random effects, as defined by Equation 3.20. The to_table() function from RiskMap produces the presentation-ready parameter and interval estimates reported in Table 3.7.
to_table(italy_fit, digits = 3)| Estimate | Lower limit | Upper limit | |
|---|---|---|---|
| (Intercept) | -1.074 | -2.032 | -0.117 |
| log(pop_dens) | 1.373 | 1.331 | 1.415 |
| Spatial process var. | 0.917 | 0.250 | 3.362 |
| Spatial corr. scale (km) | 214.520 | 56.396 | 815.994 |
| region (random eff. var.) | 0.241 | 0.012 | 4.904 |
| province (random eff. var.) | 1.457 | 0.989 | 2.146 |
| Measurement error var. | 0.196 | 0.184 | 0.209 |
3.3 Generalized linear geostatistical models
We now consider the modelling of count outcomes, focusing our attention on the Binomial and Poisson geostatistical models.
The use of a Binomial geostatistical model is appropriate whenever the reported counts have a known upper bound. For example, when analysing aggregated disease counts from cross-sectional survey data, the total number of people sampled at a location defines the largest number of possible reported disease cases. Recalling the model specification of Table 1.3, we specify the Binomial geostatistical model by stating that the reported disease counts \(Y_i\) out of \(n_i\) (the maximum value that \(Y_i\) can take) at location \(x_i\), conditionally on the spatial random effects \(S(x_i)\), follow a Binomial distribution with probability \(p(x_i)\) and logit-linear predictor
\[ \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = d(x_i)^\top \beta + S(x_i). \tag{3.21}\]
In this model, \(p(x_i)\) is commonly understood as representing the prevalence of the disease under consideration, defined as the proportion of individuals with the disease at a specific moment. While \(p(x_i)\) is not strictly a proportion, its interpretation as disease prevalence is justified by the fact that if we were to sample a sufficiently large population at location \(x_i\), the proportion of individuals with the disease would converge closely to the value of \(p(x_i)\).
A Binomial geostatistical model is also used when analysing data at individual level where the outcome is binary and usually indicates the positive or negative result of a test for a disease under investigation. In this analysis, a mix of individual-level variables and location-level variables could be introduced to model the probability of a positive test. To formally write the model, we now use \(Y_{ij}\) to indicate the binary outcome taking value 1, if the test was positive, and 0, if negative, for the \(j\)-th individual sampled at location \(x_i\). The linear predictor now takes the form
\[ \log\left\{\frac{p_{ij}}{1-p_{ij}}\right\} = e_{ij}^\top \gamma + d(x_{i})^\top \beta + S(x_i). \tag{3.22}\]
In the model above we have distinguished between covariates \(e_{ij}\), that express the effect on \(p_{ij}\) due to individual traits (such as age and gender), and covariates \(d(x_i)\), that instead quantify the effect on \(p_{ij}\) resulting from the environmental attributes of location \(x_i\) (such as elevation, temperature, and precipitation). In this model, \(p_{ij}\) is the probability of a positive test and can once again be interpreted as the disease prevalence for a population characterised by the individual traits \(e_{ij}\) and residing at location \(x_i\).
Finally, when the counts arise from a population whose size is unknown or when the occurrence of a disease is very small relative to the size of the target population, one should consider the use of a Poisson geostatistical model. We then use \(Y_i\) to denote the reported counts at location \(x_i\) and assume that this, conditionally on a Gaussian process \(S(x_i)\), follows a Poisson distribution with mean \(\lambda(x_i)\) and log-linear predictor
\[ \log\left\{\lambda(x_i)\right\} = d(x_i)^\top \beta + S(x_i). \tag{3.23}\]
For example, if \(Y_i\) corresponds to the number of new cases reported at a health facility at location \(x_i\) over a given period, \(\lambda(x_i)\) is interpreted as the expected number of reported cases over that period. If the size of the population at risk is available, it can be incorporated as an exposure. Denoting the exposure by \(m_i\), we specify the mean of \(Y_i\) as \(m_i\lambda(x_i)\), or equivalently add \(\log(m_i)\) as an offset to the linear predictor in Equation 3.23. The quantity \(\lambda(x_i)\) can then be interpreted as an incidence rate per unit of exposure.
Poisson geostatistical models are also used to map the abundance of disease vectors. In this setting, an exposure can often be defined from the sampling effort, such as the number of trap nights, the area sampled, or the duration of sampling. Including this exposure allows \(\lambda(x_i)\) to be interpreted as the expected number of vectors per unit of sampling effort. When comparable effort information is unavailable, the model instead describes the expected observed count under the sampling design used.
In all three models, namely Equation 3.21, Equation 3.22 and Equation 3.23, one might consider including an additional independent random effect \(Z_i\) to represent unstructured heterogeneity on the link scale. Such heterogeneity can induce overdispersion on the outcome scale, but its interpretation depends on the model and it does not identify the underlying cause. We shall further elaborate on this point in the following examples.
3.3.1 Example: river blindness in Liberia
In section Section 3.1.1.1, we carried out the exploratory analysis of the river blindness prevalence data to assess how elevation could be introduced into the model. In Section 3.1.3.1, we then tested for residual spatial correlation using the empirical variogram and log-transformed elevation as a covariate. Based on these results, we can now formulate a Binomial geostatistical model which we then fit to the data. Hence, we extend the model presented in Equation 3.13 through the inclusion of a spatial Gaussian process, \(S(x)\), and write \[ \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_{0} + \beta_{1}\log\{e(x_i)\} + S(x_i). \tag{3.24}\]
We shall assume that \(S(x_i)\) is a stationary isotropic Gaussian process with an exponential correlation function \(\rho(u) = \exp\{-u/\phi\}\), with scale parameter \(\phi\).
We apply the Monte Carlo Maximum Likelihood (MCML) method for parameter estimation, which is implemented via the glgpm() function. In this section, we set aside the technical details of the estimation process to focus exclusively on interpreting the results. A more detailed explanation of the MCML theory can be found in Section 3.5.2, with a concise summary of the practical aspects provided in Section 3.5.2.3.
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")to_table(fit_liberia, digits = 3)| Estimate | Lower limit | Upper limit | |
|---|---|---|---|
| (Intercept) | -3.525 | -3.997 | -3.054 |
| log(elevation) | 0.399 | 0.308 | 0.490 |
| Spatial process var. | 0.242 | 0.167 | 0.352 |
| Spatial corr. scale (km) | 28.174 | 16.169 | 49.092 |
When using the glgpm() function to fit a Binomial model, we need to specify the Binomial denominators, i.e. the number of people tested at a location in this case, through the argument denominator. Other arguments of glgpm() that are used for carrying out MCML include control_mcmc, which receives the settings created by set_control_mcmc(). Initial values for the model parameters can be supplied through start_pars; when it is omitted, RiskMap constructs default starting values. Here, return_samples = TRUE retains the MCMC samples so that they can be examined later with plot_mcmc(). The seed recorded in set_control_mcmc() makes those samples reproducible. For more information on the default values of these arguments, see the help page for glgpm(). To understand how the MCMC settings can be changed to improve the MCML estimation, see Section 3.5.2.1. The model fitted here does not include a nugget effect. This agrees with the default nugget = FALSE setting of gp().
In the above summary of the model, the interpretation of the regression coefficients is the same as in a standard logistic regression model. Hence, by taking the exponential of the estimate for the regression coefficient \(\beta_1\), we can express the effect for a unit increase in the covariate, which in this case is the log-transformed elevation, on the odds of disease prevalence. Finally, the interpretation of estimates of the covariance parameters is not affected by the type of distribution used to model the main outcome. The estimate of \(\sigma^2\) quantifies the residual spatial variation on the log-odds scale, while the estimate of \(\phi\) indicates relatively long-range spatial correlation, with a practical range close to (\(58.67 \times 3 \approx\)) 180 km.
3.3.2 Example: Anopheles gambiae mosquitoes in Cameroon
In Section 3.1.1.1, we explored the association between elevation and the number of mosquitoes and found that the data showed that on average, the number of Anopheles gambiae (An. gambiae) decreases with increasing elevation. We now fit a log-linear Poisson geostatistical model, with linear predictor
\[ \log\{\lambda(x_i)\} = \beta_0 + \beta_1 d(x_i) + S(x_i) \tag{3.25}\]
where \(d(x_i)\) corresponds to the elevation in metres at location \(x_i\). The anopheles dataset is provided as an sf object. We can then use the following script to fit the model in RiskMap.
an_fit <- glgpm(An.gambiae ~ elevation + gp(),
data = anopheles,
family = "poisson",
control_mcmc = set_control_mcmc(seed = 2026),
messages = FALSE)to_table(an_fit, digits = 3)| Estimate | Lower limit | Upper limit | |
|---|---|---|---|
| (Intercept) | 3.566 | 3.161 | 3.970 |
| elevation | -0.003 | -0.003 | -0.002 |
| Spatial process var. | 0.668 | 0.486 | 0.917 |
| Spatial corr. scale (km) | 4.311 | 2.632 | 7.061 |
The estimated regression coefficient \(\hat{\beta}_1\) is negative, which suggests that the expected number of mosquitoes decreases as elevation increases. For an increase of 100 metres in elevation, the expected count is multiplied by \(\exp(100\hat{\beta}_1)\). Conversely, a decrease of 100 metres multiplies the expected count by \(\exp(-100\hat{\beta}_1)\). The covariance-parameter estimates suggest the presence of residual spatial correlation, albeit on a rather small scale.
summarise_distance(anopheles, distance_units = "km") min max mean median
0.3840312 161.7342628 54.9685731 43.9628246 Using the summarise_distance() function, we can see that, among the locations of the dataset, the average distance is about 55 km and the minimum is about 384 metres. The estimate of \(\phi\) of 4 km is thus relatively small compared to the distances between locations in the data. In Chapter 4, we will further investigate to what extent this small-scale spatial correlation can affect the predictions of the An. Gambiae counts.
3.3.3 Using parametric bootstrap for computing confidence intervals
The bootstrap is a widely used resampling technique introduced by Efron (1979) as a general approach to estimating the variability of statistical estimates without relying on strong parametric assumptions.
While the classical non-parametric bootstrap resamples directly from the observed data, parametric bootstrap instead relies on simulating new datasets from an assumed parametric model, using plug-in estimates obtained from the original fit of the model to the data. This simulation-based method is particularly useful for complex models, such as those considered in this book, that yield an intractable likelihood function.
We now describe how to carry out parametric bootstrap for geostatistical models in RiskMap. The use of simulation-based approaches will be more extensively covered in Chapter 4, but, for the purpose of this section, it is first important to understand how data can be simulated from a geostatistical model.
The way we have formulated geostatistical models in Section 1.5 provides a hint as to what step should be followed when simulating from a geostatistical model. Using the same notation of Section 1.5, these are the steps.
- For a given set of predefined locations \(X\) compute the covariance matrix \(\Sigma\) with \((i,j)\)-th entry \(\sigma^2 \rho(u_{ij})\).
- Compute the covariate 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 from this step as \(S_{(j)}(x)\), for \(j = 1, \ldots, B\).
- Compute the linear predictor at each location of \(X\), by adding the covariate effects from step 2 to the simulated random effects from step 3, hence \(g\{\mu_{(j)}(x)\} = d(x)^\top \beta + S_{(j)}(x)\).
Note that before performing the steps above, the values of the model parameters must be fixed. For example, if simulating from a model fitted to the data, the values of the model parameters should be the maximum likelihood estimates. If the goal is to simulate actual observations for the outcome \(Y\) at the locations \(X\), then we need to add the following steps.
- If the nugget term is included in the model, draw samples \(Z_{(j)}\) for each location in \(X\) from a Gaussian distribution with mean 0 and variance \(\tau^2\). Add the simulated values \(Z_{(j)}\) to \(g\{\mu_{(j)}(x)\}\).
- Compute \(\mu_{(j)}(x)\) by applying the inverse link function \(g^{-1}(\cdot)\) to the samples of the linear predictor previously computed.
- Simulate from the distribution of \(Y\) given \(\mu_{(j)}(x)\), according to the chosen model, for \(j = 1, \ldots, B\).
These steps have been implemented in the simulate_glgpm() function. To show how to perform parametric bootstrap in RiskMap, let us consider the Liberia data on river blindness and fit the following geostatistical model.
fit_liberia_no_nugget <- glgpm(npos ~ log(elevation) + gp(),
denominator = ntest, data = liberia,
model_crs = 32629,
control_mcmc = set_control_mcmc(seed = 2026),
family = "binomial")Using the fitted Binomial geostatistical above, we can simulate a single dataset as follows:
liberia_sim <- simulate_glgpm(fit_liberia_no_nugget,
nsim = 1,
what = "data",
seed = 2026)
sim_data <- simulated_data(liberia_sim, simulation = 1)The fitted parameters are held fixed while a new spatial effect and a new outcome are generated at the observed locations. The simulated_data() function extracts one simulated dataset and retains the original covariates, denominator and geometry. We can therefore refit the same model to sim_data as follows.
fit_liberia_no_nugget_sim <- glgpm(formula = npos ~ log(elevation) + gp(),
data = sim_data,
family = "binomial",
denominator = ntest,
control_mcmc = set_control_mcmc(seed = 2027))The simulated dataset uses the modelling CRS and distance units stored in the fitted model. This ensures that the covariance parameters retain the same spatial interpretation when the model is refitted.
This idea can be extended naturally to perform parametric bootstrap. The procedure involves simulating multiple datasets from the fitted model, refitting the model to each of them, and extracting the estimates of interest. For instance, the following code simulates 100 datasets and fits the same model to each, storing the estimated regression coefficients:
# Simulate 100 datasets
n_sim <- 100
liberia_boot <- simulate_glgpm(fit_liberia_no_nugget,
nsim = n_sim,
what = "data",
seed = 2026)
# List where we will store the estimates from each of the simulated datasets
par_hat <- vector("list", n_sim)
for (i in seq_len(n_sim)) {
# Extract one simulated dataset
sim_data <- simulated_data(liberia_boot, simulation = i)
# Fitting the geostatistical model
fit_sim <- glgpm(formula = npos ~ log(elevation) + gp(),
data = sim_data,
family = "binomial",
start_pars = coef(fit_liberia_no_nugget),
denominator = ntest,
control_mcmc = set_control_mcmc(seed = 2026 + i),
messages = FALSE)
# Extracting the parameter estimates
par_hat[[i]] <- coef(fit_sim)
}After obtaining the parameter estimates, we extract them and rearrange them in a matrix format to facilitate the computation of any summary, including the confidence intervals as shown in the code below.
# Extract all parameters maintaining beta as a vector
all_params <- sapply(seq_len(n_sim), function(i) {
c(par_hat[[i]]$beta,
par_hat[[i]]$sigma2,
par_hat[[i]]$phi)
})
CI <- t(apply(all_params, 1, function(x) quantile(x, c(0.025, 0.975))))
# Combine the Wald and bootstrap intervals for comparison
fit_summary <- summary(fit_liberia_no_nugget)
wald_CI <- rbind(fit_summary$reg_coef[, c("Lower limit", "Upper limit")],
fit_summary$sp[, c("Lower limit", "Upper limit")])
parameter_names <- rownames(wald_CI)
rownames(CI) <- parameter_names
interval_comparison <- data.frame(Parameter = parameter_names,
`Wald lower` = wald_CI[, "Lower limit"],
`Wald upper` = wald_CI[, "Upper limit"],
`Bootstrap lower` = CI[, "2.5%"],
`Bootstrap upper` = CI[, "97.5%"],
check.names = FALSE)| Parameter | Wald lower | Wald upper | Bootstrap lower | Bootstrap upper |
|---|---|---|---|---|
| (Intercept) | -3.997 | -3.054 | -4.205 | -2.688 |
| log(elevation) | 0.308 | 0.490 | 0.232 | 0.537 |
| Spatial process var. | 0.167 | 0.352 | 0.134 | 0.389 |
| Spatial corr. scale (km) | 16.169 | 49.092 | 7.247 | 55.164 |
By comparing the Wald confidence intervals with those derived from the bootstrap approach in Table 3.10, we observe that the former are considerably narrower for the covariance parameters, while the intervals are fairly similar for the regression coefficients.
This discrepancy arises because the Wald-type confidence intervals from summary() rely on a Gaussian approximation of the maximum likelihood estimator for the covariance parameters, which may not be entirely reliable. RiskMap improves the accuracy of this approximation by applying it to the log-transformed parameters and then transforming the confidence intervals back to the original scale. On the other hand, bootstrap-based confidence intervals can better represent finite-sample asymmetry in the estimator, although they remain conditional on the fitted model being an adequate data-generating mechanism. They are also more computationally demanding. The use of only 100 bootstrap samples, as shown above, may be insufficient. It is generally recommended to use at least 1000 iterations, if computationally feasible, to reduce the Monte Carlo error arising from the bootstrap simulations.
3.4 The residuals of a geostatistical model: issues to consider
In geostatistical models, the concept of residuals, which are widely used in regression modelling to carry out model diagnostics, becomes more nuanced than in standard regression models (i.e., models that do not include random effects). Generally, residuals aim to quantify the discrepancy between observed outcomes and what the model predicts.
However, in GLGMs, there are multiple valid definitions of residuals. If we define the residuals as the component of the model that expresses the variation in the outcome not explained by the covariates, then we are effectively referring to the random effects component of the linear predictor, namely \(S(x_i)\), or \(S(x_i) + Z_i\) if a nugget effect \(Z_i\) is included. We refer to these as random effects residuals.
Let \(W_i\) denote the sum of all the random effects in the model, whether \(W_i = S(x_i)\), \(W_i = S(x_i) + Z_i\), or any other random effects structure. In addition to RERs, other types of residuals more common in standard regression modelling can also be used, however, because of the presence of random effects \(W_i\), their definition becomes more convoluted.
Let us consider the Pearson residuals. We define these using the predictive mean and variance as
\[ \operatorname{PR}_i = \frac{y_i - \widehat{\operatorname{E}}[Y_i^*]} {\sqrt{\widehat{\operatorname{Var}}[Y_i^*]}} \]
where \(\widehat{\operatorname{E}}[Y_i^*]\) and \(\widehat{\operatorname{Var}}[Y_i^*]\) are estimates of the predictive expectation and variance of a replicate \(Y_i^*\) given the observed data, with the model parameters fixed at their maximum likelihood estimates. Using Monte Carlo samples \(W_{i,(j)}\), for \(j = 1, \ldots, B\), from the distribution of \(W_i\) given the data (see Section 3.5.2.2 for how these samples are generated in RiskMap), let
\[ m_{ij} = \widehat{\operatorname{E}}[Y_i^* \mid W_{i,(j)}], \qquad v_{ij} = \widehat{\operatorname{Var}}[Y_i^* \mid W_{i,(j)}] \]
and \(\bar m_i = B^{-1}\sum_{j=1}^B m_{ij}\). The predictive moments can then be estimated as
\[ \widehat{\operatorname{E}}[Y_i^*] = \bar m_i, \qquad \widehat{\operatorname{Var}}[Y_i^*] = \frac{1}{B}\sum_{j=1}^B v_{ij} + \frac{1}{B-1}\sum_{j=1}^B(m_{ij}-\bar m_i)^2. \]
The two terms in the variance estimator follow from the law of total variance. The first averages the conditional variances, while the second accounts for variation in the conditional means across the sampled random effects.
An alternative, especially useful when \(Y_i\) is binary, is the deviance residual: \[ {\rm DR}_i = \text{sign}(y_i - \hat{y}_i) \sqrt{2 \left( \ell_i^{S} - \ell_i^{F} \right)} \] where \(\ell_i^{S}\) is the log-likelihood of the so-called saturated model, and \(\ell_i^{F}\) is the log-likelihood of the fitted model evaluated at \(\widehat W_i = \operatorname{E}(W_i \mid y)\). The saturated model is a hypothetical model that fits the data perfectly, meaning it has as many parameters as data points. In such a model, the predicted values are exactly equal to the observed values, and therefore the log-likelihood reaches its maximum possible value. This model serves as a benchmark or upper bound against which the fit of any other model, including our fitted model, is compared. Hence, the difference \(\ell_i^{S} - \ell_i^{F}\) reflects how much worse our model fits the data compared to this idealised scenario.
While the residuals discussed above can provide useful insights into model fit, they are often less informative in the context of count data, especially for binary data and low counts. In a linear geostatistical model, the distinction between these residuals can be expressed explicitly. Let \(r = y-D\widehat\beta\) and \(V=\Sigma+\tau^2I\). The predictors of the structured spatial effect and the remaining response error are
\[ \widehat W=\Sigma V^{-1}r, \qquad \widehat e=\tau^2V^{-1}r. \]
The response residual, Pearson residual based on the conditional variance, and Gaussian deviance residual are proportional to \(\widehat e\). The random effects residual \(\widehat W\) is a different quantity, and \(\widehat W+\widehat e=r\). It is also easy to show that the correlation structure of the estimated response residuals is related to the spatial correlation matrix of the fitted geostatistical model (see Section 3.5.3). Whether this correlation is unimportant or poses a problem, particularly when residuals are used to assess whether the model has fully accounted for spatial dependence, requires further investigation.
Given these complexities, especially in the context of count data, we have chosen not to delve deeper into residual diagnostics in this edition of the book. Instead, model validation will be addressed more comprehensively in the context of geostatistical prediction in Chapter 4.
Nevertheless, we note that some standard diagnostic tools from regression analysis (such as scatter plots of fitted values versus residuals) can still be informative and are worth exploring as part of routine model assessment (Chapter 9, Weisberg 2014).
3.5 Theory
3.5.1 Maximum likelihood estimation for non-Gaussian mixed models with independent random effects
In this section, we provide more technical details on the derivation of the likelihood function for the mixed model specified in Equation 3.7. While not essential for understanding the application and interpretation of geostatistical models, this section provides useful insights into computational challenges inherent to the models presented in this book. In this section, we restrict our attention to the case when the outcome \(Y_i\) is a count.
Let \(\theta = (\beta, \tau^2)\) denote the vector of unknown parameters, where \(\beta\) is the vector of regression coefficients and \(\tau^2\) is the variance of the random effects \(Z = (Z_1, \ldots, Z_n)\). The likelihood function for a set of observations \(Y = (Y_1, \ldots, Y_n)\), for \(i=1,\ldots,n\), under the model in Equation 3.7, is obtained by deriving the marginal distribution of \(Y_i\), hence
\[ L(\theta) = [Y] = \int_{R^n} \: [Z] [Y \:|\: Z] \: d Z \]
Since the elements \(Z_i\) in \(Z\) are mutually independent, we can write
\[ [Z] = \prod_{i=1}^n [Z_i] \]
and, using the conditional independence among the \(Y\) given \(Z\), we obtain
\[ [Y \:|\: Z] = \prod_{i=1}^n [Y_i | Z_i]. \]
From these, it then follows that
\[ L(\theta) = \int_{R^n} \: \prod_{i=1}^n [Z_i] \prod_{i=1}^n [Y_i | Z_i] \: d Z = \prod_{i=1}^n \int_{R} \: [Z_i] [Y_i \:|\: Z_i] \: d Z_i. \tag{3.26}\]
This shows that under the model in Equation 3.7, whilst accounting for overdispersion, this model does not account for any source of correlation between the observations. This is because the likelihood function can be expressed, as shown in the equation above, as the product of \(n\) likelihoods associated with each of the observations.
Another important observation is that the integrals that express the likelihood in Equation 3.26, cannot be solved analytically and need to be approximated. Several solutions are available that are either based on maximum likelihood estimation or Bayesian methods of inference. The function glmer() used in Section 3.1.3.2 provides both the Laplace and Gaussian quadrature approximations options. The former consists of finding the peak of the integrand - i.e. the function that is being integrated - and approximating that with a Gaussian distribution centred at that peak, using the curvature of the original distribution to determine the width of the Gaussian distribution. The Gaussian quadrature is another numerical integration technique used to approximate the integral of a function. It involves choosing specific points (nodes) and weights within the integration interval to approximate the integral. By strategically selecting these points and weights based on a Gaussian distribution, Gaussian quadrature can provide highly accurate results compared to other numerical integration methods, such as the Laplace approximation. For a comprehensive overview of these and other methods of approximations of intractable integrals, we consider Bishop (2006) to be an accessible textbook on this topic.
3.5.2 Maximum likelihood estimation for non-Gaussian generalized linear geostatistical models
To derive the likelihood function of a generalized linear geostatistical model (GLGM), we consider the special case of a single observation per location without loss of generality. Let \(Y\) and \(Z\) denote the same vectors of random variables introduced in the previous section, with components \(Y_i\) and \(Z_i\) corresponding to the set of locations \(X = (x_1, \ldots, x_n)\). Additionally, we define \(S(X) = (S(x_1), \ldots, S(x_n))\) to represent the spatial Gaussian process at the sampled locations \(X\). In this case the vector of unknown parameters is \(\theta = (\beta, \sigma^2, \phi, \tau^2)\) and, as done in the previous section, we can obtain the likelihood function though the marginal distribution of the data, hence we write
\[ L(\theta) = [Y] = \int_{R^n} \int_{R^n} \: [Z] [S] [Y \:|\: Z, S] \: d Z d S \tag{3.27}\]
To simplify the explanation and notation in this section and the following ones, we will first rewrite the likelihood function in terms of the random effect \(W_i = S(x_i) + Z_i\) and will denote this in its vector form as \(W = (W_1, \ldots, W_n)\). Hence we write \[ L(\theta) = [Y] = \int_{R^n} \: [W] [Y \:|\: W] \: d W \tag{3.28}\] In this model formulation, we assume that \(Y\) conditionally on \(W\) is a set of mutually independent outcomes, hence
\[ [Y | W] = \prod_{i=1}^n [Y_i | W_i]. \]
However, unlike in the previous section, the distribution of \([W]\) is a Multivariate Gaussian distribution with correlated components, hence we cannot split the integral into separate univariate ones as before.
More specifically, \(W\) follows a Multivariate Gaussian distribution with mean vector \(0\) and covariance matrix \(\Omega\), whose \((i, j)\)-th entry is
\[ [\Omega]_{ij} = \begin{cases} \sigma^2 \rho(u_{ij}; \phi, \kappa) & \mbox{if }i\neq j \\ \sigma^2 + \tau^2 & \mbox{if }i = j \end{cases} \]
where \(\rho(u;\phi,\kappa)\) is the Matern correlation function with scale parameter \(\phi\) and smoothness parameter \(\kappa\); see Equation 1.5.
In the next two sections we provide a brief overview of Monte Carlo maximum likelihood and the use of Markov chain Monte Carlo to approximate Equation 3.28.
3.5.2.1 Monte Carlo maximum likelihood
In maximum likelihood estimation, we aim to maximise the likelihood function to find the parameters that make our data the most likely to have been observed. Monte Carlo Maximum Likelihood (MCML) (see the seminal work by Geyer (1991)) is designed to handle situations where the likelihood function is intractable. In the case of a non-Gaussian generalized linear geostatistical model, its likelihood function involves a high-dimensional integral which can be approximated using MCML.
One way to apply Monte Carlo methods in this context is through importance sampling. Importance sampling allows us to approximate the likelihood function by rewriting the likelihood in a form that can be estimated using samples from an easier distribution, known as the importance sampling distribution. To this end, we then rewrite the likelihood ratio as
\[ \begin{aligned} \frac{L(\theta)}{L(\theta_0)} &= \int \frac{[W, Y; \theta]}{[W, Y; \theta_0]} [W \mid Y; \theta_0] \: dW \\ &= E\left\{\frac{[W, Y; \theta]}{[W, Y; \theta_0]}\right\}, \end{aligned} \tag{3.29}\]
where \(\theta_0\) represents an initial guess for \(\theta\). The expectation is taken with respect to \([W \mid Y; \theta_0]\), which is the importance sampling distribution. Because \(L(\theta_0)\) does not depend on \(\theta\), maximising this likelihood ratio is equivalent to maximising \(L(\theta)\).
Let \(W_{(j)}\) represent the \(j\)-th sample out of \(B\) drawn from the distribution \([W \mid Y; \theta_0]\). To simulate samples from \([W \mid Y; \theta_0]\), Markov chain Monte Carlo (MCMC) algorithms are required and we explain more in detail in the next section.
We approximate the likelihood ratio in Equation 3.29 as follows:
\[ R_{MC}(\theta, \theta_0) = \frac{1}{B} \sum_{j=1}^B \frac{[W_{(j)}, Y; \theta]}{[W_{(j)}, Y; \theta_0]}. \tag{3.30}\]
As \(B \rightarrow \infty\), the Monte Carlo approximation \(R_{MC}(\theta, \theta_0)\) converges to the exact likelihood ratio in Equation 3.29. It therefore estimates the likelihood up to the constant \(L(\theta_0)\). The maximisation of this approximation is performed using an optimisation algorithm that uses the analytical expressions of the first and second derivatives. This approach is especially useful for exploring the likelihood surface, which can be relatively flat, particularly in the space of covariance parameters. To improve the approximation of the likelihood ratio in Equation 3.30, this process can be repeated iteratively, replacing \(\theta_0\) after each step with the estimate obtained from equation Equation 3.30, until convergence is achieved.
However, due to Monte Carlo error, the estimates after each iteration may still vary considerably, making convergence difficult to reach if very strict criteria are used. In practice, relaxed convergence criteria or a fixed maximum number of iterations (typically 3-5) are employed with a sufficiently large number of Monte Carlo samples (e.g., \(10,000\)). Under these conditions, we find that the procedure provides reliable parameter estimates, and the impact of Monte Carlo error on spatial prediction has been observed to be negligible.
3.5.2.2 Markov chain Monte Carlo: sampling from the predictive distribution of the spatial Gaussian process
Building on the outline of the MCML method from the previous section, we now address the task of sampling from the distribution \([W \mid Y]\) to generate the samples \(W_{(j)}\) used in equation Equation 3.30. To simplify the notation, we will use \([W \mid Y]\) in place of \([W \mid Y; \theta_0]\), but it is important to emphasise that the parameters \(\theta_0\) are considered fixed when sampling from \([W \mid Y]\). Since direct simulation from \([W \mid Y]\) is not feasible, we must rely on Markov chain Monte Carlo (MCMC) methods, which we will briefly introduce next.
The Metropolis-Hastings algorithm is a foundational technique in the field of MCMC methods, which are widely used for sampling from complex probability distributions. The primary goal of MCMC methods is to generate a sequence of samples from a target distribution where direct sampling is challenging, such as the distribution of \([W \mid Y]\). The Metropolis-Hastings algorithm, first introduced by Metropolis et al. (1953) and later generalized by Hastings (1970), achieves this by constructing a Markov chain that has \([W \mid Y]\) as its stationary distribution. In what follows, we shall use \(\pi(W)\) to denote the density function of \([W \mid Y]\) up to a constant of proportionality.
Given a current realisation \(W_{(k)}\) of the distribution \([W \mid Y]\) in the Markov chain, the algorithm proceeds as follows.
- Propose a new realisation \(W_{(k+1)}\) from \(q(W_{(k+1)} \mid W_{(k)})\), where \(q(\cdot \mid \cdot)\) is the density function of a chosen proposal distribution.
- Compute the acceptance probability:
\[ \alpha = \min \left(1, \frac{\pi(W_{(k+1)}) q(W_{(k)} \mid W_{(k+1)})}{\pi(W_{(k)}) q(W_{(k+1)} \mid W_{(k)})} \right). \] - Accept the proposed realisation \(W_{(k+1)}\) with probability \(\alpha\). Otherwise, retain the current realisation \(W_{(k)}\).
The algorithm ensures that, under certain conditions, the chain converges to the target distribution \(\pi(W)\) as the number of steps increases Roberts and Tweedie (1996).
The functions of the RiskMap package for fitting Binomial and Poisson geostatistical models implement the Metropolis-Adjusted Langevin Algorithm (MALA), which is an improvement over traditional Metropolis-Hastings MCMC methods. Introduced by Rossky et al. (1978) and further developed by Roberts and Tweedie (1996), MALA incorporates gradient information from the target distribution to propose more informed moves, potentially leading to faster convergence and more efficient exploration of the target distribution.
MALA combines the Metropolis-Hastings framework with Langevin dynamics, where proposals \(q(\cdot \mid \cdot)\) are defined as \[ W_{(k+1)} = W_{(k)} + \frac{\epsilon^2}{2} \nabla \log \pi(W_{(k)}) + \epsilon \, \eta, \] where \(\epsilon\) is a step size parameter, \(\nabla \log \pi(W_{(k)})\) is the gradient of the log of the target density with respect to \(W_{(k)}\), and \(\eta \sim \mathcal{N}(0, I)\) is a Gaussian random variable with mean 0 and independent components with unit variance.
By leveraging the gradient, MALA can propose moves that are more likely to be accepted, especially in high-dimensional spaces or regions of low probability density.
In the context of Laplace sampling for geostatistical models, the Metropolis-Adjusted Langevin Algorithm (MALA) is used to construct an efficient proposal distribution. The primary objective is to generate samples from a conditional distribution of latent Gaussian processes, which is challenging due to high dimensionality and complex dependence structures.
The Laplace sampling method, as implemented in the RiskMap package, uses a specific form of MALA. Here, the proposal distribution is derived from a second-order Taylor expansion around the mode of the target distribution \([W \mid Y]\). Let \(\hat{w}\) be the mode of this conditional distribution, and let \(\hat{\Omega}\) be the inverse of the negative Hessian matrix evaluated at \(\hat{w}\), which approximates the local covariance structure of the distribution. We then define the random variable \(\tilde{W}\) as \[ \tilde{W} = \hat{\Omega}^{-1/2} (W - \hat{w}) \]
This linear transformation decorrelates the components of \(W\), making the MCMC updates more efficient (Geyer 1994). RiskMap applies MALA to \(\tilde{W}\) and transforms the resulting samples back to the original scale of \(W\). This approach adapts the proposal to the local geometry of the target distribution, which can improve mixing and reduce the correlation between successive samples (see Christensen (2004)).
Once the MCMC algorithm has completed all iterations, it is important to ensure that it has converged to the target distribution. The diagnostic methods listed below can be used to assess convergence using the obtained MCMC samples.
Trace plots: Plotting the sampled values against the iterations can help visually inspect if the chains have reached a stationary distribution. A well-mixed chain should fluctuate around a stable level without a persistent trend or abrupt change.
Autocorrelation plots: These plots measure the correlation between samples as a function of the lag between them. Low autocorrelation across lags indicates good mixing, suggesting convergence.
-
Effective sample size (ESS): The ESS measures the number of independent samples that an MCMC chain is equivalent to, accounting for the autocorrelation in the chain. It is calculated as:
\[ \text{ESS} = \frac{N}{1 + 2 \sum_{k=1}^{\infty} \rho_k}, \]
where \(N\) is the total number of MCMC samples and \(\rho_k\) is the autocorrelation at lag \(k\). A high ESS indicates that the chain has good mixing properties and the samples are nearly independent. The summation of autocorrelations up to a sufficiently large lag is used to estimate the reduction in effective sample size due to the dependence between samples. A small ESS suggests high autocorrelation and that more samples are needed to achieve the same precision as a set of independent samples.
Other diagnostics can provide complementary information. The Gelman–Rubin diagnostic compares variation between and within multiple independently initialised chains, so it cannot be calculated from a single chain alone (Gelman and Rubin 1992). The Geweke diagnostic compares the means of an early and a late segment of one chain (Geweke 1992). No single diagnostic proves convergence, and several complementary checks should be considered together.
The quality and computational cost of an MCMC run are influenced by three control parameters.
- Number of simulations: Increasing the number of simulations allows the algorithm to explore the posterior distribution more thoroughly, reducing the variability in estimates. A higher number of simulations typically leads to more accurate results but also requires more computational time.
- Burn-in: The burn-in phase refers to the initial set of iterations where the Markov chain has not yet converged to its stationary distribution. These initial samples are discarded to ensure that the remaining samples represent the target distribution more accurately. Properly setting the burn-in period helps avoid bias due to initial conditions.
- Thinning: Thinning retains every n-th sample. It may reduce storage requirements and the autocorrelation between retained samples, but it also discards samples and does not create additional information. The retained samples need not be independent. The number of iterations should therefore be chosen by considering the resulting ESS and the precision required for the analysis.
Initial settings can be based on prior experience and then revised in light of the diagnostics. The next section shows how these settings are specified in RiskMap.
For a more thorough and complete overview of MCMC algorithms, chapters 10 to 13 of Gelman et al. (2013) are recommended reading.
3.5.2.3 Summary with example
In summary, here are the steps required to perform Monte Carlo Maximum Likelihood (MCML):
- Start by selecting an initial guess for the model parameters, denoted as \(\theta_0\).
- Simulate samples from \([W \mid Y; \theta_0]\) using MCMC methods. Perform a check on the MCMC samples to assess convergence of the MCMC algorithm.
- Compute and maximise the approximate likelihood function in Equation 3.30.
- Set the newly estimated parameter as \(\theta_0\) for the next iteration.
- Repeat steps 2 to 4 until convergence is achieved.
- Monitor the estimates for convergence. If strict convergence criteria are set, Monte Carlo error may cause fluctuations, making it harder to achieve convergence.
- Perform additional iterations, if necessary. For inference on model parameters, continue iterating to reduce Monte Carlo error. For spatial prediction, a pragmatic approach may involve stopping after a second iteration with 10,000 samples, as prediction is less affected by Monte Carlo error.
RiskMap automates these steps through glgpm(). The MCMC sample should still be inspected by the user. Repeated MCML updates are configured with set_control_mcml() and passed through control_mcml; its default of one update preserves the usual single-stage fit. With repeated updates, convergence requires both a sufficiently small change in the parameter estimates, measured relative to their standard errors, and adequate overlap of the importance weights. The max_iterations argument sets an upper limit rather than guaranteeing convergence. In Section 3.3.1, estimation is performed using a geostatistical model fitted to prevalence data. The start_pars argument allows initial values for the model parameters to be supplied. If it is omitted, as in Section 3.3.1, the following defaults are used:
- The initial guess for the regression coefficients \(\beta\) is set to the maximum likelihood estimates obtained from a standard generalized linear model (GLM), i.e., a GLM without random effects or overdispersion.
- The initial values for both \(\sigma^2\) (the variance of the spatial Gaussian process) and \(\tau^2\) (the variance of the nugget term, if included) are set to 1.
- The initial guess for the scale parameter of the spatial correlation function \(\phi\) is based on the 0.1 quantile of the distances between all pairs of locations.
For Binomial and Poisson models with their canonical links, RiskMap also fits an equivalent transformed-Gaussian model with the same spatial, nugget and random-effect structure. It compares this candidate with the default values above using a Laplace approximation and retains the better starting point. Values supplied explicitly through start_pars remain authoritative and bypass this automatic comparison.
In our experience, this automatic procedure provides useful starting values across a range of settings. Further MCML updates or longer MCMC samples may still be needed, depending on the parameter stability, importance-weight overlap and MCMC diagnostics.
An alternative approach for setting initial guesses for the covariance parameters \(\sigma^2\), \(\tau^2\), and \(\phi\) is to apply an interpolation technique to the variogram, as proposed in the seminal work by Cressie (1985). However, in our experience, this approach does not consistently provide better approximations of the likelihood function compared to the method mentioned above. This may be due to the instability in estimating the variogram based on residuals from mixed models, as well as the sensitivity of empirical variogram estimation to the chosen distance breaks.
To better illustrate some of the points above, we continue the analysis of the Liberia data from Section 3.3.1. Diagnostic summaries on the convergence of the MCMC algorithm can be obtained with the plot_mcmc() function by passing the fitted model fit_liberia to the function. The fitted object must contain retained MCMC samples, as requested earlier with return_samples = TRUE. These diagnostic plots are generated by RiskMap using base graphics. They are retained here as package diagnostic output rather than as finished explanatory figures.
plot_mcmc(fit_liberia)
The plots shown in Figure 3.12 are for the MCMC samples of the averaged spatial random effect, that is \(\bar{S} = \sum_{i=1}^n S(x_i) /n\). Inspecting \(\bar{S}\) avoids checking every individual component \(S(x_i)\). A problem in a small number of components may be obscured by the average, so individual components should also be inspected when there is a reason for concern. This can be done as follows.
plot_mcmc(fit_liberia, check_mean = FALSE, component = 10)
Both Figure 3.12 and Figure 3.13 show reasonably good mixing and relatively low autocorrelation. The MCMC settings for glgpm() are created by set_control_mcmc() and passed through the control_mcmc argument. By default, set_control_mcmc() uses 12,000 iterations, a burn-in of 2,000 and thinning of 10, giving 1,000 retained samples. This configuration may be sufficient for the first MCML run. More iterations may be needed to obtain a larger ESS and improve the precision of the Monte Carlo approximation.
# Set starting values for the model with a nugget effect
start_liberia <- coef(fit_liberia)
start_liberia$tau2 <- 0.1
# Refit the model and increase the MCMC samples
fit_liberia2 <- glgpm(npos ~ log(elevation) + gp(nugget = TRUE),
denominator = ntest, data = liberia,
model_crs = 32629,
start_pars = start_liberia,
control_mcmc = set_control_mcmc(n_sim = 110000,
burnin = 10000,
thin = 10,
seed = 2026),
control_mcml = set_control_mcml(max_iterations = 2),
return_samples = TRUE,
family = "binomial", messages = FALSE)plot_mcmc(fit_liberia2)
The code above allows up to two MCML updates, regenerating the importance sample around the estimate from the preceding update. Each update uses 10,000 retained MCMC samples. As shown in Figure 3.14, these MCMC samples exhibit good mixing and low autocorrelation, which indicates that they adequately approximate the likelihood function of the Binomial geostatistical model. Although additional MCMC iterations could be run to obtain estimates closer to the true maximum likelihood estimate (since we are approximating the likelihood function and not obtaining the exact MLE), we choose not to do so here. This decision is due to both the computational time required and the minimal impact of Monte Carlo error on spatial predictions in this case.
3.5.2.4 Alternative approaches
In the sections above, we focused specifically on the implementation of the Monte Carlo Maximum Likelihood (MCML) method in RiskMap. However, other approaches are available, some of which are computationally more efficient than the MCML method we have illustrated. Our preference for this estimation method is based on two main reasons: 1) it delivers robust performance across a wide range of data scenarios, and 2) it does not require the specification of arbitrary prior distributions when prior information on the model parameters is unavailable.
Alternative methods worth mentioning include the expectation-maximisation algorithm for maximum likelihood estimation (Zhang 2002), integrated nested Laplace approximation (INLA) (Rue et al. 2009), and implementations of the Metropolis-Adjusted Langevin Algorithm (MALA) for Bayesian inference (Christensen et al. 2006). While these methods are valuable, explaining and comparing them is not our focus here, as it would detract from the practical orientation of this book. The RiskMap package itself will continue to evolve as computational improvements are made, and future editions of this book will update these sections without affecting the applied analysis content.
3.5.3 The Hat matrix in linear geostatistical models
We consider the geostatistical model \[ Y_i = d(x_i)^\top \beta + S(x_i) + Z_i \] where \(Y_i\) is the observed response at location \(x_i\), \(d(x_i)\) is a \(p\)-dimensional vector of covariates, \(\beta\) is a \(p \times 1\) vector of regression coefficients, \(S(x_i)\) is a spatially structured random effect, and \(Z_i \sim N(0, \tau^2)\) is independent Gaussian noise. The model can be rewritten as \[ Y = D\beta + S + Z \] where \(Y = (Y_1, \ldots, Y_n)^\top\) is the \(n \times 1\) vector of observations, \(S = (S(x_1), \ldots, S(x_n))^\top\) is the vector of spatially dependent random effects, \(Z = (Z_1, \ldots, Z_n)\) is the noise, and \(D\) is an \(n \times p\) matrix of covariates given by
\[ D = (d_1(X), \ldots, d_p(X)) \]
with each column \(d_j(X) = (d_j(x_1),\ldots,d_j(x_n))^\top\) representing the values of the \(j\)-th covariate across all locations.
This formulation accounts for both fixed effects through \(D\beta\) and spatial correlation through \(S(x_i)\). The covariance matrix of \(Y\) is \[ V = \Sigma + \tau^2 I \] where \(\Sigma\) is the spatial covariance matrix of \(S\). Assuming the covariance parameters are known, the generalized least squares estimator for \(\beta\) is
\[ \hat{\beta} = (D^\top V^{-1} D)^{-1} D^\top V^{-1} Y. \]
Treating the covariance parameters as known, or fixed at their estimates, the empirical best linear unbiased predictor of \(S\) is
\[ \hat{S} = \Sigma V^{-1} (Y - D\hat{\beta}). \]
The fitted values can therefore be written as
\[ \hat{Y} = D\hat{\beta} + \hat{S}. \]
Writing \(\widehat Y=HY\) gives the smoother matrix
\[ H = \Sigma V^{-1} + \tau^2 V^{-1}D(D^\top V^{-1}D)^{-1}D^\top V^{-1}, \]
where we have used \(I-\Sigma V^{-1}=\tau^2V^{-1}\). The matrix \(H\) is symmetric because \(\Sigma\) and \(V\) commute. Unlike the hat matrix in ordinary least squares, however, it is not idempotent and is therefore a shrinkage smoother rather than a projection. Under the usual positive-definite covariance assumptions, \(H\) has full rank. It also satisfies \(HD=D\), and its effective degrees of freedom \(\operatorname{tr}(H)\) lie between \(p\) and \(n\). The diagonal element \(h_{ii}\) describes the influence of \(Y_i\) on its own fitted value.
Under these assumptions, the difference between an observation and its fitted value is the empirical best linear unbiased predictor of the nugget effect:
\[ \widehat Z = Y - \widehat Y = \tau^2V^{-1}(Y-D\widehat\beta) = (I-H)Y. \]
This does not mean that \(\widehat Z\) is the unobserved realisation of \(Z\). Rather, it is the model-based predictor obtained after estimating \(\beta\) by generalized least squares. Defining
\[ Q=V^{-1}-V^{-1}D(D^\top V^{-1}D)^{-1}D^\top V^{-1}, \]
we have \(I-H=\tau^2Q\) and hence
\[ \operatorname{Var}(\widehat Z)=\tau^4Q. \]
The estimated nugget residuals are therefore correlated, with a covariance structure determined by both the fitted spatial covariance \(\Sigma\) and the fixed-effects design matrix \(D\).
3.6 Review questions
- Describe how you would perform exploratory analysis between counts data and covariates, distinguishing between unbounded counts (e.g. Poisson), bounded counts (e.g. Binomial) and binary outcomes.
- Define the theoretical and empirical variogram, making appropriate use of equations.
- How do we use the empirical variogram to test for spatial correlation? What are the limitations of this approach?
- Define the class of generalized linear geostatistical models and list all the assumptions.
- Explain intuitively what are the computational challenges related to the fit of geostatistical models and what approaches can be used for parameter estimation.
3.7 Exercises
Consider the Binomial mixed model with linear predictor as defined in Equation 3.7. By editing the code for the simulation shown in Section 3.1.2, generate a graph as in Figure 3.6 under the two following scenarios: \(i\)) \(\tau^2 = 0.2\) and \(n_i = 100\); \(ii\)) \(\tau^2 = 0.1\) and \(n_i = 1\). How does the variance of \(Y_i\) change under \(i\)) and \(ii\)) in comparison to Figure 3.6? How do you explain the differences?
-
Similarly to the previous exercise, consider a Poisson mixed model with linear predictor
\[ \log\left\{\mu_i\right\} = \beta_0 + Z_i, \]
where \(Z_i\) are a set of mutually independent Gaussian variables with mean 0 and variance \(\tau^2\). Using the code shown in Section 3.1.2, carry out a simulation study to compute the variance of \(Y_i\) and generate a graph similar to Figure 3.6 to compare the variance of the Poisson mixed model with that of a standard Poisson model. Generate the graph for different values of \(\tau^2\) and summarize your findings. NOTE: In this simulation the offset \(n_i\) can be set to 1.
Create an R function that computes the cloud variogram. As explained in Section 3.1.3, the cloud variogram is obtained by plotting \(\hat{V}_{ij}\) (see Equation 3.11) against the distances \(u_{ij}\). The function should take as input a dataset with three columns: the variable for which the cloud variogram is to be computed; and two columns corresponding to the location of the data. Then, use this function to create the cloud variogram for the model for river blindness in Equation 3.13. How does this compare to the empirical variogram that takes the averages within predefined distance classes, as shown in Figure 3.7?
-
Fit a Binomial mixed model to the Liberia dataset on river blindness without any covariates, i.e.
\[ \log\left\{\frac{p(x_i)}{1-p(x_i)}\right\} = \beta_0 + Z_i. \]
Making use of the R code presented in Section 3.1.3.1, use the function
variogram()to generate the empirical variogram for this model and compare this to the empirical variogram of Figure 3.7. What differences do you observe? -
Fit a Poisson mixed model to the anopheles mosquito data using elevation as a covariate, i.e.
\[ \log\{\mu(x_i)\} = \beta_{0} + \beta_{1} d(x_i) + Z_i \]
where \(d(x_i)\) is the measured elevation at location \(x_i\). How strong is the overdispersion in the data? After fitting the model, extract the estimates of the random effects \(Z_i\) and compute the empirical variogram with the 95% envelope for spatial independence. Repeat this for different classes of distances by changing the input passed to
breaksin thevariogram()function. How do the different specifications ofbreaksaffect the results? Parametric bootstrap. Using the model fitted to the Anopheles gambiae mosquitoes in Section 3.3.2, simulate 1000 datasets using the
simulate_glgpm()function and fit the model to each simulated dataset. Use the resulting 1000 estimates to compute 95\(\%\) confidence intervals for each of the model parameters and compare this with the confidence intervals obtained from the model summary.-
In Section Section 3.2.1, we explained that for a linear model with one observation per location, we cannot simultaneously estimate both the measurement error variance \(\omega^2\) and the nugget effect variance \(\tau^2\). However, if \(\omega^2\) is known, we can fix this value and estimate the nugget effect. Suppose that for the Galicia data, the measurement device has a known standard error of \(\omega = 0.15\) on the log-scale (i.e., \(\omega^2 = 0.0225\)). Fit the model
\[ Y_i = \mu + S(x_i) + Z_i + U_i \]
using
glgpm()by settingnugget = TRUEandfix_var_me = 0.0225in thegp()function, with \(\kappa = 1.5\). Compare the estimates of \(\sigma^2\) and \(\phi\) with those obtained from the model without the nugget effect (Equation 3.16). What is the estimated nugget variance \(\hat{\tau}^2\)? Repeat this analysis for different values of fixed measurement error variance: \(\omega^2 \in \{0.01, 0.05, 0.10\}\). How do the estimates of \(\sigma^2\), \(\phi\), and \(\tau^2\) change as the assumed measurement error variance increases?