Prerpints.org logo
Preprint
Article

This version is not peer-reviewed.

An Age-Structured Spatially Varying Coefficient Model for High Resolution Mapping of Vaccination Coverage

Submitted:

16 July 2025

Posted:

17 July 2025

You are already at the latest version

Abstract
High-resolution maps of vaccination coverage are valuable for uncovering heterogeneities in coverage to inform vaccine delivery strategies. Coverage maps stratified by age can reveal additional heterogeneities in the timeliness of vaccination and critical immunity gaps among birth cohorts. Here, we propose a spatially varying coefficient model relying on a Bayesian approach for age-structured mapping of vaccination coverage using geolocated individual level household survey and geospatial covariate data. Our flexible modelling framework includes parameterizations capturing spatial (non-)stationarity in differences in coverage between age groups, as well as a modification to allow coverage mapping for singe age points through the inclusion of a smoother over age. The proposed models are fitted using the INLA-SPDE approach implemented in the inlabru package in R. We choose between competing model parameterizations by examining their out-of-sample predictive performance via cross-validation and using Bayesian model choice criteria. The methodology is applied to age-structured mapping of measles vaccination coverage in Cote d’Ivoire using the 2021 Demographic and Health Survey. Our results reveal a significant delay in measles vaccination in the first year of life and substantial spatial differences in coverage by age, highlighting the need for targeted interventions to achieve equity and attain vaccine-derived immunity goals.
Keywords: 
vaccination coverage
;  
measles
;  
spatially varying coefficient model
;  
INLA-SPDE
;  
demographic and health survey
Subject: 
Public Health and Healthcare  -   Primary Health Care

Introduction

Vaccination coverage is an important metric for monitoring the performance of immunization programmes [1,2]. Accurate estimates of vaccination coverage are essential at global, national, and subnational levels to assess the relative effectiveness of vaccine delivery strategies such as routine immunization and campaigns [3,4], to guide and evaluate targeted interventions [5], and to track progress toward national and global immunization goals [1,6,7], among other applications. For certain vaccines such as the measles-containing vaccine, estimates of coverage by age are particularly valuable since achieving coverage targets across all age groups and geographic areas is integral to disease elimination strategies [8]. Such estimates enable the monitoring of the timeliness of vaccination [9,10], identification of key immunity gaps across birth cohorts and the design of catch-up vaccination activities [11,12,13], thereby equipping national immunization programmes with the evidence needed to develop and implement effective approaches to achieve and sustain disease elimination.
Despite their importance, reliable estimates of vaccination coverage remain challenging to obtain, particularly at fine geographic scales. Data for estimating coverage are usually obtained from administrative sources and household surveys. The quality of data available from these sources can vary greatly due to inherent limitations such as incomplete reporting, recall bias and inaccuracies in population denominators, with survey estimates of coverage often taken to be the “gold standard” [2]. However, survey estimates are mostly produced for large administrative areas such as the national and first administrative levels due to the operational costs and logistical challenges of conducting more intensive sampling to produce estimates for more granular levels. These large area estimates mask epidemiologically and programmatically important heterogeneities in coverage. In line with increasing emphasis on more granular coverage estimates, as encapsulated in global policy frameworks such as the Immunization Agenda 2030 [1], there have been growing efforts to estimate coverage at more operationally relevant geographic scales to better inform vaccination programming [14,15,16,17,18]. Furthermore, survey-based estimates are typically reported for standard age groups, as determined by relevant vaccination schedules [7]; and when surveys are conducted infrequently, these estimates often fail to capture critical immunity gaps in other eligible age groups.
Few published studies such as Utazi et al [17] and Takahashi et al [15] have developed approaches for age-structured, high-resolution mapping of vaccination coverage using geolocated survey data. Utazi et al [17] employed a birth cohort approach, fitting independent geostatistical models to predict measles vaccination coverage in children aged 9-11 months, 12-23 months and 24-35 months. However, their methodology provides no mechanism to exploit shared patterns - such as the underlying spatial structures - across the age groups, which could potentially improve coverage estimation. Moreover, because separate models were fitted for each age group, their approach may perform poorly in contexts where overall or cluster-level sample sizes for a given age group are small (see, for example, the data section), since this can result in greater uncertainties in the modelled estimates for the age group [19]. Takahashi et al [15] adopted a single age approach, using a generalized additive model to map measles vaccination coverage by age in the African Great Lakes region. Their method incorporates a smooth function of age, enabling estimation of coverage at single-month age points but not for broader age groups. Additionally, the inclusion of a smooth function of survey cluster-level geocoordinates in their model primarily captures broad spatial trends and may be inadequate for modelling structured spatial dependence [20].
Spatially varying coefficient (SVC) models [21,22] are particularly useful in spatial regression contexts when relationships between an outcome and certain covariates exhibit spatial non-stationarity. These models allow the regression coefficients associated with these covariates to vary across space, yielding spatial surfaces/maps for these, and are also suitable for capturing differences in space between the levels of a given categorical covariate in relation to an outcome (see methods section). Gelfand et al [21] introduced a Bayesian framework for SVC models and these have been used in ecological settings in both discrete [23] and continuous [22] spatial domains to uncover nuanced spatial variation in ecological processes such as bird species abundance and forest biomass. In particular, Finley [22] compared different SVC model specifications, including geographically weighted regression [24], using both simulated and real-world datasets. The study accounted for directional (anisotropic) and nonstationary dependence structures, demonstrating the flexibility of these models to capture complex spatial patterns in the data.
In this study, we develop a novel and flexible methodology for age-structured mapping of vaccination coverage. The proposed approach is based on a spatially varying coefficient model that captures spatial non-stationarity in the differences in coverage between age groups, with alternative model specifications characterizing different forms of age-related variation, and which also enable the estimation of coverage using the single age approach. Unlike previous work [17,19,25], the proposed models are fitted to individual level data, enabling more robust estimation of covariate-coverage relationships and potentially improving estimation in areas with few cluster-level observations. Model selection is based on Bayesian model choice criteria and using cross-validation to evaluate out-of-sample predictive performance. The methodology is applied to age-disaggregated mapping of the coverage of the first dose of the measles-containing vaccine (MCV1) in Cote d’Ivoire, using data from the 2021 Demographic and Health Survey [26]. Cote d’Ivoire was chosen as a case study due to being a priority country for an implementing/donor organization and has therefore been the focus of recent work [27].

Methods

Vaccination Coverage Data

Geolocated individual-level data on the coverage of MCV1 were obtained from the 2021 Cote d’Ivoire DHS (CDHS) [26]. The survey was designed to be representative at the national and first administrative (14 areas known as districts – supplementary Figure 1) levels, and for urban and rural areas. The sampling design involved a random, stratified two-stage selection, with urban and rural areas of each of the 14 districts constituting separate sampling strata, yielding a total of 28 strata. The first stage involved selecting the primary sampling units known as clusters (or enumeration areas) from a sampling frame, while the second stage involved selecting households within each sampled cluster, with 28 households sampled per cluster. In all, a total of 539 clusters were included in the survey [26].
Routine MCV1 is recommended for administration at 9 months of age in Cote d’Ivoire [28]. Hence, for each eligible child aged 9 to 35 months, we extracted information on their age, vaccination status (based on health records or caregiver recall) and the (latitude/longitude) geocoordinates of the cluster they were sampled from. Each individual level record was classified into one of three age groups - 9-11 months, 12-23 months or 24-35 months - representing different birth cohorts, and aggregated to the survey cluster level for validation purposes. The data extracted contained a total of 4,236 children, of whom 507, 1,920 and 1,809 belonged to the respective age groups. A summary of the extracted data is displayed in Figure 1 and supplementary Figure 2. This includes national survey weighted and unweighted estimates of MCV1 coverage for single age points and age groups, unweighted coverage estimates at the cluster level for age groups and corresponding distribution of cluster level sample sizes for each age group.
These plots generally reveal a sharp increase in coverage at lower ages, which then appears to stabilize and exhibit random fluctuations randomly in the later ages. The relatively lower coverage among children aged 9-11 months is evident both at the national and cluster levels. Cluster level sample sizes are much smaller for age 9-11 months, limiting the feasibility of a separate analysis for this age group. For the other age groups, cluster level sample sizes are also small (median = 3, interquartile range = 3). The distribution of these cluster-level sample sizes is particularly important for accurately estimating the underlying spatial structure in the data and for obtaining more robust point or grid level predictions of vaccination coverage across age groups.

Geospatial Covariates, Gridded Population and Boundary Data

Following previous studies [3,17,19,25], we assembled a wide range of geospatial covariate datasets for the analysis. These include geographic, socioeconomic, remoteness, climate, environmental, conflict and health-related variables, which have either direct or proximate associations with vaccination coverage. While most of the covariates were externally sourced, others such as ownership of a health card/document were obtained from the 2021 CDHS. These CDHS-derived covariates were first processed at the cluster level using the definitions provided in supplementary Table 1. Kriging interpolation was then used to create their interpolated surfaces at 1x1 km resolution as outlined in previous work [19,29]. Covariate selection was undertaken in a non-spatial framework, using a similar approach as outlined in Utazi et al [17]. This involved examining the relationships between the covariates and vaccination coverage, and applying log transformations to the covariates where appropriate; fitting single-covariate models and ranking the covariates based on their predictive ability (using predictive R2 values); assessing multicollinearity and choosing between highly correlated covariates (with correlation > 0.8 or variance inflation factor > 4.0) based on their ranks; and finally, using stepwise regression (backward elimination based on the Akaike Information Criterion) to identify the best model or combination of covariates for modelling the outcome in a non-spatial framework with binomial regression models. A total of 13 covariates were selected for the analysis from an initial set of 38 covariates, including urban and rural areas, elevation, vegetation index (veg_index, average number of wet days (wetdays), urban accessibility (urban_access), elevation, maximum temperature (max_temp), distance to conflict areas (dist_conf), walking travel time to the nearest health facility (walking_tt), malaria prevalence (mal_prev), maternal education (mat_educ), health card/document ownership (health_card), access to media (media) and household wealth (wealth) (see supplementary Table 1 for details). These selected covariates are displayed in supplementary Figure 3 and Figure 4.
Gridded population data for children aged under 5 years in 2021 were obtained from WorldPop [30]. These were used as weights when aggregating the grid level estimates to the administrative level [19]. All boundary data sets used in the analysis were obtained from geoBoundaries [31].

Model Development

We begin by introducing the spatially varying coefficient model and its parameterizations for estimating vaccination coverage using the birth cohort approach, and then alternative specifications that enable estimation at single age points and for all eligible individuals.
Let z i ( s j ) denote the age (in months) of the i th eligible individual/child sampled from location s j , and a 0 , a 1 ,   a 2 , … , a K , the different K + 1 age groups (or birth cohorts), one of which includes z i ( s j ) . For example, in our application, K = 2 , and a 0 = 9 ,   10 ,   11   m o n t h s ; a 1 = { 12 , … ,   23   m o n t h s } and a 2 = { 24 , … ,   35   m o n t h s } . Also, let Y i s j denote the corresponding binary outcome representing the individual’s vaccination status i . e . ,   Y i s j = 1   i f   v a c c i n a t e d   a n d   0   o t h e r w i s e , and p i s j the (unobserved) true probability of being vaccinated. The spatially varying coefficient model can be expressed as
Y i s j | p i s j ∼ B i n o m i a l 1 , p i s j ,   i = 1 , … , n j ; j = 1 , … , J ,
l o g i t p i s j = β 0 + ∑ k = 1 K β k + β k ( s j ) I ( z i ( s j ) ∈ a k ) + x s j ' α + ω s j + ϵ ( s j ) ,
where β 0 is the overall intercept, β k and β k ( s j ) are fixed and spatially varying regression coefficients for age cohort a k ( k > 0 ) , I . is an indicator function that returns a value of 1 if the specified condition is true and 0 otherwise, and x ( s j ) is a vector of geospatial covariates for location s j , with α being the corresponding fixed regression coefficients. Thus, β k ( s j ) can be viewed as a random spatial adjustment to the average change in the log-odds of vaccination for the k th age cohort, β k , at location s j . Note that the reference age group, a 0 , is omitted from equation (1) as this is absorbed into the intercept terms. In practice, the choice of reference age group may vary depending on the objective of the analysis. Also, the covariate vector, x s j ,   includes an urban-rural variable which helps to account for the urban-rural stratification used in the survey design. Only cluster-level geospatial covariates are included in the model since the goal of inference is to obtain age-specific coverage estimates at the grid level, typically at 1x1 km resolution.
Furthermore, the terms ω s j and ϵ ( s j ) are modelled as Gaussian random effects, capturing residual spatial and non-spatial (or between-cluster) variation, respectively. Specifically, ω s j represents a random spatial adjustment to the overall intercept, β 0 , such that β 0 + ω s j defines a spatially varying intercept process (for the reference age group). We model ϵ ( s j ) as an independent and identically distributed (iid) process with mean 0 and variance σ ϵ 2 . The spatial random effect ω s j is specified as a zero-mean Gaussian process with covariance function Σ ω , that is, ω = ω s 1 , … , ω ( s J ) ' ∼ N ( 0 ,   Σ ω ) , where Σ ω is assumed to follow the Matérn covariance function [32] given by
Σ ω ( s j , s j ' ) = σ ω 2 2 ν − 1 Γ ν κ ∥ s j − s j ' ∥ ν   K ν   ( κ ∥ s j − s j ' ∥ ) .
Here, the notation ∥ . ∥ denotes the Euclidean distance ,    σ ω 2 > 0 is the marginal variance of ω s j , κ > 0 is a scaling parameter related to the range r ( r = 8 ν κ ) – the distance at which spatial correlation is close to 0.1 - and K ν is the modified Bessel function of the second kind and order ν > 0 . Following Lindgren et al [33], it is common to fix the smoothness parameter ν to one to ensure identifiability.
For each of the spatially varying regression coefficients, β k s j ,   k = 1 , … K , we also assume zero-mean Gaussian processes with Matérn covariance functions denoted by Σ β k . We denote the corresponding range and marginal variance parameters using r k and σ β k 2 , respectively.
The general model in equation (1) can be further parameterized to yield simpler, more parsimonious models, depending on the nature of the dependencies assumed in the data. For example, if the spatially varying coefficients and the geospatial covariates are deemed sufficient to capture all the spatial dependencies, or more specifically, if the geospatial covariates are assumed to fully account for the spatial variation in the outcome for the reference age group, then model (1) reduces to
l o g i t p i s j = β 0 + ∑ k = 1 K β k + β k ( s j ) I ( z i ( s j ) ∈ a k ) + x s j ' α + ϵ ( s j )
which excludes the spatially correlated residual term, ω s j . Other variants of model (2) are also possible. For instance, if prior evidence suggests that step changes captured by β k are sufficient to characterise the differences between certain age groups and the reference age group, then the corresponding β k ( s j ) terms may also be excluded, while retaining ω s j . Such specifications may be particularly advantageous in contexts involving a large number of age groups.
Furthermore, if the differences in the odds of vaccination between each age group and the reference age group are assumed not vary spatially and can be adequately characterised using fixed step changes, then the spatially varying coefficients can be removed entirely from model (1), yielding
l o g i t p i s j = β 0 + ∑ k = 1 K β k I ( z i ( s j ) ∈ a k ) + x s j ' α + ω s j + ϵ ( s j )
In practice, models (2) and (3) provide alternative model specifications against which to compare model (1) to determine the most appropriate modelling approach for capturing spatial dependencies in the data.
Additionally, when the goal is to predict coverage for single ages (in months), model (3) can be modified to include a smooth function of age, f ( z i ( s j ) ) :
l o g i t p i s j = β 0 + f ( z i ( s j ) ) + x s j ' α + ω s j + ϵ ( s j ) ,
where f ( . ) is modelled using a second-order random walk: f z r z r − 1 , z r − 2 ∼ N ( 2 z r − 1 − z r − 2 , σ z 2 ) , which corresponds to the Bayesian equivalent of a cubic smoothing spline [34]. A sum-to-zero constraint was imposed on f ( . ) for identifiability, since the model includes an intercept term [34]. We note that, as in the birth cohort model, it is possible to model the effect of age in equation (4) using a spatially varying coefficient, such that f z i s j = β s j z i ( s j ) . However, such parameterization is not ideal for estimating the overall smoothed effect of coverage by age as intended in model (4). Hence, model (4) captures the non-linear effect of age on vaccination but assumes constant residual spatial dependence for all ages as in model (3). Similar to model (3), model (4) can also be used to obtain predictions for age groups if desirable (see discussion section).
Finally, to predict coverage for all eligible age groups combined, the age-specific terms in model (1) can simply be dropped, yielding
l o g i t p i s j = β 0 + x s j ' α + ω s j + ϵ ( s j ) .
Model (5) aligns with approaches commonly used in previous work to map vaccination coverage [16,19,25], except that it is estimated using individual rather than cluster level data. Note that in equation (5), as in all the models considered here, individual level variation is already accounted for in the binomial likelihood, with the variance of the logistic distribution fixed at π 2 3 ≅ 3.29 [35].
For conciseness, we subsequently refer to the candidate birth cohort model specifications in equations (1) to (3) as MODsvc1, MODsvc2 and MODnosvc, respectively. Similarly, we refer to the single-age and all-age models in equations (4) and (5) as MODsmooth and MODall, respectively.

Bayesian Inference, Model Fitting and Prediction

The proposed models are implemented in a fully Bayesian framework. We place the following prior distributions on the parameters:
β = ( β 0 ,   β 1 , … , β K ) ' ∼ N 0 , 10 3 I ;   α ∼ N ( 0 , 10 3 I ) ;   p r < r ~ = 0.01 ;   p r k < r ~ = 0.01 ,
p σ ϵ > 3 = 0.01 ;   p σ ω > 3 = 0.01 ;   p σ β k > 3 = 0.01 ;   σ z − 2 ∼ G a m m a ( 1 , 0.01 )
These comprise non-informative priors on the fixed regression parameters and a weakly informative prior on the precision of the second-order random walk in MODsmooth. The priors on the range parameters and the remaining variance parameters are penalized complexity priors, introduced by Simpson et al [36] and used in Fuglstad et al [37] in a spatial modelling context, which penalize deviations from a base (simpler) model, e.g., a constant field. For the range parameters, r ~ is chosen to be 5% of the extent of the study area in the longer direction (north-south direction for the Cote d’Ivoire application).
Let θ = ( β 0 , β 1 , … , β K , α , r , r 1 , … , r K , σ ϵ 2 , σ ω 2 , σ β k 2 ) ' denote all the fixed parameters of MODsvc1 in equation (1), β k s = β k s 1 , … , β k ( s J ) ' , ϵ = ϵ s 1 , … , ϵ ( s J ) ' and y denote all the observed data. Its joint posterior distribution can be expressed as
π θ , ω , β 1 s , … , β K s , ϵ y ∝ π θ × N ω | 0 ,   Σ ω × N ϵ | 0 ,   σ ϵ 2 I × N β 0 | 0 ,   10 3 × ∏ k = 1 K N β k s | 0 ,   Σ β k × N β k | 0 ,   10 3     × ∏ j = 1 J ∏ i = 1 n j B i n o m i a l ( y i s j ;   p i s j , z i ( s j ) , y )
where π ( θ ) is the joint prior distribution on θ . The posterior distributions of the other models can be straightforwardly derived from equation (7).
For parameter estimation, we use the integrated nested Laplace approximation-stochastic partial differential equation (INLA-SPDE) approach [33,38] implemented in the inlabru package in R [39]. The INLA approach performs approximate Bayesian inference for latent Gaussian models by using the Laplace method to approximate the marginal posterior distributions of all the unknown quantities in the model, including the parameters and latent variables. The SPDE approach facilitates the estimation of the spatial terms in the models (e.g., ω and β k s in model (1)). The approach works by representing each of these terms as a Gaussian Markov random field to reduce the computational burden inherent in the estimation of their respective covariance matrices [33].
Of particular interest in the proposed models is the prediction of vaccination coverage for the age groups or single age points of interest at unsampled locations, e.g. 1 × 1 km grids, throughout the study area. In a Bayesian framework, the prediction of vaccination coverage for age group k ∈ { 0,1 , … , K } at a new location, s ~ , proceeds by evaluating the posterior predictive distribution of p k s ~ given by
π p k s ~ y = ∫ l o g i t − 1 π η k s ~ ω , ϵ , β k s , θ , y × π ω , ϵ , β k s , θ y   d ω d ϵ d β k s d θ ,
where η k s ~ is the linear predictor for MODsvc1, such that
η k s ~ = β 0 + x s ~ ' α + ω s ~ + ϵ s ~                              i f   k = 0 β 0 + β k + β k s ~ + x s ~ ' α + ω s ~ + ϵ s ~    i f   k ≥ 1 .
In the inlabru package, equation (8) can be computed post model fitting by using either the g e n e r a t e ( ) or p r e d i c t ( ) function. We obtained 1,000 posterior samples from equation (8) for each prediction grid location, which were then summarized to produce the grid level estimates and associated uncertainties. Administrative level estimates were also obtained as population-weighted averages over all the grid cells falling within each administrative unit [19], as noted previously, using these posterior samples. Individual level predictions (in- and out-of-sample prediction for validation purposes) and predictions using models (2), (3), (4) and (5) can also be obtained in a similar manner.

Model Choice and Validation

To compare the fits of the competing birth cohort models in equations (1) to (3) (and those of MODsmooth and MODall), we employed the Watanabe-Akaike Information Criterion (WAIC) [40,41] which uses the log of the predictive density of each data point to evaluate a model’s out-of-sample predictive performance. The WAIC can be expressed as
W A I C = − 2 ∑ i = 1 n log ⁡ π y i y + 2 p W A I C ,
where π y i y = ∫ − ∞ ∞ π y i θ π θ y d θ is the posterior predictive density of data point y i (out of n data points), which can be evaluated in practice as 1 L ∑ l = 1 L π y i θ l with θ l denoting the l -th posterior draw from a total of L posterior samples. The term p W A I C represents the effective number of parameters, computed as p W A I C = ∑ i = 1 n V a r θ | y [ l o g ⁡ ( π ( y i | θ ) ) ] . The WAIC balances model fit and complexity by penalizing overfitting through p W A I C , which captures the sensitivity of the model’s fit to individual data points. A lower WAIC indicates better predictive accuracy.
To further evaluate the out-of-sample predictive performance of the proposed models both at the cluster and individual levels, we conducted a k -fold cross-validation exercise, creating the cross-validation folds as random and spatially stratified splits of J cluster locations and setting k = 10 . At the cluster level, we used the following metrics to evaluate predictive performance:
Root mean square error, R M S E = 1 m ∑ j = 1 m p ^ s j − p s j 2 ;
Mean absolute error, M A E = 1 m ∑ j = 1 m p ^ s j − p ( s j ) ;
Average bias, A V G _ B I A S = 1 m ∑ j = 1 m p ^ s j − p s j ;
Continuous ranked probability score,
C R P S = 1 m ∑ j = 1 m 1 S ∑ s = 1 S p ^ ( s ) s j − p ( s j ) − 1 2 S 2 ∑ s = 1 S ∑ t = 1 S p ^ ( s ) s j − p ^ ( t ) ( s j ) ;
where m is the number of cluster-level observations in the k th cross-validation set, p s j denotes the observed coverage and p ^ s j the corresponding predicted values (i. e., the posterior means). With the CRPS [42], the entire S posterior samples are used to measure the discrepancies between the observed values and the predicted values for each data point rather than using a particular summary as in the other metrics. The smaller the values of these metrics, the better the predictive performance.
At the individual level, we evaluated out-of-sample predictive performance by computing the Brier score, B S = 1 M ∑ j = 1 m ∑ i = 1 n j p ^ i s j − y i s j 2 , where M = ∑ j = 1 m n j is the number of individual level observations in m validation set locations and y i ( s j ) denotes the binary outcomes [43]. A lower BS indicates better model performance. All the metrics were averaged over all k cross-validation folds in each case.

Results

Model Choice

Table 1 presents the WAIC values of the candidate birth cohort models (MODsvc1, MODsvc2 and MODnosvc), along with those of MODsmooth and MODall for comparison. Among the three birth cohort models, MODsvc1 had the best predictive performance. All models that account for age clearly outperformed MODall, demonstrating better predictive accuracy at the individual level and substantial variation in coverage by age within the data. Notably, the better performance of MODsmooth relative to MODnosvc suggests that, in the absence of the spatially varying coefficients, the smooth age function in MODsmooth more effectively captured the gradients in age-specific coverage than the step changes in MODnosvc. Also, we observe that although MODsvc1 and MODsvc2 - both incorporating spatially varying coefficients - had the lowest WAIC, these also had the largest effective numbers of parameters relative to other models as expected given their greater flexibility.
The patterns observed in the WAIC results are corroborated by the BS values reported in supplementary Table 2. These results show that for individual level out-of-sample prediction under the random cross-validation scheme, MODsvc1 and MODsvc2 had the best predictive performance, although the differences between the models were modest. Under stratified cross-validation, the differences between the models were even smaller, with MODsvc1 still ranking among the best-performing models. Consistent with previous findings, for both cross-validation schemes, all models accounting for variation by age (except MODsvc2 under stratified cross-validation) outperformed MODall.
For out-of-sample prediction at the cluster level, Table 2 reports the validation metrics for MODsvc1, MODnosvc and MODsvc2, while supplementary Table 3 provides corresponding results for MODall for comparison. Due to sample size limitations, it was not feasible to evaluate cluster-level predictive performance for MODsmooth. Across the three models, out-of-sample predictive performance was broadly similar within each age group, though MODsvc1 performed slightly better overall, particularly for the 9–11 and 24–35 month age groups. For the 12-23 month age group, the differences between the models were marginal, with no model clearly outperforming the others. Importantly, all three models performed better for age groups 12-23 months and 24-35 months compared to 9-11 months (except when considering AVG_BIAS, which poorly differentiated between the models), likely reflecting the relatively larger cluster level sample sizes for both older age groups (Figure 1 and supplementary Figure 2). Indeed, for the aggregated 9-35 month age group, the validation metrics reported in supplementary Table 3 are consistently lower than those in Table 2, indicating improved predictive performance associated with larger cluster level sample sizes. As expected, predictive performance across all the models was generally poorer under the stratified cross-validation scheme, reflecting the increased difficulty of prediction in this setting.
Overall, our evaluation of the out-of-sample predictive performance of the competing birth cohort models, at both the individual and cluster levels, indicates that MODsvc1 offers the best performance. Accordingly, in the following sections, we discuss the results from MODsvc1 when using the birth cohort approach to estimate vaccination coverage.

Parameter Estimates

In Table 3, we present the estimates of the parameters of the best-fitting birth cohort model, MODsvc1. The covariates significantly associated with the odds of vaccination include: urban versus rural residence, distance to conflict areas, walking travel time to the nearest health facility, maternal education, health card/document ownership, access to media and age group. As expected, the odds of vaccination increased with greater distance to conflict areas, higher proportions of educated mothers, higher proportions of children who owned a health card and greater access to the media within households. Conversely, the odds of vaccination decreased as walking travel time to the nearest health facility increased.
Considering the estimated odds ratios, a unit increase in the proportion of children who owned a health card is associated with a sixfold increase in the odds of vaccination, for example. Notably, children residing in urban areas had 46% lower odds of being vaccinated compared to those in rural areas – a pattern that appears to be specific to the study country [26]. Furthermore, on average, children aged 12-23 months and 24-35 months had approximately twice the odds of being vaccinated compared to those aged 9-11 months, suggesting significant delays in MCV1 vaccination during the first year of life.
The spatial variation in the differences in odds of vaccination between the age groups is illustrated by the maps of β 1 ( s ) and β 2 ( s ) in supplementary Figure 5. The estimated surface for β 2 ( s ) appears much smoother than that for β 1 ( s ) , suggesting greater spatial variability in the differences (in the log-odds of vaccination) between age groups 9-11 months and 12-23 months than between 9-11 months and 24-35 months. This provides a strong justification for the use of spatially varying coefficients in our modelling approach, but it also indicates that to preserve parsimony in the model, the difference in the likelihood of vaccination between ages 9-11 months and 24-35 months could be modelled as a step change, as discussed previously in the methods section. We discuss the spatial patterns in the differences in coverage between the age groups in detail below.
The estimated spatial ranges for the spatially varying coefficients - β 1 ( s ) and β 2 ( s ) - are 78 km and 743 km, respectively, while that of the spatial random effect, ω , is 94 km. Again, these estimates suggest much greater spatial heterogeneities in the differences in coverage between age groups 9-11 months and 12-23 months than between 9-11 months and 24-35 months. Nevertheless, the corresponding marginal variance estimates for these terms are very close, and the iid residual term, ϵ , accounts for nearly twice as much variation in the model as ω .
Similar patterns are observed in the results of MODsmooth reported in supplementary Table 4 and supplementary Figure 6, both in terms of significant covariate effects and the effect of age on vaccination. In particular, the delay in the receipt of MCV1 is evident in the estimated nonlinear effect of age on vaccination, which shows a markedly lower probability of vaccination among children aged < 15 months relative to older ages. We observe a steep increase in the probability of vaccination from age 9 months to about 15 months, which then stabilizes at a high level with minor fluctuations beyond 23 months, suggesting immunity gaps in the 24-35 month birth cohort. The credible intervals for the smooth functions indicate greater uncertainty beyond age 12 months due to greater fluctuations in coverage at older ages (see Figure 1a).

Predicted 1x1 km Maps of Vaccination Coverage for Age Groups and Single age Points

Geographically, considerable heterogeneities in MCV1 coverage are observed across all age groups (Figure 2), with a pronounced effect of walking travel time to the nearest health facility (see supplementary Figure 4). The 1x1 km maps generally reveal lower coverage in the 9-11 month age group relative to other age groups, particularly concentrated in the western and northern districts of the country. The spatial distributions of coverage in these age groups show some interesting similarities, but there are also striking differences. Across all the age groups, coverage levels are consistently higher in and around Lacs, Vallée du Bandama and Yamoussoukro districts. Whilst overall coverage levels do not differ substantially between ages 12-23 months and 24-35 months, the spatial patterns in coverage appear to be more similar between ages 9-11 months and 24-35 months than 12-23 months, with the latter showing more pronounced spatial gradients, as noted previously. The modelled estimates at the national level are 54.1%, 67.9% and 70.9% for the respective age groups. Compared to corresponding direct survey estimates of 51.2%, 60.3% and 65.1%, the model slightly overestimated coverage, especially in the 12-23 month age group. For the combined 9-35 month age group, similar spatial patterns are observed, although aggregating over this wider age range obscures differences in coverage between birth cohorts and is more suitable when interest is in estimating for broader age bands.
The widths of the 95% credible intervals (CIs) associated with the coverage maps indicate greater uncertainties in areas with sparse survey data (e.g., the Zanzan region; see Figure 1) and in areas of lower coverage, especially for age groups 12-23 months and 24-35 months. Uncertainties for the 9–11 month group show less spatial variation but are slightly higher overall, likely reflecting smaller sample sizes. The median widths of the 95% CIs are 50%, 48%, and 43% for the 9–11, 12–23, and 24–35 month groups, respectively, indicating considerable uncertainties in the grid-level estimates.
At single age points (Figure 3), lower coverage at 9 months is also evident, with very similar spatial patterns observed from ages 12 to 35 months. (Note that MODsmooth does not capture spatial gradients in the differences between age groups.) The modelled national coverage estimates of 47.6%, 65%, 68.3%, 70%, 70.3% and 68.7% at 9, 12, 15, 18, 24 and 35 months, respectively, further confirm that coverage patterns stabilize and become more similar from age 15 months and above. The corresponding uncertainties, reported in supplementary Figure 7, show similar patterns and ranges as those in Figure 2.

Estimation of Vaccination Coverage at the Administrative Level

At the departmental level, differences in both coverage and spatial prioritization for the 9–11, 12–23, and 24–35 month age groups are shown in Figure 4 and Supplementary Figure 8. Consistent with the national-level estimates reported earlier, coverage remains suboptimal at the departmental level, with no department reaching the target of 95% coverage in any age group. In particular, for the 9-11 month age group, coverage is estimated to be below 60% in about 70% of departments. For the 24–35 month age group, clusters of departments with lower coverage (<60%) are observed in the neighboring Savanes (M’Bengué, Tengrela, Boundiali, and Kouto) and Denguélé (Kaniasso, Madinani, and Séguélon) districts, despite overall higher coverage at older ages - potentially contributing to the dip in coverage after 23 months noted in supplementary Figure 6.
Among 9-11 month olds, the highest coverage levels (>67%) are estimated in M’Bahiakro, Daoukro, and Prikro (all in the Lacs district), while the lowest (<44%) are in Dianra, Daloa, and Séguélon departments. For 12-23 months, coverage is highest (>81%) in M’Bahiakro, Tanda, and Transua (Lacs and Zanzan districts), and lowest (<47%) in Daloa, Dianra, and Guitry. In the 24–35 month group, coverage peaks (>85%) in Tanda, Koun-Fao, M’Bahiakro, and Transua (Zanzan and Lacs districts), and is lowest in Séguélon, Kouto, and Dianra. Although there are similarities in the locations of higher and lower coverage across age groups, differences also emerge, as reflected in the departmental rankings shown in Figure 4b. For example, even within Lacs district - where coverage is consistently higher across all age groups - the ranks of departments differ when comparing the 12–23 month group with the 9–11 and 24–35 month groups. The rankings of the 9–11 and 24–35 month groups are more similar, corroborating earlier findings about the resemblance in their spatial coverage patterns. The median absolute differences in departmental ranks between age groups are: 9 (IQR=12) for 9–11 vs. 12–23 months, 8 (IQR=11) for 9–11 vs. 24–35 months, and 14 (IQR=18) for 12–23 vs. 24–35 months. Furthermore, the rankings appear to vary most in the Zanzan, Denguélé, and Montagnes districts (based on the larger IQRs of ranks across all age groups within these districts).
The uncertainties associated with the departmental-level estimates (supplementary Figure 9) are considerably smaller than those of the grid-level estimates, with median 95% CI widths of 23%, 22%, and 18% for the respective age groups.
We also estimated coverage at the district level for the 12–23 month age group to enable comparison with direct survey estimates at this level (supplementary Figure 10). Overall, there is good agreement between the modelled and direct survey estimates; however, the relative overestimation of coverage observed earlier is most pronounced in the Abidjan, Comoé, and Lagunes districts.

Discussion

Our study proposed and evaluated a spatially varying coefficient model for age-structured mapping of vaccination coverage, along with alternative specifications offering the flexibility to map coverage for both age groups and single age points. The methodology is also applicable to mapping other health and development indicators where understanding age-related differences is of interest. The proposed models can be readily implemented using the INLA-SPDE approach, although predictions at a high resolution (e.g., 1 × 1 km) and over large study areas may require high-memory computing environments.
Among the three competing birth cohort models, our validation results based on the 2021 CDHS data provided sufficient support for the spatially varying coefficient model, MODsvc1. We also found that explicitly modelling age improved predictive performance at the individual level. As expected, however, out-of-sample prediction across multiple age categories is constrained by smaller sample sizes compared to predictions for all ages combined. The impact of cluster-level sample sizes on out-of-sample prediction has been explored in detail elsewhere [19]. For larger age groups, more parsimonious parameterizations of MODsvc1 - incorporating fixed step changes between certain age groups - can be informed by the data, guided by model selection criteria, or refined through inspection of the predicted surfaces from the fully parameterized model and subsequent adjustments in later analyses.
The estimated covariate effects in our analysis highlight the need for interventions aiming to improve coverage in urban areas, conflict areas, areas with limited access to health facilities/vaccination services, as well as those targeted towards improved maternal/caregiver education, improved vaccination card retention (also an indication of improved access to vaccination services) and improved vaccine information communication in Cote d’Ivoire. These findings are corroborated by previous studies [44,45] which analysed similar indicators of vaccination coverage in Cote d’Ivoire. Our study also demonstrated a profound effect of age on vaccination. In the birth cohort model MODsvc1, we estimated that children aged 12–23 months and 24–35 months had approximately twice the odds of being vaccinated compared to those aged 9–11 months. We also observed substantial spatial heterogeneity in the differences in coverage between the 9–11 and 12–23 month age groups. In the single-age model MODsmooth, we found markedly lower vaccination probabilities among children younger than 15 months, as well as slight dips in coverage within the 24–35 month cohort. Together, these findings provide strong evidence of significant delays in MCV1 vaccination during the first year of life. If these patterns persist, interventions targeting eligible children during the first year of life and up to 15 months (e.g., improved maternal/caregiver education and improved access to vaccination services) are likely to be most impactful.
Our predicted maps clearly showed that MCV1 coverage is lowest at younger ages (around the first year of life) in the northern and western districts, and further corroborated the relatively higher and more uniform coverage at older ages. These spatial patterns in coverage are also consistent with Ahoussou et al’s [46] study, which analysed measles case-based surveillance data in Cote d’Ivoire and found relatively higher odds of detection of positive measles-specific immunoglobin M (IgM) in health regions in the west of the country. At the grid, departmental, and district levels, estimated coverage remained well below the 95% target across the country. Specifically, in the 9–11 month age group, coverage was estimated to be below 60% in about 70% of departments. In the 24–35 month age group, despite generally higher coverage, areas of lower coverage (<60%) were concentrated in the neighboring Savanes (M’Bengué, Tengrela, Boundiali, and Kouto) and Denguélé (Kaniasso, Madinani, and Séguélon) districts. Across all age groups, we identified notable differences in departmental rankings, highlighting distinct priority areas for each age group, although some consistent patterns also emerged. For example, the departments of Dianra (Woroba district), Daloa (Sassandra–Marahoué), and Séguélon (Denguélé) were consistently among the three lowest-coverage areas in at least two age groups, whereas Guitry and Kouto were among the lowest-coverage areas only for the 12–23 and 24–35 month groups, respectively. Departmental rankings were most variable in the Zanzan, Denguélé, and Montagnes districts, indicating greater heterogeneity in age-specific coverage within these regions. The uncertainties in these modelled estimates were higher at the grid level (with median 95% CI widths ranging between 43% and 50%) than at the departmental level (with median 95% CI widths ranging between 18% and 23%), suggesting that departmental-level estimates are more actionable for programmatic purposes. Nevertheless, producing estimates at the grid level first provides the flexibility to aggregate to operationally relevant spatial scales where precision improves, while retaining the spatial granularity necessary for targeted interventions.
Our results have several implications for measles control and elimination efforts in Cote d’Ivoire. WHO and UNICEF estimates indicate that national MCV1 coverage remained suboptimal at 70% in 2023 [47]. In 2024 alone, Cote d’Ivoire reported1,872 measles cases to WHO - a 37.5% increase from the previous year [48]. Addressing existing immunity gaps is therefore critical for progress towards measles control and elimination, and our study provides evidence to help guide intervention strategies. First, our modelled estimates for the 9–11 month age group can support routine immunization interventions aimed at improving coverage in the first year of life, or in children up to 15 months of age, by identifying and prioritizing the areas most at risk. Second, these estimates can inform the planning of catch-up vaccination activities, particularly for children currently aged 5–9 years, which encompasses two of the birth cohorts analyzed in our study (12–23 months and 24–35 months at the time of the survey). For spatial prioritization, the rankings of subnational areas, as demonstrated in Figure 4, could be enhanced by using estimates of the number of zero-dose children, which may better capture areas of highest disease risk. These zero-dose estimates can also feed into microplanning for vaccination activities, while field data collected during such activities can in turn be used to validate and refine the modelled estimates.
Our findings could be further strengthened through triangulation with measles case-based surveillance data or serosurveillance data that reflect the national immunity profile, helping to build a more robust evidence base for targeted interventions [25,46,49]. On the technical side, the spatially varying coefficient model developed here can be extended to age-dependent analyses of these aggregate programmatic data at the district or provincial level [23]. Such extensions will typically involve modelling the spatially-varying regression coefficients and residual spatial dependence in the model using conditional autoregressive priors [50] and considering alternative probability distributions for the outcome. Additionally, it is possible to adapt MODsmooth to produce estimates over a specified age range. This will typically involve computing 1 z 2 − z 1 ∫ z 1 z 2 f z d z by defining a fine grid within the interval [ z 1 ,    z 2 ] , computing the posterior mean of f ( z ) at each grid point and then averaging over these points. In R-INLA, the integration points over [ z 1 ,    z 2 ] can be included in f . during model fitting and then the marginal posterior distributions of f ( . ) can be obtained post model-fitting and averaged to compute the integral. However, as noted earlier, this approach does not account for spatial non-stationarity when estimating differences in coverage between the age groups. Moreover, estimating age effects at single age points using a spatially varying coefficient model is not well suited for capturing the overall smoothed effect of age (see Supplementary Figure 11). Indeed, in our 2021 CDHS application, such a specification underperformed (WAIC = 5232) relative to MODsmooth in equation (4). In future work, we plan to explore spatiotemporal extensions of the proposed models to allow age- and time-specific estimation of coverage. We anticipate that this analysis will yield additional programmatic insights into how the timeliness of vaccination and differences in coverage between the age groups have changed over space and time within and across countries.
Our analysis has some limitations. In the CDHS application, small cluster-level sample sizes - particularly in the 9–11 month age group - resulted in greater uncertainty in the predictions. We also observed that the fitted model slightly overestimated coverage in some southeastern areas of the country, potentially due to limited sample sizes, poor spatial coverage of clusters in these areas, or oversmoothing by the model. Furthermore, our analysis relied on information from both caregiver recall and vaccination records, which could introduce recall bias. The age range considered in our application was also constrained by data availability – DHS surveys only collect information on vaccination coverage for children aged 0-35 months at the time of the survey. As noted earlier, analyses covering broader age ranges are likely to yield greater insights for program planning and decision-making.
In summary, our study demonstrates the value of the proposed spatially varying coefficient model and alternative specifications for producing high-resolution maps of vaccination coverage by age and identifying priority areas for interventions. By highlighting substantial geographic and age-related heterogeneities in MCV1 coverage in our application, our findings underscore the need for targeted strategies to close immunity gaps, particularly in younger children and underserved areas. The proposed methodological framework is flexible, scalable, and applicable to other health indicators, offering a useful tool to support evidence-based decision-making for immunization programs and beyond. Future extensions to the spatiotemporal setting will further enhance its utility.

Supplementary Materials

The following supporting information can be downloaded at the website of this paper posted on Preprints.org.

Data and code availability

The data used in this study are available from the cited sources. The authors do not have the permission to redistribute these data. All R scripts used in the analyses are available on GitHub: https://github.com/EdsonUtazi/SVC_mod_paper_code.

Author Contributions

C. Edson Utazi: Conceptualization, Data Curation, Methodology, Formal analysis, Software, Supervision, Visualization, Funding acquisition, Writing – original draft, Writing – review and editing. Somnath Chaudhuri: Data Curation, Methodology, Formal analysis, Software, Visualization, Writing – review and editing. Oghenebrume Wariri: Methodology, Writing – review and editing. Iyanuloluwa D. Olowe: Data Curation, Visualization, Writing – review and editing. Mohamed Megheib: Formal analysis, Writing – review and editing. Andrew J. Tatem: Supervision, Funding acquisition, Writing – review and editing.

Funding

This study was supported by funding from UNICEF [grant number: 43387656] and Gavi, the Vaccine Alliance.

Acknowledgements

This work represents an output of the VaxPop project led by CEU. The authors would like to thank colleagues from the WorldPop program at the University of Southampton for their support during the study.

Conflicts of Interest

The authors declare no competing interests.

References

  1. World Health Organization. Immunization Agenda 2030: A global strategy to leave no one behind; 2020. Available online: https://www.who.int/immunization/immunization_agenda_2030/en/ (accessed on 25 June 2020).
  2. Cutts, F.T.; Claquin, P.; Danovaro-Holliday, M.C.; Rhoda, D.A. Monitoring vaccination coverage: Defining the role of surveys. Vaccine 2016, 34, 4103–9. [Google Scholar] [CrossRef] [PubMed]
  3. Utazi, C.E.; Thorley, J.; Alegana, V.A.; Ferrari, M.J.; Takahashi, S.; Metcalf, C.J.E.; et al. Mapping vaccination coverage to explore the effects of delivery mechanisms and inform vaccination strategies. Nature Communications 2019, 10, 1633. [Google Scholar] [CrossRef] [PubMed]
  4. Auzenbergs, M.; Fu, H.; Abbas, K.; Procter, S.R.; Cutts, F.T.; Jit, M. Health effects of routine measles vaccination and supplementary immunisation activities in 14 high-burden countries: a Dynamic Measles Immunization Calculation Engine (DynaMICE) modelling study. The Lancet Global Health 2023, 11, e1194–e204. [Google Scholar] [CrossRef] [PubMed]
  5. Nelson, K.N.; Wallace, A.S.; Sodha, S.V.; Daniels, D.; Dietz, V. Assessing strategies for increasing urban routine immunization coverage of childhood vaccines in low and middle-income countries: A systematic review of peer-reviewed literature. Vaccine 2016, 34, 5495–503. [Google Scholar] [CrossRef] [PubMed]
  6. Gavi The Vaccine Alliance. Gavi Strategy 5.0, 2021-2025; 2020. Available online: https://www.gavi.org/our-alliance/strategy/phase-5-2021-2025 (accessed on 25 June 2021).
  7. Haeuser, E.; Byrne, S.; Nguyen, J.; Raggi, C.; McLaughlin, S.A.; Bisignano, C.; et al. Global, regional, and national trends in routine childhood vaccination coverage from 1980 to 2023 with forecasts to 2030: a systematic analysis for the Global Burden of Disease Study 2023. The Lancet 2025. [Google Scholar] [CrossRef] [PubMed]
  8. Cutts, F.T.; Ferrari, M.J.; Krause, L.K.; Tatem, A.J.; Mosser, J.F. Vaccination strategies for measles control and elimination: time to strengthen local initiatives. BMC Medicine 2021, 19, 2. [Google Scholar] [CrossRef] [PubMed]
  9. Wariri, O.; Utazi, C.E.; Okomo, U.; Metcalf, C.J.E.; Sogur, M.; Fofana, S.; et al. Mapping the timeliness of routine childhood vaccination in The Gambia: A spatial modelling study. Vaccine 2023, 41, 5696–705. [Google Scholar] [CrossRef] [PubMed]
  10. Manaksha, M.; Danya Arif, S.; Vijay Kumar, D.; Mubarak Taighoon, S.; Sundus, I.; Hamidreza, S.; et al. Coverage, timeliness of measles immunisation and its predictors in Pakistan: an analysis of 6.2 million children enrolled in the Provincial Electronic Immunisation Registry. BMJ Global Health 2025, 10, e016717. [Google Scholar] [CrossRef] [PubMed]
  11. World Health Organization, United Nations Children’s Fund (UNICEF), Gavi the Vaccine Alliance. The Big Catch-Up: An Essential Immunization Recovery Plan for 2023 and Beyond; 2023. Available online: https://www.who.int/publications/i/item/9789240075511 (accessed on 12 July 2025).
  12. O'Brien, K.L.; Lemango, E. The big catch-up in immunisation coverage after the COVID-19 pandemic: progress and challenges to achieving equitable recovery. The Lancet 2023, 402, 510–2. [Google Scholar] [CrossRef] [PubMed]
  13. Nambiar, D.; Hosseinpoor, A.R.; Bergen, N.; Danovaro-Holliday, M.C.; Wallace, A.; Johnson, H.L. Inequality in Immunization: Holding on to Equity as We ‘Catch Up’. Vaccines 2023, 11, 913. [Google Scholar] [CrossRef] [PubMed]
  14. Sbarra, A.N.; Rolfe, S.; Nguyen, J.Q.; Earl, L.; Galles, N.C.; Marks, A.; et al. Mapping routine measles vaccination in low- and middle-income countries. Nature 2021, 589, 415–9. [Google Scholar]
  15. Takahashi, S.; Metcalf, C.J.E.; Ferrari, M.J.; Tatem, A.J.; Lessler, J. The geography of measles vaccination in the African Great Lakes region. Nature Communications 2017, 8, 15585. [Google Scholar] [CrossRef] [PubMed]
  16. Utazi, C.E.; Nilsen, K.; Pannell, O.; Dotse-Gborgbortsi, W.; Tatem, A.J. District-level estimation of vaccination coverage: Discrete vs continuous spatial models. Stat Med. 2021, 40, 2197–211. [Google Scholar] [CrossRef]
  17. Utazi, C.E.; Thorley, J.; Alegana, V.A.; Ferrari, M.J.; Takahashi, S.; Metcalf, C.J.E.; et al. High resolution age-structured mapping of childhood vaccination coverage in low and middle income countries. Vaccine 2018, 36, 1583–91. [Google Scholar] [CrossRef] [PubMed]
  18. Utazi, C.E.; Thorley, J.; Alegana, V.A.; Ferrari, M.J.; Nilsen, K.; Takahashi, S.; et al. A spatial regression model for the disaggregation of areal unit based data to high-resolution grids with application to vaccination coverage mapping. Statistical Methods in Medical Research 2018, 28, 3226–41. [Google Scholar] [CrossRef] [PubMed]
  19. Utazi, C.E.; Aheto, J.M.K.; Chan, H.M.T.; Tatem, A.J.; Sahu, S.K. Conditional probability and ratio-based approaches for mapping the coverage of multi-dose vaccines. Statistics in Medicine 2022, 41, 5662–78. [Google Scholar] [CrossRef] [PubMed]
  20. Utazi, C.; Yankey, O.; Chaudhuri, S.; Olowe, I.; Danovaro-Holliday, C.; Lazar, A.; et al. Geostatistical and Machine Learning Approaches for High-Resolution Mapping of Vaccination Coverage. Preprints 2024. [Google Scholar]
  21. Gelfand, A.E.; Kim, H.-J.; Sirmans, C.F.; Banerjee, S. Spatial Modeling With Spatially Varying Coefficient Processes. Journal of the American Statistical Association 2003, 98, 387–96. [Google Scholar] [CrossRef] [PubMed]
  22. Finley, A.O. Comparing spatially-varying coefficients models for analysis of ecological data with non-stationary and anisotropic residual dependence. Methods in Ecology and Evolution 2011, 2, 143–54. [Google Scholar] [CrossRef]
  23. Meehan, T.D.; Michel, N.L.; Rue, H. Spatial modeling of Audubon Christmas Bird Counts reveals fine-scale patterns and drivers of relative abundance trends. Ecosphere 2019, 10, e02707. [Google Scholar] [CrossRef]
  24. Fotheringham, A.S.; Brunsdon, C.; Charlton, M. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Chichester: Wiley; 2003.
  25. Utazi, C.E.; Aheto, J.M.K.; Wigley, A.; Tejedor-Garavito, N.; Bonnie, A.; Nnanatu, C.C.; et al. Mapping the distribution of zero-dose children to assess the performance of vaccine delivery strategies and their relationships with measles incidence in Nigeria. Vaccine 2023, 41, 170–81. [Google Scholar] [CrossRef] [PubMed]
  26. Institut National de la Statistique, I.C.F. Côte d’Ivoire Enquête Démographique et de Santé 2021 Rapport final. Rockville, Maryland, USA et la Côte d’Ivoire: INS et ICF; 2023.
  27. UNICEF Reaching the unreached with life-saving vaccines through data science and geospatial technologies; 2025. Available online: https://data.unicef.org/resources/reaching-the-unreached-with-life-saving-vaccines-through-data-science-and-geospatial-technologies/ (accessed on 12 July 2025).
  28. World Health Organisation. Vaccination schedule for Measles (Cote d'Ivoire); 2025. 2025. Available online: https://immunizationdata.who.int/global/wiise-detail-page/vaccination-schedule-for-measles?ISO_3_CODE=CIV&TARGETPOP_GENERAL=. (accessed on 12 July 2025).
  29. Cressie, N. Statistics for Spatial Data: Wiley; 2015.
  30. Tatem, A.J. WorldPop, open data for spatial demography. Scientific Data 2017, 4, 170004. [Google Scholar] [CrossRef] [PubMed]
  31. Runfola, D.; Anderson, A.; Baier, H.; Crittenden, M.; Dowker, E.; Fuhrig, S.; et al. geoBoundaries: A global database of political administrative boundaries. PLoS One 2020, 15, e0231866. [Google Scholar] [CrossRef] [PubMed]
  32. Matérn, B. Spatial Variation. 2nd ed. Berlin, Germany: Springer-Verlag; 1960.
  33. Lindgren, F.; Rue, H.; Lindström, J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J Roy Stat Soc Series B (Stat Methodol) 2011, 73, 423–98. [Google Scholar] [CrossRef]
  34. Wang, X.; Yue, Y.R.; Faraway, J.J. Bayesian regression modeling with INLA: Chapman and Hall/CRC; 2018.
  35. Browne, W.J.; Subramanian, S.V.; Jones, K.; Goldstein, H. Variance partitioning in multilevel logistic models that exhibit overdispersion. Journal of the Royal Statistical Society: Series A (Statistics in Society) 2005, 168, 599–613. [Google Scholar] [CrossRef]
  36. Simpson, D.; Rue, H.; Riebler, A.; Martins, T.G.; Sørbye, S.H. Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science 2017, 32, 1–28. [Google Scholar] [CrossRef]
  37. Fuglstad, G.-A.; Simpson, D.; Lindgren, F.; Rue, H. Constructing priors that penalize the complexity of Gaussian random fields. Journal of the American Statistical Association 2019, 114, 445–52. [Google Scholar] [CrossRef]
  38. Rue, H.; Martino, S.; Chopin, N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2009, 71, 319–92. [Google Scholar] [CrossRef]
  39. Lindgren, F.; Bachl, F.; Illian, J.; Suen, M.H.; Rue, H.; Seaton, A.E. inlabru: software for fitting latent Gaussian models with non-linear predictors. arXiv preprint, 2024; arXiv:240700791. [Google Scholar]
  40. Watanabe, S. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of machine learning research 2010, 11, 3571–94. [Google Scholar]
  41. Gelman, A.; Hwang, J.; Vehtari, A. Understanding predictive information criteria for Bayesian models. Statistics and Computing 2014, 24, 997–1016. [Google Scholar] [CrossRef]
  42. Gneiting, T.; Raftery, A.E. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association 2007, 102, 359–78. [Google Scholar] [CrossRef]
  43. Brier, G.W. Verification of forecasts expressed in terms of probability. Monthly Weather Review 1950, 78, 1–3. [Google Scholar] [CrossRef]
  44. Ogundele, O.A.; Ogunwemimo, H.S.; Fehintola, F.O.; Ogundele, T.; Olorunsola, A.; Bello, O.E.; et al. Predictors of incomplete childhood vaccination in four West African countries: a population based cross-sectional study. Sci Rep. 2025, 15, 17119. [Google Scholar] [CrossRef] [PubMed]
  45. Douba, A.; Aka Lepri Bernadin, N.; Attoh-Toure, H.; Bangaman Akani, C.; Yao, G.H.A.; Konan Ng et, a.l. An analysis of risk factors for incomplete immunization for children in Côte d’Ivoire: Examination of 1998-1999 and 2011-2012 Demographic and Health Survey. HEALTH SCIENCES AND DISEASE 2016, 17. [Google Scholar]
  46. Ahoussou, E.M.K.; Ekra, K.D.; Dananche, C.; Aka, L.B.N.; Rabilloud, M.; Vanhems, P. Factors associated with detection of measles-specific immunoglobulin M (IgM-positive) in Côte d'Ivoire from 2015 to 2019: A nested case-control study from a measles national surveillance programme. Journal of Public Health and Epidemiology 2023, 15, 97–105. [Google Scholar]
  47. World Health Organization. WHO/UNICEF Estimates of National Immunization Coverage (WUENIC) 2024 revision; 2024. Available online: https://www.who.int/teams/immunization-vaccines-and-biologicals/immunization-analysis-and-insights/global-monitoring/immunization-coverage/who-unicef-estimates-of-national-immunization-coverage (accessed on 5 July 2025).
  48. World Health Organisation Côte d'Ivoire Reported cases of vaccine-preventable diseases, (.V.P.D.s.).; 2025 Available online: http://www.w3.org/1999/xlink" xlink:href="h.t.t.p.s.:././.i.m.m.u.n.i.z.a.t.i.o.n.d.a.t.a.w.h.o.i.n.t./.d.a.s.h.b.o.a.r.d./.r.e.g.i.o.n.s./.a.f.r.i.c.a.n.-r.e.g.i.o.n./.C.I.V (accessed on 12 July 2025).
  49. Carcelen, A.C.; Winter, A.K.; Moss, W.J.; Chilumba, I.; Mutale, I.; Chongwe, G.; et al. Leveraging a national biorepository in Zambia to assess measles and rubella immunity gaps across age and space. Scientific Reports 2022, 12, 10217. [Google Scholar] [CrossRef] [PubMed]
  50. Lee, D. A comparison of conditional autoregressive models used in Bayesian disease mapping. Spat Spatiotemporal Epidemiol. 2011, 2, 79–89. [Google Scholar] [CrossRef] [PubMed]
Preprints 168458 g001
Figure 1. Empirical estimates of measles vaccination coverage for children aged 9 to 35 months at the national (a, b) and survey cluster (c) levels for (a) single age points and (b, c) age groups, according to the 2021 Cote d’Ivoire Demographic and Health Survey.
Figure 1. Empirical estimates of measles vaccination coverage for children aged 9 to 35 months at the national (a, b) and survey cluster (c) levels for (a) single age points and (b, c) age groups, according to the 2021 Cote d’Ivoire Demographic and Health Survey.
Preprints 168458 g001
Preprints 168458 g002
Figure 2. 1x1 km maps of MCV1 coverage and associated uncertainties for age groups 9-11 months, 12-23 months, 24-35 months and 9-35 months.
Figure 2. 1x1 km maps of MCV1 coverage and associated uncertainties for age groups 9-11 months, 12-23 months, 24-35 months and 9-35 months.
Preprints 168458 g002
Preprints 168458 g003
Figure 3. 1x1 km maps of MCV1 coverage for single age points at 9, 12, 15, 18, 24 and 35 months. Corresponding uncertainty estimates, obtained as the widths of the 95% credible intervals, are shown in the supplementary file.
Figure 3. 1x1 km maps of MCV1 coverage for single age points at 9, 12, 15, 18, 24 and 35 months. Corresponding uncertainty estimates, obtained as the widths of the 95% credible intervals, are shown in the supplementary file.
Preprints 168458 g003
Preprints 168458 g004
Figure 4. (a) Estimates of MCV1 coverage and (b) corresponding ranks at the departmental level for each age group. The vertical lines differentiate between departments in each of the 14 districts.
Figure 4. (a) Estimates of MCV1 coverage and (b) corresponding ranks at the departmental level for each age group. The vertical lines differentiate between departments in each of the 14 districts.
Preprints 168458 g004
Table 1. WAIC values of the proposed models.
Table 1. WAIC values of the proposed models.
Model p W A I C WAIC
MODsvc1 210 5191
MODsvc2 219 5197
MODnosvc 186 5215
MODsmooth 187 5204
MODall 184 5262
Table 2. Results of cluster level k-fold cross-validation.
Table 2. Results of cluster level k-fold cross-validation.
Model Age group RMSE MAE AVG_BIAS CRPS
RANDOM
MODsvc1 9 -11 0.439 0.403 0.006 0.332
MODsvc2 9 -11 0.440 0.404 0.006 0.335
MODnosvc 9 -11 0.440 0.403 0.008 0.334
MODsvc1 12 - 23 0.315 0.263 0.008 0.206
MODsvc2 12 - 23 0.315 0.262 0.009 0.203
MODnosvc 12 - 23 0.314 0.262 0.012 0.208
MODsvc1 24 - 35 0.298 0.249 -0.015 0.197
MODsvc2 24 - 35 0.299 0.249 -0.015 0.196
MODnosvc 24 - 35 0.300 0.251 -0.016 0.198
STRATIFIED
MODsvc1 9 - 11 0.444 0.408 0.004 0.337
MODsvc2 9 - 11 0.445 0.408 0.004 0.341
MODnosvc 9 - 11 0.445 0.409 0.004 0.338
MODsvc1 12 - 23 0.318 0.266 0.003 0.206
MODsvc2 12 - 23 0.320 0.267 0.002 0.207
MODnosvc 12 - 23 0.319 0.265 0.009 0.209
MODsvc1 24 - 35 0.304 0.252 -0.017 0.196
MODsvc2 24 - 35 0.302 0.253 -0.017 0.198
MODnosvc 24 - 35 0.304 0.255 -0.018 0.200
Table 3. Table of parameter estimates. Estimates corresponding to significant regression coefficients are bolded.
Table 3. Table of parameter estimates. Estimates corresponding to significant regression coefficients are bolded.
Parameter Mean Odds ratio Std. Dev. 2.5% 97.5%
β ^ 0 -3.835 0.022 5.557 -14.727 7.057
Urban -0.611 0.543 0.143 -0.892 -0.33
Veg_index -3.384 0.034 2.143 -7.585 0.816
Wetdays -0.102 0.903 0.102 -0.303 0.099
Dist_conf 0.214 1.239 0.091 0.035 0.393
Elevation − 1 × 10 − 5 1.000 1 × 10 − 5 − 3 × 10 − 5 1 × 10 − 5
Urban_access − 2 × 10 − 4 1.000 5 × 10 − 4 -0.001 0.001
Walking_tt -0.002 0.998 0.001 -0.004 -0.001
Mal_prev -0.558 0.572 0.834 -2.192 1.077
Max_temp 0.117 1.124 0.144 -0.166 0.4
Mat_educ 1.035 2.815 0.236 0.572 1.498
Health_card 1.824 6.197 0.302 1.231 2.417
Media 0.552 1.737 0.272 0.020 1.085
Wealth -0.135 0.874 0.247 -0.619 0.349
β ^ 1 ( 12 − 23   m o n t h s ) 0.716 2.046 0.124 0.473 0.959
β ^ 2 ( 24 − 35   m o n t h s ) 0.827 2.286 0.305 0.228 1.425
r ^ 0.844 - 0.456 0.293 2.028
σ ^ ω 0.367 - 0.094 0.216 0.583
r ^ β 1 0.702 - 0.344 0.283 1.595
σ ^ β 1 0.399 - 0.112 0.214 0.649
r ^ β 2 6.697 - 5.659 1.561 21.783
σ ^ β 2 0.400 - 0.17 0.171 0.828
σ ^ ϵ − 2 4.056 - 1.228 2.151 6.94
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
bsky
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2025 MDPI (Basel, Switzerland) unless otherwise stated