Impact of jittering on raster- and distance-based geostatistical analyses of DHS data

Abstract

Fine-scale covariate rasters are routinely used in geostatistical models for mapping demographic and health indicators based on household surveys from the Demographic and Health Surveys (DHS) program. However, the geostatistical analyses ignore the fact that GPS coordinates in DHS surveys are jittered for privacy purposes. We demonstrate the need to account for this jittering, and we propose a computationally efficient approach that can be routinely applied. We use the new method to analyse the prevalence of completion of secondary education for 20-49 year old women in Nigeria in 2018 based on the 2018 DHS survey. The analysis demonstrates substantial changes in the estimates of spatial range and fixed effects compared to when we ignore jittering. Through a simulation study that mimics the dataset, we demonstrate that accounting for jittering reduces attenuation in the estimated coefficients for covariates and improves predictions. The results also show that the common approach of averaging covariate values in windows around the observed locations does not lead to the same improvements as accounting for jittering.

1 Introduction

Fine-scale spatial estimation of demographic and health indicators has become commonplace (Burstein et al., 2019; Utazi et al., 2019; Local Burden of Disease Vaccine Coverage Collaborators, 2021). This article is focused on prevalences, which include many important indicators such as completion of secondary education, neonatal mortality, and vaccination coverage (Fuglstad et al., 2021). For low- and middle-income countries (LMICs), the household surveys conducted by the Demographic and Health Surveys (DHS) Program are a crucial data source. Geographic information in DHS data is given through GPS coordinates, which describe centres of clusters of households. However, cluster centres are randomly displaced by up to 10 km before being published in order to protect participants' privacy (Burgert et al., 2013). We refer to such small random displacements as jittering of the GPS coordinates.

In global health, it is common practice to ignore jittering and estimate risk using a standard geostatistical model with a binomial likelihood. The latent spatial variation in risk is modelled as the combination of raster- and distance-based covariates and a Gaussian random field (GRF). However, covariate values extracted from rasters can vary widely on the distance scale of jittering. Using the covariate value at the jittered location instead of the original location induces a non-standard form of measurement error (Gustafson, 2003). This may in turn lead to attenuation of effect estimates and errors in uncertainty. Furthermore, not accounting for the positional uncertainty for the GRF, artificially reduces estimated spatial dependency and may reduce predictive power as well (Cressie and Kornak, 2003; Fanshawe and Diggle, 2011; Fronterrè et al., 2018).

To address uncertainty in covariates, Perez-Heydrich et al. (2013, 2016) suggested (a) to use regression calibration in the context of distance-based covariates (Warren et al., 2016), and (b) to average spatial covariates within a 5 km buffer zone for continuous and categorical rasters. However, this approach does not address the issue of attenuation of associations. Fanshawe and Diggle (2011) proposed a Bayesian approach to account for positional uncertainty for the GRF, but did not propagate uncertainty in the covariates, and only used Gaussian likelihoods that are not applicable to prevalences. The approach was also computationally expensive, but Fronterrè et al. (2018) made the approach computationally efficient and demonstrated its applicability to analyse malnutrition based on DHS data.

Recently, Wilson and Wakefield (2021) formulated a full geostatistical model for DHS data that includes an observation model for the jittered GPS coordinates, and estimated the model with integrated nested Laplace approximations (INLA) (Rue et al., 2009) within Markov chain Monte Carlo (MCMC) (Gómez-Rubio and Rue, 2018). Their approach addresses the effect of positional uncertainty on both the spatial covariates and the GRF, but was computationally expensive with 1000 MCMC iterations requiring 52 hours in their simulation study. Altay et al. (2022) proposed a similar model as Wilson and Wakefield (2021), but used a more efficient inference scheme with computation time being measured in minutes instead of hours. Their approach was made possible through an approximation of the likelihood, the stochastic partial differential equations (SPDE) approach (Lindgren et al., 2011), and Laplace approximations through template model builder (TMB) (Kristensen et al., 2016).

The simulation study in Altay et al. (2022) revealed that small spatial ranges for the GRF or larger jittering than the DHS scheme were required to see substantial improvements with the new approach over ignoring jittering. Altay et al. (2022) focused on the impact of jittering on the GRF, and lacked raster- and distance-based covariates. Such covariates are far more variable at small spatial scales than a smoothly varying GRF. The aim of this article is to extend the approach in Altay et al. (2022) to a full generalized geostatistical model for prevalence, and to demonstrate that ignoring jittering can lead to attenuation of associations and reduced predictive power when analysing DHS data. We show this via a spatial analysis of the prevalence of secondary education completion among women aged 20-49 in 2018 based on the 2018 Nigeria DHS (NDHS2018) (National Population Commission - NPC and ICF, 2019).

In addition to the new approach, which adjusts for jittering, and the standard approach, which ignores jittering, we consider the common approach of averaging covariates in 5 km × 5 km windows around the provided GPS coordinates (Perez-Heydrich et al., 2013, 2016). The three methods cannot be compared with cross-validation since the true coordinates of the clusters are not known. Therefore, we construct a simulation study that mimics the NDHS2018 dataset to compare them in terms of their ability to estimate parameters and to predict risk at unobserved locations. We use bias and root mean square error (RMSE) to assess parameter estimation, and RMSE and continuous rank probability score (CRPS) (Gneiting and Raftery, 2007) to assess predictive ability.

We introduce the datasets and variables of interest in Section 2. Then, in Section 3, we describe the new approach that adjusts for jittering, and discuss its implementation. In Section 4, we evaluate parameter estimation and prediction with the different methods through a simulation study that mimics the prevalence of secondary education completion. Then, in Section 5, we demonstrate the differences between adjusting and not adjusting for jittering when analysing the prevalence of secondary education completion for women aged 20-49 in Nigeria. The article ends with discussion and conclusions in Section 6. The code used in the article can be found in the GitHub repository https://github.com/umut-altay/GeoAdjust.

2 Data sources and variables of interest

Our outcome of interest is completion of secondary education, which is an indicator of social wellbeing and life outcome (Lewin, 2008). Rates vary strongly between women and men, but also between urban and rural areas. According to UNESCO (2019), only 1 % of the poorest girls in low income countries will complete secondary education. If a girl completes secondary education, the risk of HIV infection is reduced by about 50 % (UNAIDS, 2022).

We consider the prevalence of secondary education completion for women aged 20-49 years in Nigeria in 2018. The lower bound was chosen because younger women may not have completed secondary education yet, and the upper limit was chosen since older women are not available in the DHS surveys. The year 2018 was chosen since this corresponds to the most recent DHS survey in Nigeria, NDHS2018.

Nigeria is an LMIC with a population of more than 200 million, and NDHS2018 has data collected from 1389 clusters with responses from 33,398 women aged 20-49, where 15,621 of them reported that they had completed their secondary education. In this article, we use the 1380 clusters with valid GPS coordinates, which have responses from 33,193 women aged 20-49 where 15,490 reported that they had completed their secondary education.

Figure 1(b) shows the direct estimates (Rao and Molina, 2015) computed based on the data and the survey design for the 36 states and one federal capital territory. The corresponding uncertainty is expressed through the coefficient of variation (CV). These 37 areas are the first administrative level (admin1), and the results show considerable variation between areas. The CVs increase when direct estimates are calculated at finer spatial scales due to observations needing to be distributed amongst more areas, resulting in smaller sample sizes per area. The second administrative level (admin2) for Nigeria, for example, consists of the 774 local government areas shown in Figure 1(a). The red dots in the figure indicate the 1380 clusters with GPS coordinates within Nigeria available in NHDS2018. The national boundary, admin1 boundaries and admin2 boundaries are based on GADM version 4.0 (GADM).

Figure 1

Maps of (a) admin2 areas and cluster locations for NDHS2018, and (b) admin1 areas with direct estimates.

We expect the prevalence of completion of secondary education to be closely related to the access to educational resources, such as technological infrastructure, schools and teachers. We consider five spatial covariates: population count (PopD) (World Pop, 2022), travel time to nearest city (CityA) (Weiss et al., 2018), elevation (Elev) (National Oceanic and Atmospheric Administration, 2022), distance to nearest river or lake (DistW) (Natural Earth, 2012), and urbanicity ratio (Pesaresi et al., 2016). For UrbR, we scale the original covariate to be on a zero to one scale, and for the other four covariates, we use a log(1 + x)-transformation and then center and standardize the covariate rasters across the pixels. The information about the covariate rasters and figures is summarized in Table 1. The covariate rasters are available at different resolutions, and have not been resampled or aggregated to the same resolution since this would involve an extra preprocessing step that would induce misalignment error and potentially add ecological bias (Greenland and Morgenstern, 1989).

Table 1

Summary of covariate rasters providing name, description and figure. CityA, Elev, DistW and PopD are transformed, while UrbR is not.

Name	Description	Figure
PopD	Population count (250 m × 250 m)	2(a)
CityA	Travel time in minutes (1 km × km)	2(b)
Elev	Elevation in meters (1 km × 1 km)	2(c)
DistW	Distance to nearest river or lake in degrees (1 km × 1 km)	2(d)
UrbR	Urbanicity ratio (250 m × 250 m)	2(e)

We aim to map spatial variation at a continuous spatial scale through a geostatistical model. A key challenge with using the above covariates is that the locations shown in Figure 1(a) are not the true locations. The true locations have been randomly displaced while respecting the admin 2 borders. This means that we can only imprecisely extract covariates from the rasters shown in Figure 2. The new method to account for the uncertainty in locations is described in Section 3.

Figure 2

Transformed covariate rasters for Nigeria.

3 Adjusting for jittering in a geostatistical model

3.1 Notation for DHS data

For a given country and DHS household survey, $C$ clusters are visited. These clusters constitute small geographic areas and are collections of households. A total of $n_c$ people at risk are observed and $y_c\le n_c$ individuals have positive outcomes for clusters $c=1,2,\dots ,C$ . The reported GPS coordinates of the cluster centres are $s_c\in {ℝ}^2,c=1,\dots ,C$ . These locations are not the true GPS coordinates, but the jittered GPS coordinates. Additionally, the urban/rural designation is known for each visited cluster.

3.2 Geostatistical model

3.2.1 Model for spatial variation in risk

We envision a spatially varying risk, $r\left(\cdot \right)$ , for the country of interest $\mathcal{D}\subset {ℝ}^2$ modelled through

$$r\left(s\right)={\text{logit}}^{-1}\left(\eta \left(s\right)\right)={\text{logit}}^{-1}\left(x{(s)}^{\text{T}}\beta +u\left(s\right)\right),s\in \mathcal{D},$$

where $x\left(\cdot \right)$ is a $p$ -dimensional vector of covariates, $\beta$ is a $p$ -dimensional vector of covariate effect sizes, and $u\left(\cdot \right)$ is a Matérn GRF. The Matérn covariance function with smoothness $v=1$ is parametrized as

$$C_{\text{M}}\left(s_1,s_2;{\sigma }_{\text{S}}^2,{\rho }_{\text{S}}\right)={\sigma }_{\text{S}}^2\left(\sqrt{8}\frac{s_1-s_2}{{\rho }_{\text{s}}}\right){\text{K}}_1\left(\sqrt{8}\frac{s_1-s_2}{{\rho }_{\text{s}}}\right),s_1,s_2\in \mathcal{D},$$

where ${\sigma }_{\text{S}}^2$ is the marginal variance, ${\rho }_{\text{S}}$ is the spatial range, and ${\text{K}}_1$ is the modified Bessel function of the second kind, order 1.

3.2.2 Unadjusted model

When jittering is ignored, the reported cluster locations $s_1,\dots ,s_C$ are treated as the true locations. This gives the unadjusted observation model:

$$\begin{array}{l}{y}_c\mid r_c,n_c \sim \text{Binomial}\left(n_c,r_c\right),\\ r_c =r\left(s_c\right)={\text{logit}}^{-1}\left(\eta \left(s_c\right)\right),\end{array}$$

(3.1)

where $r_c$ is the risk in cluster $c$ , for $c=1,\dots ,C$ .

3.2.3 Adjusted model

Let $s_1^{\text{*}},\dots ,s_C^{\text{*}}\in \mathcal{D}$ denote the true locations corresponding to the jittered locations $s_1,\dots ,s_C$ . The adjusted observation model is,

$$\begin{array}{c}{y}_c\mid r_c,n_c \sim \text{Binomial}\left(n_c,r_c\right),s_c\mid s_c^{\text{*}}\sim {\pi }_{\text{Urb}\left[c\right]}\left(s_c\mid s_c^{\text{*}}\right),\\ r_c\mid s_c^{\text{*}} =r\left(s_c^{\text{*}}\right)={\text{logit}}^{-1}\left(\eta \left(s_c^{\text{*}}\right)\right),\end{array}$$

(3.2)

where $r_c$ is the risk in cluster $c$ , and $\text{Urb}\left[c\right]\in \left\{\text{U},\text{R}\right\}$ corresponds to the cluster's urban (U) or rural (R) designation, for $c=1,\dots ,C$ . In this observation model, both $y_c$ and $s_c$ are treated as observed quantities. The unobserved true locations $s_c^{\text{*}}$ are treated as random quantities and assigned a uniform prior $s_c^{\text{*}}\sim \mathcal{U}\left(A\left(s_c\right)\right.$ ), where $A\left(s\right)\in \left\{1,\dots ,K\right\}$ denotes the administrative region containing location $s\in \mathcal{D}$ that a cluster location must be jittered within and $K$ is the number of such administrative regions. For example, the NDHS2018 jittered cluster locations must be in the same admin2 area as the associated true cluster locations. This implies that we treat all true cluster locations $s_c^{\text{*}}$ in the corresponding admin 2 area and within the maximum jittering distance of $s_c^{\text{*}}$ as equally likely a priori.

The jittering distributions ${\pi }_U$ and ${\pi }_R$ follow from the (known) DHS jittering scheme. Then for an urban cluster $c$ , which can be jittered up to $2\text{ km}$ , the jittering distribution is

$${\pi }_{\text{U}}\left(s_c\mid s_c^{\text{*}}\right)\propto \frac{\mathbb{I}\left(A\left(s_c\right)=A\left(s_c^{\text{*}}\right)\right)\cdot \mathbb{I}\left(d\left(s_c,s_c^{\text{*}}\right)<2\right)}{d\left(s_c,s_c^{\text{*}}\right)},s_c\in \mathcal{D},$$

where $d\left(s_c,s_c^{\text{*}}\right)$ is the distance in kilometers between $s_c$ and $s_c^{\text{*}}$ , and $\mathbb{I}$ is the indicator function. For a rural cluster $c$ , which can be jittered up to $5\text{ km}$ except for $1\text{\%}$ of clusters jittered up to $10\text{ km}$ , the jittering distribution is:

$${\pi }_{\text{R}}\left(s_c\mid s_c^{\text{*}}\right)\propto \frac{\mathbb{I}\left(A\left(s_c\right)=A\left(s_c^{\text{*}}\right)\right)}{d\left(s_c,s_c^{\text{*}}\right)}\left[\frac{99\mathbb{I}\left(d\left(s_c,s_c^{\text{*}}\right)<5\right)}{100}+\frac{\mathbb{I}\left(d\left(s_c,s_c^{\text{*}}\right)<10\right)}{100}\right],s_c\in \mathcal{D}.$$

3.2.4 Priors

We assume linear covariate associations, and use the prior $\beta \sim {\mathcal{N}}_p\left(0,25I_p\right)$ . The range, ${\rho }_{\text{S}}$ , and marginal variance, ${\sigma }_{\text{S}}^2$ , of the Matérn GRF is assigned a penalized complexity (PC) prior (Fuglstad et al., 2019). This requires selecting two hyperparameters: the a priori median of range $R_0$ , and the a priori 95th percentile of marginal standard deviation $S_0$ .

3.3 Implementation

3.3.1 Inference scheme

The observation model in Equation (3.2) can be written as,

$$\begin{array}{c}\pi \left(y_c,s_c\mid \eta \left(\cdot \right)\right) ={{\displaystyle \int }}_{{ℝ}^2}\pi \left(y_c,s_c\mid \eta \left(\cdot \right),s_c^{\text{*}}\right)\pi \left(s_c^{\text{*}}\right)\text{d}{s}_c^{\text{*}}\\ ={{\displaystyle \int }}_{{ℝ}^2}\pi \left(y_c\mid \eta \left(s_c^{\text{*}}\right)\right){\pi }_{\text{Urb}\left[c\right]}\left(s_c\mid s_c^{\text{*}}\right)\pi \left(s_c^{\text{*}}\right)\text{d}{s}_c^{\text{*}},\end{array}$$

(3.3)

for $c=1,\dots ,C$ . Let $\theta =\left(\text{log}\left({\sigma }_{\text{S}}^2\right),\text{log}\left({\rho }_{\text{S}}\right)\right)$ . We propose an empirical Bayes approach:

• Step 1: Calculate the maximum a posteriori (MAP) estimate, $\hat{\theta }$ , of $\theta$ using $\pi \left(\theta \mid y_1,\dots ,y_C,s_1,\dots ,s_C\right)$ .

• Step 2: Extract inference about $\beta$ from $\pi \left(\beta \mid y_1,\dots ,y_C,s_1,\dots ,s_C,\theta =\hat{\theta }\right)$ .

• Step 3: Estimate risk $r\left(s\right)$ at location $s$ using

  $$\pi \left(r\left(s\right)\mid y_1,\dots ,y_C,s_1,\dots ,s_C,\theta =\hat{\theta }\right).$$

Two key components are combined for rapid inference: the SPDE approach to approximate the Matérn GRF (Lindgren et al., 2011), and Template Model Builder (TMB) for empirical Bayesian inference (Kristensen et al., 2016).

3.3.2 SPDE approach

For each cluster $c$ , the true location $s_c^{\text{*}}$ is not known, and the observation model in Equation (3.3) involves the spatial field $u\left(\cdot \right)$ at all locations that are compatible with the jittered location $s_c$ . If we replace the integral in Equation (3.3) by a integration scheme using $N_{\text{Int }}$ integration points, we need to evaluate the spatial field at $C\cdot N_{\text{Int }}$ locations. A standard implementation of the Matérn model would result in a dense $C\cdot N_{\text{Int }}\times C\cdot N_{\text{Int }}$ matrix and make computations infeasible even for a few locations.

The SPDE approach (Lindgren et al., 2011) overcomes this issue by approximating the Matérn GRF that results in a sparse precision matrix. First, the area of interest is triangulated with a triangulation consisting of $m$ nodes. Then the GRF $u\left(\cdot \right)$ is approximated by

$$\tilde{u}\left(s\right)={\displaystyle \sum }_{i=1}^m{w}_i{\varphi }_i\left(s\right)$$

(3.4)

where ${\varphi }_i\left(\cdot \right)$ are pyramidal basis functions and $w={\left(w_1\dots w_m\right)}^{\text{T}}$ are weights for the basis functions. The SPDE approach results in a distribution $w\sim {\mathcal{N}}_m\left(0,Q{(\theta )}^{-1}\right)$ , where $Q\left(\theta \right)$ is sparse.

From Equation (3.4), the value at any location is a linear transformation $\tilde{u}\left(s\right)=a{(s)}^{\text{T}}w,s\in \mathcal{D}$ , where $a\left(s\right)\in {ℝ}^m$ is sparse with at most three nonzero elements depending on the location $s$ . This means that the spatial field can be evaluated at a large number of locations quickly.

The SPDE is given by

$$\left({\kappa }^2-\text{Δ}\right)\left(\tau u\left(s\right)\right)=\mathcal{W}\left(s\right),s\in \tilde{\mathcal{D}},$$

where $\kappa >0$ and $\tau >0$ are related to marginal variance and range, $\text{Δ}$ is the Laplacian, $\mathcal{W}\left(\cdot \right)$ is standard Gaussian white noise, and $\tilde{\mathcal{D}}\supset \mathcal{D}$ is an extended domain to reduce boundary effects. We use Neumann boundary conditions to make the problem well defined, and following Lindgren et al. (2011), the effective range and marginal variance are calculated from the SPDE parameters as

$${\rho }_{\text{S}}=\frac{\sqrt{8}}{\kappa }\text{ and }{\sigma }_{\text{S}}^2=\frac{1}{4\pi {\tau }^2{\kappa }^2}.$$

3.3.3 Template model builder

We implement the empirical Bayesian inference scheme by employing the built-in auto-differentiation and Laplace approximations of the TMB R package. Unlike sampling based MCMC methods, TMB uses numerical integration, through Laplace approximations, to perform inference. The auto-differentiation is used to speed up the Laplace approximations. See Kristensen et al. (2016) for the details.

In TMB, one can compute arbitrary likelihoods, and we can approximate the likelihood in Equation (3.3) through the integration scheme

$$\pi \left(y_c,s_c\mid \eta \left(\cdot \right)\right)\propto {\displaystyle \sum }_{k=1}^K{\alpha }_k\pi \left(y_c\mid \eta \left(s_{c,k}^{\text{*}}\right)\right){\pi }_{\text{Urb}\left[c\right]}\left(s_c\mid s_{c,k}^{\text{*}}\right)\pi \left(s_{c,k}^{\text{*}},\right)$$

(3.5)

where ${\alpha }_1,\dots ,{\alpha }_K$ are integration weights. More details are available in Altay et al. (2022). Critically, the integration scheme in Equation (3.5) involves

$$\eta \left(s_{c,k}^{\text{*}}\right)=x{\left(s_{c,k}^{\text{*}}\right)}^{\text{T}}\beta +u\left(s_{c,k}^{\text{*}}\right),k=1,\dots ,K,c=1,\dots ,C.$$

Based on the known jittering distribution, we construct the integration scheme with rings of integration points around each cluster center. The observed cluster center is the first integration point, and we use 5 and 10 rings for the clusters that are located within urban and rural administrative areas, respectively. Each ring consists of 15 angularly equidistant integration points.

Throughout this article we consider the three approaches shown in Table 2. UnAdj denotes the traditional model that does not adjust for jittering, Smoothed denotes a model where the covariates have been averaged over $5\text{ km}\times 5\text{ km}$ windows around each location (Perez-Heydrich et al., 2013, 2016) and FullAdj denotes the new geostatistical model that fully adjust for the positional uncertainty.

Table 2

The three approaches considered in the article.

Approach	Description
UnAdj	Traditional geostatistical model.
Smoothed	UnAdj with covariates averaged over 5 km × 5 km windows.
FullAdj	Geostatistical model adjusting for jittering.

4 Simulation study

In this section, we evaluate the three models illustrated in Table 2 using a known data generating model, where the parameters correspond to FullAdj estimated based on completion of secondary education in Section 5, see Table 4. The data generating model is chosen to achieve realistic scenarios, and the aim is to evaluate accuracy in parameter estimation, and accuracy in predictive distributions. The key interest is how these accuracies vary across different strengths of the signal from the covariates.

We assume that the true spatial risk varies as

$$r\left(s\right)={\text{logit}}^{-1}\left(x{(s)}^{\text{T}}\beta +u\left(s\right)\right),s\in \mathcal{D},$$

(4.1)

where $\mathcal{D}$ is Nigeria, $x\left(\cdot \right)$ is a 6-dimensional spatially varying vector with 1 as the first element denoting the intercept, and the covariates DistW, CityA, Elev, PopD, and UrbR as the five last elements, $\beta ={\left(\mu ,{\beta }_{\text{DistW }},{\beta }_{\text{City }},{\beta }_{\text{Elev }},{\beta }_{\text{PopD }},{\beta }_{\text{UrbR }}\right)}^{\text{T}}$ , and $u\left(\cdot \right)$ is a Matérn GRF. Spatial range ${\rho }_{\text{S}}$ and the marginal variance ${\sigma }_{\text{S}}^2$ of the GRF are set to the values given in Table 4 for FullAdj. Then we construct three scenarios according to weaker, the same, and stronger association to the covariates compared to estimated coefficients under FullAdj in Table 4:

1. SignalLow: $\beta$ is set to 0.5 times the estimated values.

2. SignalMed: $\beta$ is set to 1.0 times the estimated values.

3. SignalHigh: $\beta$ is set to 1.5 times the estimated values.

For each of the three scenarios, we generate $n_{\text{sim }}=50$ true risk surfaces, $r\left(\cdot \right)$ , according to Equation (4.1).

For observations, we follow the design of the NDHS2018 survey. We fix the number of clusters to $C=1,380$ , and for each cluster $c$ , and fix 568 urban locations and 812 rural locations. We also fix the number-at-risk, $n_c$ , according to NDHS2018 for $c=1,\dots ,C$ . For each of the 150 true risk surfaces, the true locations $s_1^{\text{*}},\dots ,s_C^{\text{*}}\in \mathcal{D}$ are simulated as a Poisson point process with intensity proportional to population density. This ensures higher likelihood to place observations at locations with more people, and no locations at locations with no people. Then, for each cluster $c=1,\dots ,C$ , we simulate response $y_c\mid r\left(s_c^{\text{*}}\right),n_c\sim \text{Binomial}\left(n_c,r\left(s_c^{\text{*}}\right)\right)$ and observed location $s_c\mid s_c^{\text{*}}$ according to the DHS jittering scheme.

We fit the models UnAdj, Smoothed and FullAdj described in Section 3 using an intercept and the five covariates described above. The PC prior on the Matérn GRF is specified by $\text{P}\left({\sigma }_{\text{S}}>S_0\right)=$ 0.05 and $\text{P}\left({\rho }_{\text{S}}>R_0\right)=0.50$ , where $S_0=1$ is the 95th percentile of the marginal standard deviation, and $R_0=160\text{ km}$ is the median range. Inference is performed as described in Section 3.3.

Parameter estimation is evaluated by computing the RMSE, $\frac{1}{n_{\text{sim }}}{\sum }_{b=1}^{n_{\text{sim }}}{\left({\widehat{\theta }}^{\left(b\right)}-\theta \right)}^2$ , and the bias, $\frac{1}{n_{\text{sim }}}{\sum }_{b=1}^{n_{\text{sim }}}({\widehat{\theta }}^{(b)}-\theta )$ , where ${\hat{\theta }}^{\left(b\right)}$ is the posterior mean (or MAP in the case of ${\rho }_{\text{S}}$ and ${\sigma }_{\text{S}}^2$ ) for dataset $b$ and $\theta$ is the true value of the coefficient. Predictions are evaluated on a fixed set of 1,000 randomly selected locations within Nigeria, where we predict $\eta \left(s\right)=\text{logit}\left(r\left(s\right)\right)$ with the posterior median. These predictions are evaluated by the average RMSE and CRPS defined by ${\int }_{{ℝ}^2}{(F\left(x\right)-\mathbb{I}\left(y\le x\right))}^2\text{ d}x$ , where $y$ is the true value and $F\left(\cdot \right)$ is the predictive distribution.

The simulation study was run on a computing server that operates on Linux (Ubuntu 20.04). The computing server has 28 cores (2X14-core Xeon 2.6GHz) and provides 256GB memory limit per user. It took on average 4 minutes to estimate the parameters ${\sigma }_{\text{S}}$ and ${\rho }_{\text{S}}$ for FullAdj for each one of the 150 datasets, which is the most time-consuming part of the inference described in Section 3.3.1.

Table 3 shows the bias and RMSE for each of the parameters and covariate coefficients for SignalMed. The results show that UnAdj and Smoothed performs almost the same for ${\sigma }_{\text{S}}^2$ and ${\rho }_{\text{S}}$ , but that Smoothed gives higher bias and RMSE for the coefficients of the covariates. The latter is a consequence of the fact that smoothing the covariates changes the interpretation of the coefficients.

Table 3

Bias and RMSE of parameter estimation for SignalMed.

		Parameter
	Model	${\rho }_{\mathrm{S}}$	${\sigma }_{\mathrm{S}}^2$	$\mu$	${\beta }_{\text{DistW }}$	${\beta }_{\text{CityA }}$	${\beta }_{\text{Elev }}$	${\beta }_{\text{PopD }}$	${\beta }_{\text{UrbR }}$
Bias	UnAdj	-21.32	0.05	0.32	0.00	0.08	0.01	-0.16	0.84
	Smoothed	-19.97	0.05	0.21	0.14	-0.22	0.04	0.15	-1.03
	FullAdj	7.51	0.009	0.27	0.01	0.02	0.02	-0.09	0.37
RMSE	UnAdj	22.48	0.09	0.42	0.27	0.09	0.19	0.16	0.85
	Smoothed	21.37	0.09	0.34	0.37	0.24	0.29	0.17	1.17
	FullAdj	11.96	0.09	0.37	0.26	0.04	0.17	0.09	0.42

For example, the association with the urbanicity ratio of a 250 m x 250 m pixel and the urbanicity ratio of a 5 km x 5 km pixel are two different things. Additional tables for SignalLow and SignalHigh are found in Section 3 in the Supplementary Material.

Figure 3 shows boxplots of estimated ${\rho }_{\text{S}}$ and ${\beta }_{\text{UrbR}}$ underlining that difference between models is also clear when you consider variation between simulations within the same scenario. There is a clear trend to underestimate the spatial range for UnAdj and Smoothed, and the association with the covariate is too weak for UnAdj. For Smoothed, the model is estimating a coefficient with different interpretation than for the other models due to the averaging of covariate rasters over $5\text{ km}\times 5\text{ km}$ windows. The differences are larger for stronger covariate signal levels.

Figure 3

Box plots of estimated ${\rho }_{\mathrm{S}}$ and ${\beta }_{\mathrm{U}\mathrm{r}\mathrm{b}\mathrm{R}}$ for SignalLow, SignalMed and SignalHigh. The horizontal red lines show the true parameter value.

Figure 4 shows the variation in RMSE and CRPS across datasets for predictions. FullAdj and UnAdj perform almost the same in both predictive metrics for SignalLow, but FullAdj is slightly better for SignalMed, and substantially better for SignalHigh. This indicates that the stronger the signal of the spatial covariates, the larger the gain from adjusting for jittering. The figure also indicates that Smoothed does not have better predictive ability than UnAdj.

Figure 4

Box plots of CRPS and RMSE for predictions for SignalLow, SignalMed and SignalHigh.

In addition to the simulation study described above, we performed a simulation study with the true locations fixed to the observed locations of NDHS2018. The goal was to examine the predictive metrics for the specific spatial design of NDHS2018. There were only minor differences between this simulation study, and the one described above, and we, therefore, show these results in Section 2 of the Supplementary Material. This suggests that we should expect the unadjusted model to underestimate range, and lose predictive accuracy due to reduced association with the covariates.

5 Analysis of completion of secondary education

In this section, we analyse the spatial variation in prevalence of completion of secondary education among 20-49 year old women in Nigeria in 2018 based on the NDHS2018 using the data sources described in Section 2. Our analysis has two aims. First, to map the fine-scale spatial variation in completion of secondary education for 5 km x 5 km pixels, and for the admin 1 areas. Second, to determine the associations between the spatial variation in risk and a set of explanatory spatial covariates.

The NDHS2018 has C = 1,380 clusters with jittered GPS coordinates available under the same jittering distribution as in Section 3.2. For all clusters, the jittering was restricted to stay within the correct admin2 area. In total, 33,193 women aged 20-49 years were interviewed and 15,490 of these had completed secondary education. We use the notation $n_c$ individuals-at-risk, $y_c$ successes, and jittered GPS coordinate ${\mathit{s}}_c$ for $c=1,\dots ,C$ . The five covariates of interest are introduced in Section 2.

We fit the models UnAdj, Smoothed and FullAdj described in Section 3 using an intercept and the five covariates described in Section 2. We use the PC prior on the Matérn GRF specified by $\mathrm{P}\left({\sigma }_{\mathrm{S}}>S_0\right)=0.05$ and $\mathrm{P}\left({\rho }_{\mathrm{S}}>R_0\right)=0.50$ , where $S_0=1$ is the 95th percentile of the marginal standard deviation, and $R_0=160\mathrm{k}\mathrm{m}$ is the median range. Inference is performed as described in Section 3.3.

Table 4 shows the estimated parameters and their corresponding credible interval lengths (except for ${\rho }_{\mathrm{S}}$ and ${\sigma }_{\mathrm{S}}^2$ , which are fixed to their MAP estimates). The results show that ${\rho }_{\mathrm{S}}$ is estimated substantially smaller for UnAdj and Smoothed than for FullAdj. Based on the results in Section 4, this indicates that spatial correlation is lost by not accounting for jittering. Additionally, ${\sigma }_{\mathrm{S}}^2$ is estimated higher for UnAdj and Smoothed than for FullAdj. This could be because the covariates are able to explain less of the spatial variation for the former two models.

Table 4

Parameter estimates and the corresponding 95% credible interval lengths in parentheses. Uncertainty is not computed for ${\rho }_{\mathrm{S}}$ and ${\sigma }_{\mathrm{S}}^2$ .

	Parameter
Model	${\rho }_{\mathrm{S}}$	${\sigma }_{\mathrm{S}}^2$	$\mu$	${\beta }_{\text{DistW }}$	${\beta }_{\text{City }}$	${\beta }_{\text{Elev }}$	${\beta }_{\text{PopD }}$	${\beta }_{\text{UrbR }}$
UnAdj	64.15	1.91	$-2.24\left(0.89\right)$	$0.91\left(1.24\right)$	$-0.40\left(0.13\right)$	$-0.14\left(0.69\right)$	$0.14\left(0.09\right)$	$-0.14\left(0.42\right)$
Smoothed	58.42	1.88	$-2.32\left(0.82\right)$	$1.08\left(1.43\right)$	$-0.78\left(0.29\right)$	$-0.40\left(0.88\right)$	$0.49\left(0.30\right)$	$-2.97\left(2.20\right)$
FullAdj	107.68	1.65	$-2.21\left(1.04\right)$	$0.62\left(1.29\right)$	$-0.43\left(0.17\right)$	$-0.02\left(0.61\right)$	$0.32\left(0.13\right)$	$-1.35\left(0.69\right)$

For PopD and UrbR, there is a strong attenuation when jittering is ignored. The credible intervals for the coefficient of UrbR, ${\beta }_{\mathrm{U}\mathrm{r}\mathrm{b}\mathrm{R}}$ , suggest that ${\beta }_{\mathrm{U}\mathrm{r}\mathrm{b}\mathrm{R}}$ is not significant at the $95\mathrm{\%}$ level for UnAdj, whereas ${\beta }_{\text{UrbR }}$ is significant for FullAdj. This suggest that not accounting for jittering can lead to misleading conclusions. As discussed in Section 4, it is hard to do direct comparisons in the estimated coefficients for Smoothed and for UnAdj and FullAdj since averaging covariates changes their meaning.

Figure 5 shows the $5\mathrm{k}\mathrm{m}\times 5\mathrm{k}\mathrm{m}$ pixel maps of predicted risk and CVs for UnAdj, Smoothed and FullAdj. The results show that some areas such as Borno (in the north-east) have up to three times the risk under the UnAdj approach as under FullAdj. Further, UnAdj tends to lead to higher uncertainty in the predictions than FullAdj.

Figure 5

Rows 1, 2 and 3 are predicted risk and the CVs for UnAdj, Smoothed, and FullAdj, respectively. Row 4 shows ratios (UnAdj/FullAdj) of predictions and CVs.

We aggregate point level predictions with respect to population density to produce areal estimates at the $37\mathrm{a}\mathrm{d}\mathrm{m}\mathrm{i}\mathrm{n}1$ areas (for more on aggregating point level predictions with respect to a population, see Paige et al. 2022). Figure 6 shows the predicted risk and associated CVs for UnAdj, Smoothed and FullAdj. From Figure 6(g), we see that the point estimates vary from a factor 0.9 to 1.2, and Figure 6(h) shows that some areas differ with a factor of up to 1.4 in CVs.

Figure 6

Rows 1 and 2 are predicted risk and CVs for UnAdj and FullAdj respectively at the admin 1 level, and row 3 shows ratios (UnAdj/FullAdj) of predictions and CVs.

In Figure 7, we compare the areal estimates from the FullAdj model against direct estimates. Figure 7(b) shows that the FullAdj reduces uncertainty compared to the direct estimates, and Figure 7(a) shows that the estimates from FullAdj and the direct estimates are compatible with no unexpected deviations.

Figure 7

Comparison direct estimates and (areal) model-based estimates for (a) the mean, and (b) the coefficient of variations.

The ability to predict risk at unobserved locations for UnAdj, CovAdj and FullAdj cannot be compared with cross validation. If data is held-out from NDHS2018, we can only evaluate the models' ability to predict risk at a new jittered cluster with unknown true location. However, the simulation study in Section 4 indicates that FullAdj leads to an improvement in prediction if the data-generating model is the one estimated in this section.

6 Discussion

Accounting for jittering substantially changed the parameter estimates for the geostatistical model for completion of secondary education among women aged 20-49 years. The simulation study demonstrated that these differences were linked to the strength of the signal of the spatial covariates when explaining the spatial variation. For strong signals, the associations were attenuated and the predictive power was reduced.

The most important aspect of jittering in the context of geostatistical models for DHS data is to account for the resulting uncertainty in covariates extracted from rasters or extracted based on distances. This induces measurement error that may lead to attenuation in associations between covariates and the responses. Some covariates such as sanitation practices and household assets can be known exactly (Burgert-Brucker et al., 2016), but these cannot be included when the goal is prediction since fine-scale rasters are not available.

In influential work such as Burstein et al. (2019) and Local Burden of Disease Vaccine Coverage Collaborators (2021), covariates are resampled to a $5\mathrm{k}\mathrm{m}\times 5\mathrm{k}\mathrm{m}$ grid. This is similar, but slightly different than the Smoothed model discussed in this article which uses $5\mathrm{k}\mathrm{m}\times 5\mathrm{k}\mathrm{m}$ windows around the observed locations. However, the results suggest that such approaches that average covariates do not address the issue of jittering, and do not improve predictions.

This article used uniform priors for the unknown true locations. One could expect including information about population density and urbanicity into the priors would produce more accurate inference. However, population density maps and urbanicity maps are also modelled surfaces with biases and uncertainties that are not well understood. This means that evaluation of the sensitivity to such maps would have to be investigated, and one would need a way to evaluate whether such a model works better.

The inference scheme uses empirical Bayes. It is possible to investigate methods such as INLA, but the implementation in the $\mathrm{R}$ ( $\mathrm{R}$ Core Team, 2022) package inla does not allow the likelihood to depend on the latent risk at multiple locations, which is necessary due to the integration points. MCMC algorithms such as STAN (Stan Development Team, 2020) has the required flexibility, but is infeasible for thousands of spatial locations.

When analysing completion of secondary education, we found that, for urbanicity, an effect size of 0 was contained in the $95\mathrm{\%}$ credible interval when not adjusting for jittering, and not contained when adjusting for jittering. This suggests that not accounting for jittering when analysing DHS data is a practice that can alter conclusions about statistical significance. Since the proposed approach is fast for spatial analysis, we suggest using the new approach routinely for analysing DHS data to avoid the risk of misleading conclusions and reduced predictive power. The code used throughout the article can be found in the GitHub repository https://github.com/umut-altay/GeoAdjust.

Supplementary material

Supplementary material for this article is available online.

Supplemental Material for Impact of jittering on raster- and distance-based geostatistical analyses of DHS data by Umut Altay, John Paige, Andrea Riebler and Geir-Arne Fuglstad, in Statistical Modelling

Appendix

Computational details

$\mathrm{R}$ ( $\mathrm{R}$ Core Team, 2022) is the main programming language used for the code that produced the results in this manuscript. The computational part of the study utilizes auto-differentiation via TMB (Kristensen et al., 2016, 2023). This means that low-level computations are done according to C++ code while high-level computations are done in $\mathrm{R}$ . The code scripts that are used for this article can be found in the GitHub repository https://github.com/umut-altay/GeoAdjust.

The application in Section 5 was run on a MacBook Pro (2.4GHz Quad-Core Intel Core i5 and 16 GB memory). For FullAdj, it took 7.8 minutes to estimate parameters, and 4.2 minutes to compute predictions.

A newly developed R-package GeoAdjust (Altay et al., 2023) is now available on CRAN for analysing DHS data while accounting for jittering in a more compact way, similar to the analysis in Section 5. The package provides a user-friendly functionality by isolating the user from $\mathrm{C}++$ and complex R code. It is flexible on allowing analysing data with Gaussian, binomial and Poisson responses as well.