key: cord-0440038-s9pomkda
authors: Dong, Tracy Qi; Wakefield, Jon
title: Space-time smoothing models for sub-national measles routine immunization coverage estimation with complex survey data
date: 2020-07-07
journal: nan
DOI: nan
sha: 7959566600c953321694abfb46e3c9d0e4ff2b74
doc_id: 440038
cord_uid: s9pomkda

Despite substantial advances in global measles vaccination, measles disease burden remains high in many low- and middle-income countries. A key public health strategy for controling measles in such high-burden settings is to conduct supplementary immunization activities (SIAs) in the form of mass vaccination campaigns, in addition to delivering scheduled vaccination through routine immunization (RI) programs. To achieve balanced implementations of RI and SIAs, robust measurement of sub-national RI-specific coverage is crucial. In this paper, we develop a space-time smoothing model for estimating RI-specific coverage of the first dose of measles-containing-vaccines (MCV1) at sub-national level using complex survey data. The application that motivated this work is estimation of the RI-specific MCV1 coverage in Nigeria's 36 states and the Federal Capital Territory. Data come from four Demographic and Health Surveys, three Multiple Indicator Cluster Surveys, and two National Nutrition and Health Surveys conducted in Nigeria between 2003 and 2018. Our method incorporates information from the SIA calendar published by the World Health Organization and accounts for the impact of SIAs on the overall MCV1 coverage, as measured by cross-sectional surveys. The model can be used to analyze data from multiple surveys with different data collection schemes and construct coverage estimates with uncertainty that reflects the various sampling designs. Implementation of our method can be done efficiently using integrated nested Laplace approximation (INLA).

1. Introduction. Measles is a highly contagious, vaccine-preventable disease. Despite tremendous advances in the global measles vaccination coverage in the past few decades, challenges remain in many low-and middle-income countries (LMICs), where disease burdens are the highest (World Health Organization, 2020a). A key public health strategy in such high-burden settings is the combination of routine immunization (RI) and supplementary immunization activities (SIAs) (The Measles and Rubella Initiative, 2012) . Although the World Health Organization (WHO) recommends two doses of measles-containing-vaccine (MCV) for all children -first at the age of 9 months and second at the age of 15 months, most LMICs only administer one dose in their RI schedule (World Health Organization, 2017) . In order to increase measles vaccination coverage, SIAs are carried out in the form of mass immunization campaigns every 2-4 years. During these campaigns, health workers run fixed-post vaccination sites and provide MCV to all children in a pre-specified target age group, regardless of whether they have been vaccinated previously. For large countries where there is substantial heterogeneity in disease burden and RI coverage across geographical regions, SIA campaigns are often designed and implemented based on each region's specific situations. For example, the SIA campaigns in Nigeria's 20 northern states have been implemented at a different schedule, and sometimes with a different target age group, compared to the SIAs campaigns in its 17 southern states, as shown in Table S2 of the Supplementary Material .

To achieve balanced implementations of RI and SIAs against measles, robust measurement of sub-national RI-specific coverage is crucial. Program officers rely on accurate estimation of local RI coverage to identify areas where the routine vaccine delivery systems are the weakest (Biellik and Orenstein, 2018; The Measles and Rubella Initiative, 2019) or the SIA campaigns are the most effective (Utazi et al., 2020) . In addition, RI-specific coverage of the first dose of MCV (MCV1) is an essential input to measles epidemiological models for estimating underlying susceptible population dynamics and informing optimal SIA schedules (Verguet et al., 2015; Thakkar et al., 2019) . Specifically, in order to calculate the rate at which the susceptible population grows over time, one needs to estimate the number of infants who are not protected by vaccination after losing their maternal immunity against measles at the age of approximately 9 months. This requires reliable estimation of the RI-specific MCV1 coverage across different birth cohorts and sub-national areas.

Despite ongoing efforts to strengthen administrative data quality (BID Initiative, 2016) , the health management information systems in many LMICs are too weak to provide data of adequate quality for assessing and guiding health programs (Hancioglu and Arnold, 2013) . Therefore, household surveys are the primary data sources for vaccination coverage estimates (Cutts et al., 2016; Dolan and MacNeil, 2017) . The Demographic and Health Surveys (DHS) (ICF, 2020) and the Multiple Indicator Cluster Surveys (MICS) (UNICEF, 2020) are two major survey programs that report sub-national measles vaccination coverage estimates. In addition, geospatial modeling approaches have been applied to the household survey data to generate MCV coverage estimates at various geographical scales (Utazi et al., 2018 (Utazi et al., , 2019a . However, all the above-mentioned estimates are not RI-specific -they reflect the overall coverage of measles vaccination received via RI, SIAs or both sources. These coverage estimates can be influenced by past vaccination campaigns, and hence are not suitable for the aforementioned RI-specific purposes. For example, using the overall MCV coverage in epidemiological models may lead to under-estimation of the rate at which the susceptible population grows, resulting in inaccurate estimation of measles transmission dynamics. To obtain RI-specific coverage estimates, the current method relies on only using data from birth cohorts that have not had any SIA opportunities yet. This restriction may result in some or all data from a household survey not being eligible for analysis -an inefficient use of already sparse data.

In this paper, we propose a space-time smoothing model for estimating sub-national RIspecific MCV1 coverage using data from complex surveys. We demonstrate the method by calculating the state-level RI-specific MCV1 coverage estimates in Nigeria using data from 8 household surveys conducted between 2003 and 2018. Inspired by the small area estimation methods developed by Mercer et al. (2015) , our model is able to combine data from multiple surveys with different sampling designs and construct sub-national coverage estimates with uncertainty that reflects the various data collection schemes. Our method allows more efficient use of data by modeling survey outcome not just from children who only had RI, but also from those who had RI and an SIA opportunity. We account for the impact of SIAs by incorporating information from the WHO SIA calendar (World Health Organization, 2020b). This paper is structured as follows. In Section 2, we describe the Nigeria data upon which estimation will be based. The details of our proposed method are provided in Section 3. The results of our modeling efforts of RI-specific MCV1 coverage in Nigeria from 2003 to 2018 are presented in Section 4. Section 5 contains concluding remarks.

2. Data sources. We focus on the state-level RI-specific MCV1 coverage in Nigeria using data from four DHS surveys conducted in 2003 , 2008 , 2013 , three MICS surveys conducted in 2007 , 2011 , and two National Nutrition and Health Surveys (NNHS) conducted in 2014 and 2015. In addition, we obtain the SIA calendar from the WHO website, which contains the start and end dates, target age groups and geographic areas of all SIAs implemented in Nigeria. More details about the data are provided below.

2.1. Household surveys. The DHS, MICS and NNHS used similar, but not identical stratified cluster sampling designs to provide vaccination coverage estimates at the national level and in each of Nigeria's 36 states and the Federal Capital Territory (FCT) (Dong, Rhoda and Mercer, 2020) . Hereafter we refer to these as Nigeria's 37 states. The DHS surveys were stratified by urban and rural areas within each state, whereas MICS and NNHS were stratified only by state. Selection of survey respondents took place in three stages: First, within each stratum, primary sampling units (PSUs or clusters) were selected from a census frame of enumeration areas (EAs). The 2003 DHS and 2007 MICS used frames based on the 1991 census and the other surveys used frames based on the 2006 census. Second, within each selected cluster, a pre-specified number of households were randomly selected from an updated listing of all residential households. Third, within each selected household, field teams attempted to interview all women aged 15-49 years and collect vaccination data for children under age 3 or 5 years.

Every child in the survey is associated with a survey weight that quantifies the relative number of children he or she represents in the population. Generally, DHS and MICS surveys calculate a base weight for each child as the inverse of the product of the probabilities of selection from each sampling stage. The base weights are then adjusted to reflect household or individual response rates. In contrast, the NNHS were conducted based on the Standardized Monitoring and Assessment of Relief and Transitions (SMART) method (SMART Technical Advisory Group, 2006 Group, , 2017 : each child was assigned the same base weight across all clusters, and the base weights were then post-stratified to match a set of state population proportions without further adjustment for the non-response rate. Table S1 of the Supplementary Material summarizes key information regarding the sampling design of each survey, including source of the sampling frame, survey strata, cluster and household sample sizes, age group with vaccination data, and features of weight calculation.

For each child sampled in a survey, the care giver was asked to provide information about specific vaccines that the child had received up until the time of interview. Specifically, care givers were first asked to present any home-based record (HBR), such as vaccination cards, to the survey interviewer. If HBR is available, the interviewer would populate a designated section of the survey questionnaire that detail the types and doses of vaccines received by the child. After this step is complete, the interviewer would ask the care giver a follow-up question of whether additional vaccine doses had been administered besides the ones listed on the HBR. This would include vaccines received at health care facilities via RI but which had not been recorded on paper, or those received via SIAs. If additional doses were administered, the interviewer will go back to the same section and record any additional doses based on the care giver's recall. If no HBR is available, the interviewer would skip the HBR section and ask the care giver to recall whether the child has received each specific vaccine at the time of survey interview, regardless of whether the vaccination is administered through RI or SIAs. As such, the surveys were designed to measure the percentage of children in various birth cohorts who received specific vaccines at any time before the survey through any source, and it is extremely difficult, if not impossible, to differentiate whether a child received MCV1 through RI or SIAs based on the answers to the survey questionnaires alone.

2.2. SIA calendar. The WHO SIA calendar is a publicly available data source where key information on the completed and planned SIAs in member states is provided (World Health Organization, 2020b). For each SIA, the calendar records the start and end dates, target age group, size of the target population and geographical areas in which the campaign is carried out. Table S2 of the Supplementary Material shows the details of the SIAs in Nigeria conducted before the end of 2018 that we consider for this analysis. Excluded from the analysis are two mop-up campaigns implemented in 2007 and one outbreak response campaign implemented in 2013. These campaigns' scales were very small in terms of target population and geographic area, and therefore were assumed to have negligible impact on the state-level MCV1 coverage.

3. Methods. The basic idea of our method is that we consider the RI-specific MCV1 coverage as the percentage of children in a birth cohort who received a first dose of measles vaccine through the routine vaccine delivery system according to the RI schedule, i.e., at the age of 9 months. Each cross-sectional survey serves as a snap-shot measurement of the MCV1 coverage within a birth cohort at the time of survey interview. We can carefully align the birth time, RI schedule, SIA schedule and survey time for each child to determine whether a birth cohort has had the RI opportunity or any additional SIA opportunity when the survey data are collected. Figure 1 illustrates the data generating mechanism we assume for a generic area. We discretize time into 6-month intervals and let the beginning of 2000 be the reference starting point. We let b index birth cohort, such that the children born between Jan. and Jun. in 2000 belong to birth cohort b = 1, and the children born between Jul. and Dec. in 2000 belong to birth cohort b = 2, and so on. We let t index survey time, such that a survey conducted between Jan. and Jun. in 2001 has t = 3, and a survey conducted between Jan. and Jun. in 2002 has t = 5, and so on. We assume each child is scheduled to receive MCV1 through RI at the age of 9 months.

As an example, suppose an SIA targeting children aged 9-59 months is implemented between Jun. and Dec. 2001, and consider two birth cohorts in detail -cohort b = 1 represented by the blue circles, and cohort b = 2 represented by the orange circles. Each circle indicates the underlying MCV1 coverage of the birth cohort at a specific time. The blue/orange shade represents the proportion of children covered by RI, i.e., the RI-specific MCV1 coverage; and the yellow shade represent the proportion covered by SIA. Because all children within the target age group can receive vaccination during an SIA regardless of whether they have been vaccinated before, some children in the birth cohorts might be covered by both RI and SIA. They are represented by the intersection of the blue/orange shade and the yellow shade. Now we consider three annual cross-sectional surveys conducted at time t = 3, 5 and 7 respectively. Let p b,t denote the proportion of children in birth cohort b who have received MCV1 by time t. At t = 3, children in cohort b = 1 should already have had an RI opportunity to receive MCV1, but those in cohort b = 2 are still too young to be vaccinated against measles. Therefore, we can use the survey data to obtain a design-based weighted estimate of the MCV1 coverage for birth cohort b = 1 at survey time t = 3, which we denotep 1,3 . In Figure 1 ,p 1,3 and its associated uncertainty interval are represented by a blue dot with a vertical line segment. Because the birth cohort has not been impacted by any SIA at the time of survey,p 1,3 gives an estimate of the RI-specific MCV1 coverage. In the next survey at t = 5, both birth cohorts should have had their RI opportunities. In addition, both birth cohorts should also have been targeted during the SIA implemented in the previous time interval. Since survey data cannot differentiate the source of the vaccination, the designbased estimates of the MCV1 coverage for birth cohort b = 1 and 2 at time t = 5, denoted p 1,5 andp 2,5 , would reflect the overall MCV1 coverage in these birth cohorts instead of the RI-specific coverage. In Figure 1 , the overall MCV1 coverage is represented by the union of the blue/orange shade and the yellow shade. The same logic applies to the survey conducted at t = 7.

FIG 1. Illustration of the data generating mechanism for a generic area. We let b index 6-month birth cohort and consider two birth cohorts -cohort b = 1 represented by the blue circles, and cohort b = 2 represented by the orange circles. Each circle shows the underlying MCV1 coverage of the birth cohort at a specific time. The blue/orange shade represents the proportion of children covered by RI, i.e., the RI-specific MCV1 coverage; and the yellow shade represents the proportion covered by SIA. Because all children within the target age group can receive vaccination during an SIA, regardless of whether they have been vaccinated before, some children in the birth cohorts might be covered by both RI and SIA. They are represented by the intersection of the blue/orange shade and the yellow shade.

As mentioned in Section 1, our method uses survey data from children who only had RI and children who had RI and one SIA opportunity. Based on the data generating mechanism described above, our proposed approach takes three steps to construct sub-national RI-specific MCV1 coverage over time: first, individual-level data from each survey is processed to identify eligible birth cohorts; second, design-based MCV1 coverage estimates and associated variances are calculated for all eligible birth cohorts identified in step one at the sub-national area level; last, a space-time smoothing model is applied to the design-based estimates from step two to predict RI-specific MCV1 coverage over time for all sub-national areas. We now describe each step in detail.

3.1. Data processing. In each survey, the time of interview for each child is recorded in Century-Month Code (CMC) format. In addition, the child's age in months at the time of interview is also recorded. Combining the two pieces of information allows us to identify the birth month, and hence the birth cohort, for each child.

Assuming every child has their RI opportunity at the age of 9 months, we use the age information of each child to see whether he/she has had the RI opportunity at the time of the survey. We also align the SIA calendar with each child's birth month to identify how many SIA opportunities he/she has had before the survey interview.

With each child's birth cohort, RI and SIA information collected, we now identify, within each sub-national area, the birth cohorts in which all children only had their RI opportunity, or all children had RI plus one SIA opportunity in each survey. Data from these eligible birth cohorts will be used in the next two steps.

Design-based overall MCV1 coverage estimation. For each eligible birth cohort identified in each survey, the design-based overall MCV1 coverage estimate is calculated at the sub-national area level. Let i, b, s and k be the indices for area, birth cohort, survey and sampled child respectively. We obtain the Horvitz-Thompson (HT) direct estimate (Horvitz and Thompson, 1952) of the overall MCV1 coverage in birth cohort b in area i from survey s, denotedp ibs , usingp

where y ibsk is the 0/1 indicator of whether child k in birth cohort b in area i has ever received MCV1 when interviewed for survey s, w ibsk is the survey weight associated to the child, and n ibs is the total number of children in birth cohort b in area i who are sampled in survey s. In addition, we also obtain the design-based estimate of the variance associated withp ibs , which we denotev ibs . All design-based estimates can be computed using functions from the survey package (Lumley, 2004) within the R computing environment (Team, R Core, 2020).

Space-time smoothing model for RI-specific MCV1 coverage estimation. We extend the Fay-Herriot approach (Fay III and Herriot, 1979) developed in Mercer et al. (2015) and specify a Bayesian hierarchical space-time smoothing model to construct sub-national RIspecific MCV1 coverage estimates over time with uncertainty that reflects the various survey sampling designs.

First, the design-based estimates for the logit-transformed overall coverage and its associated variance, denoted byθ ibs andV ibs , are obtained viâ

.

Next, we take the following asymptotic distribution as a working likelihood:

where θ ibs is the underlying logit-transformed overall coverage in cohort b in area i measured by survey s. To account for the potential impact of SIA and specific surveys on the overall coverage, we specify θ ibs as

The potential impact of SIA on the overall MCV1 coverage is captured by β 1 × x ibs , where x ibs is an 0/1 indicator of whether children in birth cohort b in area i has had one SIA opportunity in addition to RI in survey s, and β 1 is a parameter quantifying the impact. We would expect β 1 > 0. In addition, we use the survey random effects s ∼ Normal (0, σ ) to capture any survey-specific biases, although this is relative to the average of all the surveys, and does not correct for any overall bias in the surveys combined.

The underlying RI-specific MCV1 coverage in birth cohort b in area i is captured by the RI term µ ib . Note that this quantity remains unchanged across surveys. Since we consider the RIspecific MCV1 coverage to be the percentage of children in a birth cohort who received a first dose of MCV through the routine immunization system at the age of 9 months, the underlying RI-specific coverage of a birth cohort only depends on when and where the children are born.

As such, if a birth cohort b in area i only had RI by the time of a survey s, x ibs would be zero and the underlying overall coverage p ibs measured in the survey would be

If the birth cohort has had one SIA opportunity in addition to RI by the time of the survey, the underlying overall coverage p ibs measured in the survey would be

Note that ecological fallacy is not a concern here, because x ibs is constant within areas. Since we only use data from birth cohorts in which all children only had their RI opportunity (x ibs = 0), or all children had RI plus one SIA opportunity (x ibs = 1) in a survey, there is no individual-level heterogeneity in x ibs within each birth cohort.

Space-time smoothing is achieved by modeling the RI term µ ib as a combination of six components, see Equation (3). Specifically, β 0 is a common intercept for all birth cohorts and areas. There are two spatial terms that correspond to the convolution model of Besag, York and Mollié (1991) . In particular, α i ∼ ICAR(σ α ) are the intrinsic conditional autoregressive (ICAR) terms (described in Section S1.1 of the Supplementary Material) for i = 1, . . . , I sub-national areas (e.g., I = 37 for Nigeria's 37 states), and γ i ∼ iid Normal(0, σ γ ) are independent and identically distributed (IID) random effects. There are also two temporal terms, with δ b ∼ RW2(σ δ ) following an (intrinsic) random walk of order 2 model (described in Section S1.2 of the Supplementary Material) to pick up structured temporal trends in RI-specific coverage, and τ b ∼ iid Normal(0, σ τ ) being IID to pick up short-term fluctuations with no structure. Finally, the space-time interaction term φ ib can be specified as one of the four types of interactions described in Knorr-Held (2000) to capture any additional independent or structured variation in RI coverage across space over time.

To construct RI-specific MCV1 coverage estimates and associated uncertainty intervals for each sub-national area, we use the Bayesian posterior samples of the RI component of the space-time smoothing model. Specifically, the percentage of children in birth cohort b in area i who received a first dose of MCV via RI, p RI,ib , can be estimated by drawing posterior samples

where the superscript (m) denotes the mth posterior sample of the respective component. The survey specific effects s are not included because they are bias terms.

Computation in this step, including model fitting and prediction, can be efficiently carried out using the Integrated Nested Laplace Approximation (INLA) method (Rue, Martino and Chopin, 2009 ) as implemented in the INLA package in R. In Section S4 of the Supplementary Material, we present a simulation study that demonstrates how our approach performs under a number of scenarios with surveys of various sizes. The simulation results also show that the proposed model improves estimation in terms of accuracy and precision through more efficient use of data under each scenario.

4. Applying methods to household surveys in Nigeria. We apply our proposed method to estimate the state-level RI-specific MCV1 coverage in Nigeria between 2000 and 2018 using the household survey and SIA calendar data described in Section 2. We discretize time into 6-month intervals and let the beginning of 2000 be the reference starting point. During this time period, all children in Nigeria were scheduled to receive one dose of MCV via RI at the target age of 9 months. 4.1. Data processing and design-based overall MCV1 coverage estimation. We carry out the data processing step outlined in Section 3.1 to identify, within each state, the birth cohorts in which all children only had their RI opportunity, or all children had RI plus one SIA opportunity in each survey. The design-based HT estimates of the state-level overall MCV1 coverage are calculated for each eligible birth cohort in each survey using the survey package (Lumley, 2004) . Figure 2 shows the resultant estimates with the associated confidence intervals across birth cohorts for Nigeria's 37 states. The plots are arranged to approximately match the relative geographical locations of the states. The different colors represent the surveys upon which the design-based estimates are calculated, and the shape of the dots indicates whether a birth cohort had RI Only or had RI + 1 SIA at the time of survey. The northern states generally have lower MCV1 coverage estimates than the southern states. As expected, within the same survey, the birth cohorts with an additional SIA opportunity tend to have higher estimated MCV1 coverage than the RI-Only cohorts, although the differences are often small in magnitude and difficult to determine here because of the large uncertainty and lack of a modeling framework. The coverage estimates of the same birth cohort from different surveys sometimes show large differences even after accounting for vaccination activity, suggesting that there might be some survey-specific effects. As a result of small sample sizes in some areas and surveys, the design-based coverage estimates for some birth cohorts are associated with very wide confidence intervals. 4.2. RI-specific MCV1 coverage estimation via space-time smoothing. We fit six spacetime smoothing models to the design-based estimates. Each model has the same working like-TABLE 1 The space-time interaction assumptions of the models.

Space-time Interaction φ ib IID-IID IID spatial effects × IID temporal effects ICAR-IID ICAR spatial effects × IID temporal effects IID-RW2 IID spatial effects × RW2 temporal effects ICAR-RW2 ICAR spatial effects × RW2 temporal effects IID-AR1 IID spatial effects × AR1 temporal effects ICAR-AR1 ICAR spatial effects × AR1 temporal effects lihood and model components as specified in Equations (1) -(3), but a different specification of the space-time interaction term φ ib . Specifically, we consider the IID-IID, ICAR-IID, IID-RW2, ICAR-RW2, IID-AR1, ICAR-AR1 models, whose space-time interaction assumptions are detailed in Table 1 (Knorr-Held, 2000) . In particular, AR1 stands for the autoregressive (AR) model of order 1 for Gaussian random vectors. The definition of the AR1 model can found in Section S1.3 of the Supplementary Material. Implementation of the models was carried out in R using the INLA package. We used the Penalized Complexity (PC) priors of Simpson et al. (2017) for all hyper-parameters. Specifically, we applied weak PC priors on the precision parameters for the random effects, such that Pr(σ > 2) = 0.01 for the spatial and temporal terms, and Pr(σ > 1) = 0.01 for the independent terms.

In Table 2 , we report the posterior medians and 95% CIs of the model parameters. In addition, we also compute three model assessment metrics: the deviance information criterion (DIC), the Watanabe-Akaike information criterion (WAIC) and the log-score conditional predictive ordinate (LCPO) values, all of which are described in detail in Section S2 of the Supplementary Material. Based on the estimates for β 1 , the odds of receiving at least one dose of MCV is 33.6% -37.7% higher after an SIA campaign is conducted in a birth cohort. The ICAR-AR1 Model performs the best among the six models based on all three model metrics. Table 3 contains the proportions of variance explained by each of the random effects in the space-time smoothing model. The "ICAR" column corresponds to the structured spatial effects; "Space IID" to the unstructured spatial effects; "RW2" to the structured temporal trends; "Time IID" to the unstructured temporal effects; and "Space×Time" to the respective space-time interaction terms in each model. Across all models, the combined effect from the RW2 and Time IID terms accounts for an average of 4.3% of the variability, indicating that the overall trend in the RI-specific MCV1 coverage remains relatively stable during the time period of interest. The survey-specific effects and the spatial IID effects also account for relatively small proportions of the variability. The majority of the variation is explained by the structured spatial ICAR terms and the space-time interactions: the combined effect of the two terms accounts for an average of 88.5% of the variability, and is relatively consistent across all models. The relative proportions explained by the space-time interaction term depend on the specification of the interaction structure. This suggests that the distribution of variability between the spatial ICAR term and the space-time interaction can be sensitive to model specification. It also implies that spatial variation plays an important role in explaining the RI-specific MCV1 coverage.

We now focus on the ICAR-AR1 Model and visualize the estimated random effects in the model. The posterior medians of the ICAR spatial random effects α i and the IID spatial random effects γ i are shown in Figure 3 . The ICAR random effects show apparent spatial dependencies across states, with the northern states having lower estimated values than the southern states. In contrast, and as expected, the IID random effects show no apparent spatial pattern, and their magnitudes are much smaller compared to the ICAR random effects. This suggests that the overall spatial variation in MCV1 coverage is mostly structured, agreeing with the breakdown of variation discussed above.

The posterior medians of the RW2 temporal random effects δ b and the IID random effects τ b are shown in Figure 4 , along with the 95% point-wise CI envelopes. The RW2 temporal effects show slight downward trends at the beginning and the end of the time period, and a slight upward trend in the middle. Due to data sparsity, the uncertainties associated with the estimates at the beginning and the end of the time period are relatively large compared to those in the middle of the time period.

In Figure 5 and Figure 6 , we show the estimated space-time interactions across space and over time for all states. At each time point, the spatial dependency across states is moderate, with northern states having lower values. The temporal trends in most states were relatively flat without large fluctuations, except for a few northern states in which there were large changes over time.

The estimated survey-specific effects s are shown in Figure 7 . The estimates suggest that survey effects are not so large -a result that is also reflected in Table 3 , where the survey effects are the third largest contributor to the variance explained by the random effects (6.9%) in the ICAR-AR1 Model model.

Based on the fitted ICAR-AR1 Model, we calculated the state-level RI-specific MCV1 coverage estimates and the associated 95% CIs for the birth cohorts born between 2000 and 2018 using the procedure described in Section 3.3. Figure 8 shows the model results, along with the design-based overall MCV1 coverage estimates from each survey for reference. There is considerable spatial variation in the estimated RI-specific MCV1 coverage, with southern states generally having higher coverage than the northern states. In most states, the RI-specific coverage remained at a constant level without significant changes in trend between 2003 and 2015. The 95% CIs associated with the RI coverage estimates during this time period are relatively narrow. In contrast, the RI-specific coverage estimates in 2000-2002 and 2016-2018 show greater temporal variation, but they are associated with wider CIs. Therefore, although there seem to be a slight downward trend in many states, one should be mindful of the uncertainties associated with the trend, especially when considering prediction for the 2018 birth cohorts.

In summary, this analysis illustrates how our model is able to combine information from multiple household surveys to produce spatially-and temporally-smoothed RI-specific MCV1 coverage estimates while accounting for the impact of SIAs as well as survey-specific effects.

4.3. Sensitivity analysis: using 12-month birth cohorts. We chose to use 6-month birth cohorts in our model as they gave a good balance between two competing objectives: on one hand, we want larger birth cohorts so that there would be enough samples from each sub-national area in each survey to provide design-based coverage estimates with reasonably high precision; on the other hand, smaller birth cohorts are more desirable because they are more likely to have children with the same RI and SIA history, and hence be eligible for analysis. Here we present a sensitivity analysis using 12-month birth cohorts.

We discretize time into 12-month intervals and let the beginning of 2000 be the reference starting point. We carry out the same data processing step outlined in Section 3.1 and calculate the design-based HT estimates of the state-level overall MCV1 coverage for each eligible birth cohort in each survey. Figure 9 shows the resultant estimates with the associated confidence intervals across birth cohorts for Nigeria's 37 states. As expected, the estimates for 12-month birth cohorts generally have narrower 95% CIs compared to the estimates for 6-month birth cohorts shown in Figure 2 . However, only 61900 of the 86033 (71.9%) survey samples that were eligible for the previous analysis remain eligible for this analysis, and the proportion of the RI Only design-based coverage estimates (among the RI Only and RI + 1 SIA design-based coverage estimates) decreased from 608/1802 (33.7%) to 133/760 (17.5%). In particular, 35 of the 37 states do not have any RI Only design-based coverage estimate for any birth cohort born after 2012. As a result, we see that many southern states do not have any eligible data for birth cohorts born between 2006 and 2008, and a majority of the designbased overall MCV1 coverage estimates were influenced by a previous SIA.

The six space-time smoothing models specified in the main analysis were fitted to the design-based overall MCV1 coverage estimates for the 12-month birth cohorts. The posterior medians and 95% CIs of the model parameters as well as the model assessment metrics are reported in Table S3 of the Supplementary Material. Compared to the results in the main analysis, the estimates of the intercept β 0 are considerably higher and the estimates of the coefficient for the SIA indicator β 1 are lower. The 95% CIs for these two parameters are also wider across all models. These results are likely due to the fact that less data, especially RI Only design-based coverage data, is used in this analysis. Since less information is available for the model to draw inference on the parameters, the resultant estimates differ from those in the main analysis and have higher associated uncertainties.

Similar to the main analysis, the ICAR-AR1 Model performs the best among the six models based on all three model metrics considered. The proportions of variance explained by the various random effects also show similar patterns compare to the main analysis. The details are provided in Table S4 of the Supplementary Material.

We now focus on the ICAR-AR1 Model and compare the results from fitting this model to the 12-month birth cohort data with the ICAR-AR1 Model results in the main analysis. The estimated ICAR and IID spatial random effects and the survey specific effects show very FIG 9. Design-based HT estimates of the overall MCV1 coverage across birth cohorts from each survey for Nigeria's 37 states in the sensitivity analysis with 12-month birth cohort data. similar patterns between the two analyses. The estimated space-time interactions also show similar general trends, but some small-scale fluctuations shown in the main analysis are missing in this analysis, especially during time periods where RI Only design-based coverage data are not available. The estimated RW2 temporal random effects in this analysis are following a slightly different trend compared to the corresponding estimates in the main analysis: instead of showing a slight downward trend after 2012, the estimated RW2 random effects follow a slight increasing trend. In addition, the IID temporal random effects have considerably wider uncertainty intervals than those in the main analysis. Plots of these estimates and the associated uncertainties can be found in Section S5 of the Supplementary Material.

In Figure 10 , we compare the state-level RI-specific MCV1 coverage estimates produced by fitting the ICAR-AR1 Model to the 6-month and 12-month birth cohort data. The estimates from the 12-month birth cohort analysis are generally higher than those from the 6-month analysis. In particular, the differences between the two analyses tend to be large during time periods in which (1) there is no design-based coverage data (e.g., between 2006 and 2008), or (2) there is no RI Only coverage data (e.g., after 2012) for fitting the space-time smoothing model in the 12-month birth cohort analysis. The pattern we see in this comparison, as well as the difference in parameter and random effect estimates we showed earlier, are likely due to the fact that using 12-month birth cohorts resulted in considerably less survey data being eligible for the smoothing model, hence some temporal patterns that were captured in the 6-month birth cohort analysis were unable to be estimated in this analysis. 5. Discussion. In this paper, we developed a space-time smoothing model for estimating RI-specific MCV1 coverage at the sub-national level over time using data from household surveys. The proposed method is able to account for the impact of SIAs on the overall MCV1 coverage measured by cross-sectional surveys, and makes efficient use of data from not only birth cohorts who only had RI opportunities, but also those who had RI plus one SIA opportunity at the time of surveys. In addition, our model can be applied to data from multiple surveys with different sampling designs and construct coverage estimates with uncertainty that reflects the various data collection schemes. Any survey-specific effects such as sampling biases are also accounted for when predicting underlying RI coverage over time. We applied our method to analyze data from nine household surveys conducted in Nigeria and constructed the RI-specific MCV1 coverage in each of Nigeria's 37 states from 2003 to 2018.

Since the DHS and MICS regularly conduct household surveys in most LMICs, our method is readily applicable to many countries with high measles burden to provide reliable RI-specific MCV1 coverage estimates over time to track routine vaccine delivery system strength. It helps identify areas in which the RI program is the weakest and facilitates balanced implementation of RI and SIAs to improve overall measles vaccination coverage. In addition, our method can be easily modified to produce sub-national coverage estimates for other types of vaccines, including those not influenced by SIAs such as the diphtheria, pertussis, and tetanus (DPT) vaccines and the Bacille Calmette-Guerin (BCG) vaccine.

Under our framework, the RI-specific MCV1 coverage has a clear and sensible definition -the percentage of children in a birth cohort who receive a first MCV dose according to the RI schedule. This definition is particularly relevant in the context of building measles epidemiological models and evaluating SIA efficacy, where the RI-specific coverage is an essential input for calculating the number of children entering the susceptible pool over time (Verguet et al., 2015; Thakkar et al., 2019; . We emphasize here that the denominator of this coverage only includes children who are still alive at their scheduled RI date, i.e., at the age of 9 months. In other words, our coverage definition excludes those who died before they become eligible to receive MCV1 through RI, and hence, does not require consideration of the mortality rate among those under the age of 9 months.

Our method does make assumptions about the mortality rates among children above the age of 9 months. Specifically, we assume that there is negligible difference in mortality rate between the children who are vaccinated against measles and the children who are not. This allows us to assume the underlying MCV1 coverage within each birth cohort remains constant as the children grow older, unless they have an SIA opportunity. As such, we can treat each cross-sectional survey as a snap-shot measurement of the same underlying MCV1 coverage within a birth cohort at the time of survey interview, and use survey data collected from children in older age groups for our model. However, if the unvaccinated children suffer significantly higher mortality than vaccinated children, the proportion of the surviving children in a birth cohort who are vaccinated would increase as the children grow older. In this case, the underlying MCV1 coverage within each birth cohort no long remains unchanged over time, and the coverage estimates from surveys implemented at different times among the same birth cohort are measuring different underlying coverage values. Therefore, a more complex model as well as additional mortality data are required to properly account for the impact of discrepancies in mortality rate between vaccinated and unvaccinated children, which would be an important topic for future research.

In addition, our model also assumes that children have their RI opportunity according to the stated schedule, and SIAs were carried out promptly as recorded in the calendar. These assumptions may not hold in reality, especially in settings where vaccine stock outs or disruption of services lead to delayed RI and SIAs. Their impact on RI-specific coverage estimation would depend on when and where SIAs and surveys were carried out. Further simulation studies need to be conducted to investigate the issue under various scenarios. More recently, the COVID-19 pandemic has disrupted routine immunization services on an unprecedented scale, and many countries have temporarily suspended measles SIAs due to risk of transmission and the need to maintain physical distancing (World Health Organization, 2020c). Region-specific adjustment will be required to account for the impact of the pandemic when estimating measles vaccination coverage during the relevant time periods.

Last, but not least, we assume there is negligible systematic recall bias across age groups in surveys. However, it is possible that care givers would under-or over-report doses for older children due to memory recall. Previous studies gave mixed results regarding the presence, direction and magnitude of care giver recall bias in vaccination coverage surveys (Porth et al., 2019; Valadez and Weld, 1992; Ramakrishnan et al., 1999; Hu et al., 2019; Binyaruka and Borghi, 2018) , and the results vary greatly by geographical regions. Future research could focus on how to detect and account for recall bias, especially given the specific context of measles vaccination.

Supplement to "Space-time smoothing models for sub-national measles routine immunization coverage estimation with complex survey data" (doi: TBD; .pdf).

Bayesian image restoration, with two applications in spatial statistics

Strengthening routine immunization through measles-rubella elimination

Validity of parental recalls to estimate vaccination coverage: evidence from Tanzania. BMC health services research

Monitoring vaccination coverage: defining the role of surveys

Comparison of inflation of third dose diphtheria tetanus pertussis (DTP3) administrative coverage to other vaccine antigens

Impact of state weights on national vaccination coverage estimates from household surveys in Nigeria

Supplement to âȂIJSpace-time smoothing models for sub-national measles routine immunization coverage estimation with complex survey data

Estimating efficacy of measles supplementary immunization activities via discrete-time modeling of disease incidence time series

Estimates of income for small places: an application of James-Stein procedures to census data

Measuring Mortality, Nutritional Status, and Food Security in Crisis Situations: SMART Methodology Manual v1.0. SMART, Action Against Hunger Canada

Measuring Mortality, Nutritional Status, and Food Security in Crisis Situations: SMART Methodology Manual v2.0. SMART, Action Against Hunger Canada

Measuring coverage in MNCH: tracking progress in health for women and children using DHS and MICS household surveys

A generalization of sampling without replacement from a finite universe

Validity of maternal recall to assess vaccination coverage: evidence from six districts in Zhejiang province

The Demographic and Health Surveys. The DHS Program website

Empowering countries to enhance immunization and overall health service delivery through improved data collection, quality, and use

Bayesian modelling of inseparable space-time variation in disease risk

Global Measles and Rubella Strategic Plan

Routine Immunization. https:// measlesrubellainitiative.org/learn/the-impact/routine-immunization

Space-time smoothing of complex survey data: small area estimation for child mortality

Immunization Agenda 2030: A Global Strategy to Leave No One Behind

Immunization, vaccines and biologicals -Data, statistics and graphics

At least 80 million children under one at risk of diseases such as diphtheria, measles and polio as COVID-19 disrupts routine vaccination efforts, warn Gavi, WHO and UNICEF

Childhood immunization in Ethiopia: accuracy of maternal recall compared to vaccination cards

Magnitude of recall bias in the estimation of immunization coverage and its determinants

Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations

Penalising model component complexity: A principled, practical approach to constructing priors

R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing

Decreasing measles burden by optimizing campaign timing

The UNICEF Multiple Indicator Cluster Surveys

High resolution age-structured mapping of childhood vaccination coverage in low and middle income countries

A spatial regression model for the disaggregation of areal unit based data to high-resolution grids with application to vaccination coverage mapping

Mapping vaccination coverage to explore the effects of delivery mechanisms and inform vaccination strategies

Geospatial variation in measles vaccine coverage through routine and campaign strategies in Nigeria: Analysis of recent household surveys

Maternal recall error of child vaccination status in a developing nation

Controlling measles using supplemental immunization activities: a mathematical model to inform optimal policy