#
Estimating Under Five Mortality in Space and Time

in a Developing World Context

###### Abstract

Accurate estimates of the under-5 mortality rate (U5MR) in a developing world context are a key barometer of the health of a nation. This paper describes new models to analyze survey data on mortality in this context. We are interested in both spatial and temporal description, that is, wishing to estimate U5MR across regions and years, and to investigate the association between the U5MR and spatially-varying covariate surfaces. We illustrate the methodology by producing yearly estimates for subnational areas in Kenya over the period 1980–2014 using data from demographic health surveys (DHS). We use a binomial likelihood with fixed effects for the urban/rural stratification to account for the complex survey design. We carry out smoothing using Bayesian hierarchical models with continuous spatial and temporally discrete components. A key component of the model is an offset to adjust for bias due to the effects of HIV epidemics. Substantively, there has been a sharp decline in U5MR in the period 1980–2014, but large variability in estimated subnational rates remains. A priority for future research is understanding this variability. Temperature, precipitation and a measure of malaria infection prevalence were candidates for inclusion in the covariate model.

Keywords: Complex Surveys, Space-Time Smoothing, Stratified Cluster Sampling, Under-5 Mortality Rates

## 1 Introduction

Currently UNICEF estimates the under-five child mortality rate (U5MR) at the national level (which is known as Admin 0), using the Bayesian B-spline bias-reduction (B3) method (Alkema et al., 2014; Alkema and New, 2014). However, subnational variation is of great interest, and has been highlighted as such in the Sustainable Development Goals (SDGs). SDG 3.2 states, “By 2030, end preventable deaths of newborns and children under 5 years of age, with all countries aiming to reduce neonatal mortality to at least as low as 12 per 1,000 live births and under-5 mortality to at least as low as 25 per 1,000 live births.” From https://sustainabledevelopment.un.org/post2015/transformingourworld, with reference to review processes, paragraph 74.g states, “They will be rigorous and based on evidence, informed by country-led evaluations and data which is high-quality, accessible, timely, reliable and disaggregated by income, sex, age, race, ethnicity, migration status, disability and geographic location and other characteristics relevant in national contexts.”

In much of the developing world, there is no vital registration, and estimates of U5MR are based on survey data. In this paper, we carry out detailed analyses of such data from Kenya. Intervention and subnational policy are generally implemented at Admin 2, the second administrative level. This naming convention can be confusing with the 8 provinces of Kenya not appearing in the Admin hierarchy, and the 47 counties being labeled Admin 1. It is at this level that policies are implemented for Kenya, and hence is our spatial target of inference. We use data from Demographic Health Surveys (DHS). The DHS Program began in 1984 and has carried out more 300 surveys in over 90 countries. Typically stratified cluster sampling is carried out information is collected on population, health, HIV and nutrition.

We briefly review previous approaches to producing sub-national U5MR estimates. Adopting demographic notation, , so that we are interested in (note that strictly speaking U5MR is a probability rather than a rate). Dwyer-Lindgren et al. (2014) compare various spatial models for U5MR modeling in Zambia using DHS data. The logit of the U5MR is modeled as normally distributed, but with a single common variance across all studies, which is clearly inappropriate since it does not acknowledge the differing effective sample sizes in each area. Computation was carried out using the integrated nested Laplace approximation (INLA) of Rue et al. (2009). Mercer et al. (2015) analyzed DHS data from 22 regions in Tanzania and assumed a likelihood in which the logit of the weighted (design) estimator was assumed to be normally distributed with variance given by the design variance. A discrete space, discrete time (5-year intervals) model (Knorr-Held, 2000) was used to smooth the mean of this distribution, with implementation via INLA. Unfortunately in our Kenya study we require finer temporal and spatial scales, and at such scales the data are sparse and the weighted estimators are unstable, making the method of Mercer et al. (2015) unfeasible. Pezzulo et al. (2017) model across 27 countries in Sub-Saharan Africa, at the Admin 1 level. Estimation was based on the most recent DHS with the log weighted estimators assumed to be normally distributed with spatial smoothing being carried out via the model of Leroux et al. (1999). Extensive covariate modeling was carried out with potential variables being averaged within areas, and also allowing interactions by large regions (with three regions in total). The associations at the area-level cannot be transferred to the individual-level as this opens up the possibility of the ecological fallacy (Wakefield, 2008).

Burke et al. (2016) follow a different approach to modeling U5MR across sub-Saharan Africa. Kernel density estimation (KDE) is carried out with surfaces produced at a geographical scale of approximately 10km10km. This approach follows Larmarange and Bendaud (2014) who used the same method in the context of HIV prevalence estimation. Inference, including producing uncertainty surfaces, is difficult to obtain with KDE and the approach has been found to be inferior to Bayesian geostatistical modeling (Hallett et al., 2016).

More recently, Golding et al. (2017) carried out subnational estimation of U5MR for sub-Saharan Africa, with a continuous model for space. Four separate models were fitted to the age groups 0–1 months, 1–11 months, 12–35 months, 36–59 months, with the subsequent estimates being combined to give the U5MR. This combination is done by taking draws from the posteriors assuming they are independent, but they are not, since they are based on the same children. As well as full birth history data from DHS, summary birth history is also included. These latter typically consist of the number of children ever born, and the number who have died, along with the age of the mother. For full birth history the data are modeled as binomial with no explicit correction for the survey design. The summary birth histories are also assumed to be binomially distributed, with an artificial response and denominator created through an elaborate procedure. A space-time smoothing model is specified via the stochastic partial differential equations (SPDEs) formulation of Lindgren et al. (2011). The same space-time covariance parameters are assumed for the whole of Africa. Covariates are also modeled, and we give further details of the approach in Section 4. There is no adjustment for mothers lost to HIV, which can lead to serious underestimation in countries with HIV epidemics. Estimates in each spatial grid cell are adjusted so that the national total agrees with the Global Burden of Disease (GBD) estimates. The most recent GBD (GBD 2016 Mortality Collaborators, 2017) produced national estimates for 195 countries and territories over the period 1970–2016. Some of the constituent data in the study of Golding et al. (2017) do not contain GPS locations, but rather the administrative region within which the clusters were sampled. In this case, Golding et al. (2017, Supplementary Materials, Section 8) assign the data to a single point within the area, where this point is obtained as the weighted combination of representative points that are obtained through -means clustering. This approach is, at best, an approximation, since one needs to take a mixture over the likelihoods at each potential location, see Wilson and Wakefield (2017).

The rest of this paper is structured as follows. In Section 2 we describe the data that we use for analysis. Section 3 develops the method and gives the results for constructing the space-time child mortality surface, while Section 4 does the same for covariate modeling. Section 5 concludes the paper with a discussion of ways in which we would like to extend the model.

## 2 Data

### 2.1 Survey Data

To estimate child mortality in Kenya, we use data from three DHS conducted in 2003, 2008–2009 and 2014. Both the 2003 and 2008–2009 Kenya DHS were designed to give estimates for the 8 provinces, and for urban and rural regions separately. To this end, the sample was stratified by 8 provinces crossed with an urban/rural designation to yield 15 strata (Nairobi is solely urban). In each of these surveys the first sampling stage selected 400 enumeration areas (EAs) from a sampling frame constructed from the 1999 Census. In the second stage for both the 2003 and 2008–2009 surveys, 10,000 households were selected within the sampled EAs. The 2014 Kenya DHS was designed to make estimates of demographic indicators at the 47 county level, so it was stratified by the 47 counties crossed with urban/rural indicators. This yields 92 strata since Nairobi and Mombasa are both entirely urban. The first sampling stage of the 2014 survey selects 1,584 EAs (that produced data that could be used) from the 92 strata using a sampling frame developed from the 2009 Census. In the second stage, 40,300 households were sampled from the selected EAs. To estimate U5MR we use the portion of the survey devoted to retrospective birth histories. Women who slept in the house the night before, and are aged 15–49 are asked to enumerate all births with dates of birth, and for children who have died, dates of death. Birth histories are converted into person months for each child in the dataset. Using a discrete hazard model, each person month yields a Bernoulli (binary) random variable, survived/dead. Hence, we implement a discrete time event history analysis. It is important to note that each unique case can result in at most one death. Figure 1 shows the cluster locations for the three surveys along with the boundaries of the 47 counties. We see that the distribution of the sampling locations is far from uniform, reflecting population density. Reported response rates for households and women are high. Such data are potentially subject to various biases, e.g., recall bias, as the birth histories may go back many years if the woman surveyed is old. Though we have data from only three survey waves, the retrospective birth history gives us data on births over the period 1980–2014. According to the final reports from the Kenya DHS, nationally the U5MR per 1,000 births was 115 in 2003, 74 in 2008–2009 and 52 in 2014. While these estimates show a clear trend of decreasing child mortality, we would like to investigate the subnational variability across the 47 counties. Kenya provides a good test example due to large number of clusters (1,584) sampled in the 2014 DHS.

### 2.2 HIV Adjustment

Kenya has had a relatively high prevalence of HIV, and this can lead to serious bias in estimates of U5MR, particularly before antiretroviral therapy (ART) treatment became widely available. Pre-treatment HIV positive women had a high risk of dying, and such women who had given birth were therefore less likely to appear in surveys. The children of HIV positive women are also more likely to die before age 5 compared to those born to HIV negative women, and therefore we expect to underestimate U5MR if we do not adjust for the missing women, i.e., the missing data are non-ignorable.

Estimates of bias may be obtained using the cohort component projection model of Walker et al. (2012). Under this model, for a particular survey, year and province, the number of births is estimated, and these are attributed to HIV-negative and HIV-positive women, using estimates of the number of women in need of services to prevent mother-to-child transmission. The children born are then further subdivided into those that will and those that will not become infected with HIV, and survival probabilities of these children are then estimated, to produce a bias ratio. Let represent the true U5MR and the biased (unadjusted for HIV) U5MR in survey , year and province , . The Walker et al. (2012) method gives an estimate of,

(1) |

Figure 2 shows the ratios (reciprocal bias) plotted against year for each of the three surveys, and for the 8 regions of Kenya for which we have available data; we would prefer to have estimates at the 47 county level, but the constituent data are not available, and the 47 counties are nested within the 8 provinces, which eases the application of the adjustment. We see that the ratios of reported to true decrease as the HIV epidemic takes hold and then increase with the uptake of ART. Figure 3 shows maps of the ratios in 1995, and the large between-province differences are apparent. The ratios will clearly make a significant impact on our estimates, and are included in an offset in the model we describe in Section 3. A current weakness of our approach is that we do not account for the uncertainty in the manner by which the ratios were estimated.

## 3 Constructing a Space-Time Surface

### 3.1 The Space-Time Model

Survey data come from and describe a finite population. The DHS provides sampling weights for each individual that account for the selection probability and non-response. Skinner and Wakefield (2017) review the design and analysis of survey data. The design-based (or randomization) approach to inference is to couch inference in the context of repeated sampling from the fixed finite population. The word fixed is key here, the data are not viewed as random, rather the indices of the units (households, here) within the population that are sampled are the random variables. Weighted (often referred to as direct) estimators (Horvitz and Thompson, 1952) provide a design-consistent approach to estimation, but the sparsity of data in both time and space, are problematic since a greater proportion of cells with zero deaths in some age groups occur when we drill down to finer spatio-temporal units. Even with small numbers of deaths, variance estimates are unstable. This is a small area estimation (Rao and Molina, 2015) and at the scale for which inference is desired, smoothing in space and time is required.

In general, the design must be acknowledged when inference is performed, otherwise biased estimates with an incorrect measure of uncertainty will be produced. As an extreme example, in the DHS sampling is stratified by urban/rural and if a particular county only urban clusters were selected then ignoring this aspect will lead to bias in the estimation of the county level estimate, if U5MR is associated with urban/rural. As an alternative to design-based inference, a more traditional statistical approach may be employed in which a probability model for the observations is assumed, and the mean model contains terms that reflect the design, with a carefully chosen variance model. This approach is known as model-based inference; Wakefield et al. (2016) compare the two approaches via simulation in a spatial context.

As in Mercer et al. (2015) we assume a discrete hazard model, with six hazards for each of the (monthly) age bands: [0,1), [1,12), [12,24), [24,36), [36,48), [48,60]. Detailed argument in Allison (2014) show that the contributions for a generic child correspond to the product of up to 60 Bernoulli likelihoods with being a binary indicator of survival in month , , for a child in survey , in a household sampled at location in year , . For a month beginning at , the hazard is,

where links the month to the six age bands , i.e.,

The likelihood for survival from month to in survey and at location in year is,

The latent logit model consists of a part that is used for prediction, random effects to acknowledge the cluster sampling, survey and independent temporal effects, and an offset that adjusts for the bias due to HIV epidemics, given in (1). In summary,

(2) |

We now describe each of the components. More details on the HIV bias offset are given in the supplementary materials, but the adjustment is carried out at the province level, indexed by , with corresponding to the province in which the cluster at is located. The random cluster effects allowing for dependence amongst the individuals in the households at location ; hence, these effects will acknowledge the cluster design. The survey random effects allow for systematic biases in each of the three surveys (though of course this is relative to the average of the three surveys, and does not correct for any overall bias in the three surveys combined). The temporal terms allow for yearly perturbations that have no structure in time. Each of the six age bands, has its own intercept . The surveys are each stratified on an urban/rural indicator and on either 8 (years 2003 and 2008–2009) or 47 (year 2014) areas. The area-level stratification is strongly confounded with space and so we do not include a fixed effect for these strata, rather we assume the spatial field accounts for any such differences at a relatively large scale. The urban/rural classification changes far more quickly around urban centers, and for this reason we include a strata fixed effect . The temporal terms are random walks of order 2 (RW2), with one each for [0,1) and [1,12) months and then a third for the remaining period of [12,60] months. We decided on these splits based on initial analyses and on the demographic pattern in which the majority of U5MR deaths occur in the first year of life. In each, for reasons of parsimony, the same precisions were used (we investigated the use of different precision parameters for the three age groups, but there was little difference in the resultant inference), i.e., the distribution is for all three age bands. Sharing the precision parameter forces the same smoothness in the temporal evolution for the logit of the hazard in each age group, but the temporal trends are independent between age groups, conditional on the precision parameter. The RW2s have sum-to-zero constraints to make them identifiable when combined with the age-group-specific intercepts.The most complex term to explain is the space-time interaction , and we begin with a description of separable processes.

A separable spatio-temporal process has a covariance function that is a combination of a spatial dependence structure, , and a temporal dependence structure, , through

The multiplicative structure is beneficial because it is easy to construct valid spatio-temporal covariance functions by combining valid spatial and temporal covariance functions.

In our model we use a combination of a Matérn covariance structure for space, which is approximated via SPDE, and an AR(1) process in time. Inference is done using INLA with samples drawn from the approximate posterior for inference on functions of interest. The process is written as and is a combination of a temporal structure and a spatial structure, which translates to,

if the process is observed on (in which case is , is and is ).

The hazard for each age group is expected to vary spatially, but due to data sparsity the data will not support a separate spatial main effects for each of the six age bands. A parsimonious model would include a shared spatial main effect for the age groups, but since a spatio-temporal interaction is necessary to account for the yearly changes in the spatial pattern, we do not include the spatial main effect. It is too expensive to apply the necessary temporal sum-to-zero constraints that would be required to give identifiable spatial main effects alongside a spatio-temporal interaction.Therefore, the shared spatial main effect and the shared spatio-temporal interaction are both handled with a separable spatio-temporal model that combines an AR(1) structure with the Matérn covariance function. The resulting spatio-temporal covariance function can be explained through a constructive example which gives some intuition on the space-time interaction. A stable AR(1) process with marginal variance 1 can be generated by

where , for , and . The temporal process can be made spatio-temporal by replacing the starting condition and the innovations with spatial Matérn fields, to give

where , and , where is the stationary Matérn covariance function. Hence, a proportion of the marginal variance is explained by the previous time step and a proportion is arising from a new realization of a spatial field.

The joint identifiability of the six temporal trends and the spatio-temporal interaction can be achieved through integrate-to-zero constraints for each year. This integration is carried out with respect to the spatially varying population density :

where is the separable spatio-temporal process. These yearly integrate-to-zero constraints mean that the spatial average of the spatio-temporal effect is constantly equal to zero and that the temporal change in the spatial average of the logits of the hazards of each age group is explained by the corresponding temporal main effects. In particular, the RW2 trends are approximately interpretable as the change in the national level with time.

This spatio-temporal effect on a temporal resolution of 35 years is too expensive to include in the Bayesian model, but since we want the spatio-temporal process to change gradually in time, it is possible to use an approximation that changes piecewise linearly in time, a similar approach was taken in Blangiardo and Cameletti (2013, Chapter 8). We decrease the resolution of the spatio-temporal process to 8 time steps by defining for knot locations , corresponding to years , and defining

where gives the factor required for linear interpolation between the two knot locations. Note that if the integrate-to-zero constraint is satisfied for for , the integrate-to-zero constraint is also satisfied for linear combinations for .

Each of the precisions for the independent and identically distributed effects have Gamma(0.5, 0.0005) priors (which give 5%, 50%, 95% quantiles for the standard deviations of 0.016, 0.047, 0.52). The spatial part of the spatio-temporal interaction has a “penalized complexity” prior (Fuglstad et al., 2015; Simpson et al., 2017) with and ; all other parameters have default priors.

For predictions, the cluster, survey and temporal independent and identically distributed effects are not included so that the only contribution is . The predicted U5MR at location and at time is,

where , for and with given by (2).

The data and the fitted model are on a continuous spatial scale, but the aim is to produce values on a discrete scale using the 47 administrative regions. The information available is the posterior of the spatially varying U5MR and the population density . We obtained the latter from worldpop.org (Linard et al., 2012). We really should be using the births density, but such data are difficult to obtain; we examined a surface of estimated live births for one year that was available (WorldPop, 2017), and the surface for that year showed very little change. We assume that the infinite superpopulation has the same relative variation in population density as the real population and define the U5MR of region by

(3) |

where denotes administrative region . This averaging gives zero weight to areas with no population, even though the continuous surface is defined at such points.

### 3.2 Constructing a Space-Time Surface Results

We begin by summarizing inference on some of the key elements of the model, before reporting on substantive summaries. The left panel of Figure 4 shows the posterior medians of the RW2 median fits for each of the [0,1), [1,12), [12,60] age groups, along with 95% point-wise credible interval envelopes. We see that the temporal trend decreases for all three age groups. While the [0-1] age group shows a decreasing slope from 1995 onwards, a continuing strong decrease can be seen for the other two age groups, with the most prominent drop being for the 12–59 month age group. One of the major reasons for the drop in U5MR is increased vaccination (Haakenstad et al., 2016), but this is carried out once maternal antibodies are no longer present after 6–12 months. The right panel of Figure 4 shows the HIV adjusted version of this plot and the effect of the epidemic is clear to see in all three age groups.

Table 1 gives posterior summaries of key parameters in the space-time model. The standard deviations are not all comparable since for the RW2 it is a conditional standard deviation. The spatio-temporal standard deviation is relatively large indicating that there are strong spatial effects for the Kenya data and the median of the range parameter is 1.77, which is quite large (about a fifth the size of the study region). There is also strong year-to-year correlation in the AR(1) model.

Parameter | 2.5% | 50% | 97.5% |
---|---|---|---|

Standard deviation for RW2 time | 0.0097 | 0.018 | 0.031 |

Standard deviation for IID-time | 0.023 | 0.049 | 0.099 |

Range for spatio-temporal effect | 1.28 | 1.77 | 2.45 |

Standard deviation for spatio-temporal effect | 0.49 | 0.59 | 0.71 |

AR(1) parameter for spatio-temporal effect | 0.78 | 0.86 | 0.92 |

Standard deviation for IID-cluster | 0.32 | 0.36 | 0.39 |

Standard deviation for IID-survey | 0.017 | 0.045 | 0.13 |

Effect of rural vs urban | 0.011 | 0.080 | 0.15 |

Figure 5 shows a comparison between the modeled and weighted estimates at the 47 county level, and aggregated over 5 years (aggregation over years is required, otherwise the direct estimates are unstable). We see some attenuation due to shrinkage, as expected. In the Supplementary Materials we include more detailed plots and show the uncertainty in the modeled and weighted estimators. Again, as expected, the modeled estimates have much greater precision.

As mentioned in Section 1, we wish to make inference at the spatial level at which policy interventions occur. For Kenya, this is at the 47 county level, and Figure 6 shows a sequence of 9 maps of for the years (we have 35 yearly estimates, but for space reasons we take a 5-year spacing). The last two of these years are obtained by forecasting from the model. The density of hatching reflects the uncertainty. The dramatic decrease over time in is apparent, though strong subnational variation persists. The Supplementary Materials contain maps of the uncertainty.

Figure 7 shows the posterior medians of the spatio-temporal terms for the years 1980, 1985, 1990, …, 2015, 2020. The last two of these years are obtained by predicting forward the space-time field. From 1980 onwards strong spatial effects can be seen in the counties Turkana and West Pokot in the north west part, the province Nyanza in the middle west part and the counties Kilifi, Tana River and Garissa in the south east part of Kenya. While the highs in the north west and south east have almost disappeared by 2004, a higher effect in the counties Minori and Homa Bay of the province Nyanza persist and, without interventions, one would expect these trends to continue until 2020. Around 1990 to 1995 higher effects can also seen in the north east.

While it seems that the spatio-temporal trend decreases over time it should be emphasized that there is still a strong effect present in recent periods and also in the future. To illustrate this we computed the 95% and 5% points of the pixel values for each of the nine maps. In 1980 the quantile was and the quantile was leading to a ratio of . While the quantile decreases until 2005 and then increases again, the quantile decreases almost constantly. The ratio of 95% to 5% points increases until 1995 with a value of and then decreases. In 2010 the ratio is still and in the predicted years and it is and , respectively. Thus, there remains strong subnational differences in U5MR.

The Millenium Development Goals (MDG) aimed for a drop of 67% in U5MR between 1990 and 2015. In the left hand panel of Figure 8 we map the posterior median of the percentage drop at the county level, with counties in the central part of Kenya experiencing very small decreases only. In the right hand panel we plot the posterior probability that each county achieved this aim and we see that very few attained a 67% drop. Over the country the posterior median drop was 55% with 95% credible interval of (49.9%, 60.0%), and a 0% probability that Kenya achieved the goal.

To examine the performance of the space-time smoothing model, we held out some of the data and then predicted the U5MR at these points using the weighted estimate and the smoothed estimate. Specifically, we get estimates of the U5MR for all counties and periods from the model using all the 2003 and 2008–2009 DHS, along with 397 clusters from the 2014 DHS. We then calculate weighted estimates of U5MR using the remaining 1,187 clusters, and these are treated as the closest to the truth, since they are based on a large sample. Due to stability of the weighted estimates we look only at the periods 1990–1994,1995–1999, 2000–2004, 2005–2009 and 2010–2014, and form estimates for each of the 47 counties. Let denote the weighted estimator and the smoothed estimator (from our model) in county and period . We compare these estimates with the weighted estimates from the 1,187 clusters, . In particular, we calculate,

(4) |

for 1990–1994,1995–1999, 2000–2004, 2005–2009 and . Table 2 presents the MSEs and we see that in all cases the smoothing model has far superior performance.

Period | Weighted | Smoothed |
---|---|---|

1990–1994 | 55 | 34 |

1995–1999 | 36 | 13 |

2000–2004 | 22 | 6.8 |

2005–2009 | 9.1 | 3.9 |

2009–2014 | 7.4 | 2.8 |

## 4 Exploratory Covariate Modeling

### 4.1 The Covariate Model

In this section we carry out an exploratory investigation into whether any of the spatial variability we see in Kenya can be attributed to a variety of covariates. Before outlining our approach, we provide a brief literature review of suggestions for building covariate models in the setting considered here.

Gething et al. (2015) describe the use of DHS data to construct surfaces of: access to HIV testing in women, stunting in children, anemia prevalence in children and access to improved sanitation. For each outcome and each country the following procedure was carried out. A collection of 17 covariates were examined. Initially, simple linear regression was used with three versions (the original, the square and the square root) of each of the 17 variables taken. Cross-validation was then used to reduce these to a subset of 17 terms. Two-way interactions for these 17 were added to the collection to give additional terms. This complete set was reduced to 20, again via cross-validation. Then the resultant potential models, that were combinations of these 20 terms, were compared.

Bhatt et al. (2017) use an approach known as stacked generalization (Wolpert, 1992) in which multiple predicting algorithms are weighted to produce a final prediction. This approach is closely related to the more general super-learner approach (Van der Laan et al., 2007). This approach is interesting and has optimality properties for prediction but has a lack of interpretability, the model is not suitable for predictions into the future, and there are questions over whether uncertainty in the procedure can be incorporated into interval estimates for the surface. A similar approach was used by Golding et al. (2017).

The covariates we choose to examine are access (estimated travel time to cities with at least 50,000 people; Nelson, 2008), aridity (Zomer et al., 2007; Zomer et al., 2008), precipitation (Fick and Hijmans, 2017), temperature (Fick and Hijmans, 2017), enhanced vegetation index (EVI; Didan, 2015), Plasmodium falciparum parasite rate (PfPR) in children (Bhatt et al., 2015), population (Lloyd et al., 2017) and a wealth index (calculated from DHS household data for each and then spatially modeled). Further details on these covariates can be found in the Supplementary Materials. For the purposes of exploration, we model access, aridity, temperature, and precipitation as time-invariant; plots of these variables can be found in the top row of Figure 9. Data on PfPR, population, and vegetation were obtained for the years 2000–2014 and subsequently averaged within each of the three 5-year periods (2000–2004, 2005–2009 and 2010–2014) to obtain values for each period; these data are also displayed Figure 9.

In order to determine which covariates are predictive of U5MR, we will use a simplified version of the model described in Section 3.1, in which we replace the yearly model with a model over 5-year periods 2000–2004, 2005–2009, 2010–2014. The model is,

(5) |

where are the age-specific intercepts, are stratum (fixed) effects, is a temporal random effect (assumed common to all age groups) and is modeled using a RW1 (since we have three periods only), are cluster random effects and are survey random effects. In comparisons to be presented in Section 4.2 we compare four different approaches/models with referring to the direct estimates and , corresponding to choosing the “Other Variables” in (5) to be space only, covariates only and space and covariates, i.e.,

where are the spatial covariates at location and in period , and is a spatial random effect at cluster with location . The spatial model is as before, a Gaussian Markov random field with Mátern covariance function (fitted using the SPDE approach) and, for simplicity, we assume it has the same structure for every age group and period. We divide the data into training and test sets. In the training set we build the models and in the test set we compare their performance. We split the 2014 DHS into two, roughly equal-sized, groups. We use 799 clusters from the 2014 DHS as our test set (for comparison purposes). The other clusters in the 2014 DHS along with data from the 2003 and 2008/2009 DHS will be used for training the model, resulting in 1,581 clusters being used. To emphasize, the spatial model, , is fit just once, while and are fit multiple times, for each combination of covariates. For these models, we assess their performance by using the DIC (Spiegelhalter et al., 1994), CPO (Held et al., 2010), and WAIC (Watanabe, 2013) criteria. As a result exercise, we determine the best models in each of the and collections to be used to compare with the direct estimates and spatial model obtained from the training clusters. We will have a total of four final comparisons (with all estimates bases on the training data): direct estimates, a model with a spatial random effect, a model with the “best” collection of covariates, and a model with alternative “best” collection of covariates that are chosen when spatial effects are included in the model.

Under we have an estimator of the U5MR for each area and period , . Under model , the direct estimator has normal distribution , and under , we have posterior distributions with posterior means and posterior variances , . Then, with the “truth” (direct estimate from test data) ,

The best approach is that which minimizes the MSE.

### 4.2 Exploratory Covariate Modeling Results

The DIC, CPO and WAIC scores for all possible covariate combinations for models and are reproduced in the Supplementary Materials. There is good agreement between the three different assessments of model fit. For (no spatial effects and covariates), the best model was that which included temperature and PfPR. For (spatial effects and covariates), the model that included only precipitation and PfPR performed best.

The MSE and constituent squared bias and variance are shown in Figure 10. We see that , and perform much better than , with and being better than . We see in Figure 11 that the predicted surfaces are almost identical under models and . Somewhat surprisingly, the spatial standard deviation and range parameters did not change with the addition of covariates, and the strength of the associations changed little also (see Supplementary Materials). As expected, there is a strong positive association between the logit of U5MR and log PfPR. In model , the logit of U5MR also showed a positive association with temperature while in model precipitation was negatively associated.

## 5 Discussion

In this paper we have developed a continuous space/discrete time model for investigating the dynamics of U5MR in a developing world setting. We have illustrated that the model improves on the use of weighted estimates, and can provide reliable inference at the required geographical scale. However, there are a number of aspects that we aim to improve upon in future work. An adjustment for HIV epidemics is crucial, given the extent of the epidemic in Kenya (and in many other countries), and we would like to acknowledge the uncertainty in the bias correction.

The age pattern of human mortality between ages 0 and 5 years follows a regular, decreasing pattern across a wide range of overall levels. Net of level, this age pattern can be characterized by the ratio of mortality at each age compared to a reference age. Our models estimate mortality in six independent age groups, and it is possible that the age pattern that results from combining the estimates from the six models does not follow any of the regularly observed age patterns of human child mortality. In our analysis, this was not a problem (see Supplementary Materials), but we are currently working on a flexible model of the age pattern of mortality that can enforce this constraint.

We would like to include other data sources, for both Kenya and other countries. Early DHS do not contain the GPS coordinates of the sampled clusters, but rather the administrative areas within which sampling took place. We plan to extend methods presented in Wilson and Wakefield (2017) to model the location of the unknown sampling point. WAs described in Section 1, we have utilized so-called full birth history data in which the births and deaths of each child are available. Summary birth history consist of only the number of children ever born and the number who died, by age of mother. These data are easier to collect and are available in a large number of surveys and censuses. The incorporation of such data into a model-based framework is a priority for future work.

In this work we have used a continuous spatial model, whereas our major interest was to inspect results on the discrete scale for the 47 administrative regions. For this purpose we integrated over the spatial field and included the population density to produce the results at the county level. An obvious question that arises is: what advantages are there with this approach as compared to using a discrete spatial model, such as the ICAR model (Besag et al., 1991), directly? One advantage of the continuous model is that we get a smoothed estimated field giving an indication of the U5MR at a finer resolution. Furthermore, the adjustment for the survey design is implicitly integrated into the model when taking the population density into account. It is important to note that this would not be possible using a discrete spatial model. Another advantage is that when using a continuous random field we do not need to specify a neighborhood structure. The 47 administrative regions of Kenya vary widely in shape and size, and therefore in the number of neighbors, so that it is not clear how to define a sensible neighborhood dependence structure. Part of our future research will be to investigate how a discrete spatial model would perform in this setting. Here, we are particularly interested in the performance of the recently proposed model by Riebler et al. (2016) and in the comparison of the results to the continuous model presented here.

There are several limitations to the covariate modeling carried out in Section 4. For one, we use geographically-referenced covariates rather than household or individual level variables since we were interested in understanding large-scale patterns in U5MR. Therefore, we do not directly model several variables that are known to have an impact on childhood mortality such as biological factors (e.g., vaccination rates, disease prevalence), maternal demographics (e.g., age, education), and household characteristics (e.g., toilet facilities, access to water). Though spatial surfaces do exist for some of these variables (e.g., measles vaccination coverage: Takahashi et al., 2017) or surfaces could be developed based on DHS data (Gething et al., 2015), there is greater uncertainty associated with these variables, which can lead to misleading inference (Foster et al., 2012). We therefore limited the number of heavily-modeled covariates in our model. Additionally, many of these factors are associated with variables already included in our model.

The computations were run on a computing server with 32 Intel Xeon 2.7 GHz CPUs available. The full Bayesian model required around 14 hours for estimation and 19.5 hours for predictions. An empirical Bayes version of the model required around 2.5 hours for estimation and 10 hours for predictions. Supplementary materials and code to run the models described here can be found at http://faculty.washington.edu/jonno/software.html.

## Acknowledgments

Wakefield and Wilson were supported by grant R01CA095994 from the National Institutes of Health, Fuglstad and Riebler by project number 240873/F20 from the Research Council of Norway, Godwin by R01AI029168 from the National Institutes of Health and Clark by R01HD086227 from the National Institutes of Health. We would like to thank Zehang Li, Yuan Hsiao, Bryan Martin, Danzhen Yu, Lucia Hug, Leontine Alkema, Jon Pedersen, Patrick Gerland, Trevor Croft, Bruno Masquelier, Kenneth Hill, David Sharrow, Roy Burstein, Simon Hay and Jonathan Muir for providing data and helpful comments.

## References

- Alkema and New (2014) Alkema, L. and J. New (2014). Global estimation of child mortality using a Bayesian B-spline bias-reduction model. The Annals of Applied Statistics 8, 2122–2149.
- Alkema et al. (2014) Alkema, L., J. R. New, J. Pedersen, D. You, et al. (2014). Child mortality estimation 2013: an overview of updates in estimation methods by the United Nations Inter-Agency Group for Child Mortality Estimation. PloS ONE 9, e101112.
- Allison (2014) Allison, P. (2014). Event History and Survival Analysis, Second Edition, Volume 46. SAGE publications.
- Besag et al. (1991) Besag, J., J. York, and A. Mollié (1991). Bayesian image restoration with two applications in spatial statistics. Annals of the Institute of Statistics and Mathematics 43, 1–59.
- Bhatt et al. (2017) Bhatt, S., E. Cameron, S. Flaxman, D. Weiss, D. Smith, and P. Gething (2017). Improved prediction accuracy for disease risk mapping using Gaussian process stacked generalization. Journal of The Royal Society Interface 14, 20170520.
- Bhatt et al. (2015) Bhatt, S., D. Weiss, E. Cameron, D. Bisanzio, B. Mappin, U. Dalrymple, K. Battle, C. Moyes, A. Henry, P. Eckhoff, et al. (2015). The effect of malaria control on plasmodium falciparum in Africa between 2000 and 2015. Nature 526, 207–211.
- Blangiardo and Cameletti (2015) Blangiardo, M. and M. Cameletti (2015). Spatial and spatio-temporal Bayesian models with R-INLA. John Wiley and Sons.
- Burke et al. (2016) Burke, M., S. Heft-Neal, and E. Bendavid (2016). Sources of variation in under-5 mortality across sub-Saharan Africa: a spatial analysis. The Lancet Global Health 4, e936–e945.
- Didan (2015) Didan, K. (2015). MOD13A3 MODIS/Terra vegetation Indices Monthly L3 Global 1km SIN Grid V006. NASA EOSDIS Land Processes DAAC. https://doi.org/10.5067/modis/mod13a3.006.
- Dwyer-Lindgren et al. (2014) Dwyer-Lindgren, L., F. Kakungu, P. Hangoma, M. Ng, H. Wang, A. D. Flaxman, F. Masiye, and E. Gakidou (2014). Estimation of district-level under-5 mortality in Zambia using birth history data, 1980–2010. Spatial and Spatio-Temporal Epidemiology 11, 89–107.
- Fick and Hijmans (2017) Fick, S. E. and R. J. Hijmans (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37, 4302–4315.
- Foster et al. (2012) Foster, S., H. Shimadzu, and R. Darnell (2012). Uncertainty in spatially predicted covariates: is it ignorable? Journal of the Royal Statistical Society: Series C 61, 637–652.
- Fuglstad et al. (2015) Fuglstad, G.-A., D. Simpson, F. Lindgren, and H. Rue (2015). Constructing priors that penalize the complexity of Gaussian random fields. arXiv preprint arXiv:1503.00256.
- GBD 2016 Mortality Collaborators (2017) GBD 2016 Mortality Collaborators (2017). Global, regional, and national under-5 mortality, adult mortality, age-specific mortality, and life expectancy, 1970–2016: a systematic analysis for the Global Burden of Disease Study 2016. The Lancet 390, 1084–1150.
- Gething et al. (2015) Gething, P., A. Tatem, T. Bird, and C. Burgert-Brucker (2015). Creating spatial interpolation surfaces with DHS data. Technical report, ICF International. DHS Spatial Analysis Reports No. 11.
- Golding et al. (2017) Golding, N., R. Burstein, J. Longbottom, A. Browne, N. Fullman, A. Osgood-Zimmerman, L. Earl, S. Bhatt, E. Cameron, D. Casey, L. Dwyer-Lindgren, T. Farag, A. Flaxman, M. Fraser, P. Gething, H. Gibson, N. Graetz, L. Krause, X. Kulikoff, S. Lim, B. Mappin, C. Morozoff, R. Reiner, A. Sligar, D. Smith, H. Wang, D. Weiss, C. Murray, C. Moyes, and S. Hay (2017). Mapping under-5 and neonatal mortality in Africa, 2000–15: a baseline analysis for the Sustainable Development Goals. The Lancet. Available online, September 25th, 2017.
- Haakenstad et al. (2016) Haakenstad, A., M. Birger, L. Singh, P. Liu, S. Lim, M. Ng, and J. Dieleman (2016). Vaccine assistance to low-and middle-income countries increased to $3.6 billion in 2014. Health Affairs 35, 242–249.
- Hallett et al. (2016) Hallett, T., S.-J. Anderson, C. A. Asante, N. Bartlett, V. Bendaud, S. Bhatt, C. Burgert, D. F. Cuadros, J. Dzangare, D. Fecht, et al. (2016). Evaluation of geospatial methods to generate subnational HIV prevalence estimates for local level planning. AIDS 30, 1467–1474.
- Held et al. (2010) Held, L., B. Schrödle, and H. Rue (2010). Posterior and cross-validatory predictive checks: A comparison of MCMC and INLA. In T. Kneib and G. Tutz (Eds.), Statistical Modeling and Regression Structures – Festschrift in Honour of Ludwig Fahrmeir, pp. 91–110. Physica-Verlag.
- Horvitz and Thompson (1952) Horvitz, D. and D. Thompson (1952). A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association 47, 663–685.
- Knorr-Held (2000) Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine 19, 2555–2567.
- Larmarange and Bendaud (2014) Larmarange, J. and V. Bendaud (2014). HIV estimates at second subnational level from national population-based surveys. AIDS 28, S469–S476.
- Leroux et al. (1999) Leroux, B., X. Lei, and N. Breslow (1999). Estimation of disease rates in small areas: A new mixed model for spatial dependence. In M. Halloran and D. Berry (Eds.), Statistical Models in Epidemiology, the Environment and Clinical Trials, pp. 179–192. New York: Springer.
- Linard et al. (2012) Linard, C., M. Gilbert, R. W. Snow, A. M. Noor, and A. J. Tatem (2012). Population distribution, settlement patterns and accessibility across africa in 2010. PloS One 7, e31743.
- Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Linström (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic differential equation approach (with discussion). Journal of the Royal Statistical Society, Series B 73, 423–498.
- Lloyd et al. (2017) Lloyd, C., A. Sorichetta, and A. Tatem (2017). High resolution global gridded data for use in population studies. Scientific Data 4, 170001.
- Mercer et al. (2015) Mercer, L., J. Wakefield, A. Pantazis, A. Lutambi, H. Mosanja, and S. Clark (2015). Small area estimation of childhood of childhood mortality in the absence of vital registration. Annals of Applied Statistics 9, 1889–1905.
- Nelson (2008) Nelson, A. (2008). Estimated travel time to the nearest city of 50,000 or more people in year 2000. Technical report, Global Environment Monitoring Unit - Joint Research Centre of the European Commission, Ispra, Italy.
- Pezzulo et al. (2017) Pezzulo, C., E. Utazi, T. B. A. Sorichetta, A. Tatem, J. Yourkavitch, T. Pullum, and C. Burgert-Brucker (2017). Subnational modelling of child mortality and its drivers across 27 countries in Sub-Saharan Africa. Technical report, Paper presented at PAA Meeting.
- Rao and Molina (2015) Rao, J. and I. Molina (2015). Small Area Estimation, Second Edition. New York: John Wiley.
- Riebler et al. (2016) Riebler, A., S. H. Sørbye, D. Simpson, and H. Rue (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research 25(4), 1145–1165.
- Rue et al. (2009) Rue, H., S. Martino, and N. Chopin (2009). Approximate Bayesian inference for latent Gaussian models using integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B 71, 319–392.
- Simpson et al. (2017) Simpson, D., H. Rue, A. Riebler, T. Martins, and S. Sørbye (2017). Penalising model component complexity: A principled, practical approach to constructing priors (with discussion). Statistical Science 32, 1–28.
- Skinner and Wakefield (2017) Skinner, C. and J. Wakefield (2017). Introduction to the design and analysis of complex survey data. Statistical Science 32, 165–175.
- Spiegelhalter et al. (1994) Spiegelhalter, D., L. Freedman, and M. Parmar (1994). Bayesian approaches to randomized trials (with discussion). Journal of the Royal Statistical Society, Series A 157, 357–416.
- Takahashi et al. (2017) Takahashi, S., C. J. E. Metcalf, M. J. Ferrari, A. Tatem, and J. Lessler (2017). The geography of measles vaccination in the african great lakes region. Nature Communications 8.
- Van der Laan et al. (2007) Van der Laan, M. J., E. C. Polley, and A. E. Hubbard (2007). Super learner. Statistical Applications in Genetics and Molecular Biology 6.
- Wakefield (2008) Wakefield, J. (2008). Ecologic studies revisited. Annual Review of Public Health 29, 75–90.
- Wakefield et al. (2016) Wakefield, J., D. Simpson, and J. Godwin (2016). Comment: Getting into space with a weight problem. Discussion of, “Model-based geostatistics for prevalence mapping in low-resource settings”, by P.J. Diggle and E. Giorgi. Journal of the American Statistical Association 111, 1111–1119.
- Walker et al. (2012) Walker, N., K. Hill, and F. Zhao (2012). Child mortality estimation: methods used to adjust for bias due to AIDS in estimating trends in under-five mortality. PLoS Med 9, e1001298.
- Watanabe (2013) Watanabe, S. (2013). A widely applicable Bayesian information criterion. Journal of Machine Learning Research 14, 867–897.
- Wilson and Wakefield (2017) Wilson, K. and J. Wakefield (2017). Pointless continuous spatial surface reconstruction. arXiv:1709.09659.
- Wolpert (1992) Wolpert, D. (1992). Stacked generalization. Neural Networks 5, 241–259.
- WorldPop (2017) WorldPop (2017). Kenya 1km births, version 2. Technical report, University of Southampton, DOI: 10.5258/SOTON/WP00349.
- Zomer et al. (2007) Zomer, R. J., D. A. Bossio, A. Trabucco, L. Yuanjie, D. C. Gupta, and V. P. Singh (2007). Trees and water: smallholder agroforestry on irrigated lands in Northern India. Technical report, International Water Management Institute.
- Zomer et al. (2008) Zomer, R. J., A. Trabucco, D. A. Bossio, and L. V. Verchot (2008). Climate change mitigation: A spatial analysis of global land suitability for clean development mechanism afforestation and reforestation. Agriculture, Ecosystems and Environment 126, 67–80.