Space-Time Smoothing of Survey Outcomes using the R Package SUMMER

Consider a study region that is partitioned into $n$ areas, with interest focusing on estimating the prevalence of a binary indicator in each area, possibly over time. The data are collected via some complex survey design. For each individual $j$, let $y_{j}$ denote the individual’s outcome, and $w_{j}$ denote the design weight associated with this individual. Further, let $s_{it}$ represent the indexes of individuals sampled in area $i$ and in time period $t$. The design-based estimator (Horvitz and Thompson 1952; Hájek 1971) is

$$\hat{p}_{it}^{\texttt{HT}}=\frac{\sum_{j\in s_{it}}w_{j}y_{j}}{\sum_{j\in s_{it}}w_{j}}.$$

(1)

This is an example of a direct estimate. The variance of $\hat{p}_{it}^{\texttt{HT}}$ can be calculated using standard methods (Wolter 2007) and can be easily computed using the survey package. Let $V^{\texttt{HT}}_{it}$ denote the design-based variance of $\mbox{logit}(p_{it}^{\texttt{HT}})$, obtained from the design-based variance of $p_{it}^{\texttt{HT}}$ via linearization (the delta method). We take the logit transformed direct estimates as input data and estimate the true prevalence with the random effects model,

	$\displaystyle\mbox{logit}(\hat{p}_{it}^{\texttt{HT}})|\lambda_{it}$	$\displaystyle\sim\mbox{Normal}(\lambda_{it},V^{\texttt{HT}}_{it}),$		(2)
	$\displaystyle\lambda_{it}$	$\displaystyle=x_{it}^{\intercal}\beta+\alpha_{t}+\epsilon_{t}+S_{i}+e_{i}+\delta_{it}.$		(3)

In this model, which is a space-time smoothing extension of the Fay and Herriot (1979) model, $\mbox{expit}(\lambda_{it})$ is the true prevalence we aim to estimate, and $x_{it}$ are area-level covariates. The rest of the terms are normally distributed random effects including structured time trends $\alpha_{t}$, unstructured, independent and identically distributed (iid), temporal terms $\epsilon_{t}$, structured spatial trends $S_{i}$, unstructured spatial terms $e_{i}$, and space-time interaction terms $\delta_{it}$. The terms $e_{i}+S_{i}$ are implemented via the BYM2 parameterization (Riebler et al. 2016), a reparameterization of the classical BYM model (Besag, York, and Mollié 1991) that combines iid error terms with intrinsic conditional autoregressive (ICAR) random effects. Several different temporal models are implemented in SUMMER for the structured temporal trends and space-time interaction effects, including random walks of order 1 and 2, and autoregressive models (Håvard Rue and Held 2005) with additional linear trends. The interaction term $\delta_{it}$ can be one of the type I to IV interactions of the chosen temporal model and the ICAR model in space, as described in Knorr-Held (2000). In order for the model to be identifiable, we impose sum-to-zero constraints on each group of random effects. More details on the prior choices are provided in the supplementary materials.

For composite indicators such as mortality rates, the direct estimates require additional modeling. Here we focus on the estimation of the U5MR, one of the most critical and widely available population health indicator. The methodology and the functions in SUMMER are readily applicable to mortality rates of other age groups as well, but we note that modeling mortality beyond age 5 is usually more challenging in practice because death becomes rarer, and survey data alone are not sufficient for reliable inference.

The SUMMER package implements the discrete hazards model described in Mercer et al. (2015). We use discrete time survival analysis to estimate age-specific monthly probabilities of dying in user-defined age groups. We assume constant hazards within the age bands. The default choice uses the monthly age bands

$$[0,1),[1,12),[12,24),[24,36),[36,48),[48,60)$$

for U5MR and they can be easily specified by the user. The U5MR for area $i$ and time $t$ can be calculated as,

$$\hat{p}_{it}^{\texttt{HT}}={}_{60}{\hat{q}}^{it}_{0}=1-\prod_{a=1}^{6}\left(1-{}_{n_{a}}{\hat{q}}^{it}_{x_{a}}\right),$$

(4)

where $x_{a}$ and $n_{a}$ are the start and end of the $a$-th age group, and ${}_{n_{a}}{q}^{it}_{x_{a}}$ is the probability of death in age group $[x_{a},x_{a}+n_{a})$ in area $i$ and time $t$, with ${}_{n_{a}}{\hat{q}}^{it}_{x_{a}}$ the estimate of this quantity. This calculation follows the synthetic cohort life table approach in which mortality probabilities for each age segments based on real cohort mortality experience are combined. This allows the full use of the most recent data, which is especially useful when survey data are sparse and is the default approach that The DHS Program (2020) uses.

The constant one-month hazards in each age band can be estimated by a weighted logistic regression model (Binder 1983):

$$\mbox{logit}\left({}_{1}{q}^{it}_{m}\right)=\beta_{a[m]}^{it},$$

(5)

where $a[m]$ is the age band indicator for the $m$-th month, i.e.

$$a[m]=\left\{\begin{array}[]{ll}1&\mbox{if }m=0,\\ 2&\mbox{if }m=1,\dots,11,\\ 3&\mbox{if }m=12,\dots,23,\\ 4&\mbox{if }m=24,\dots,35,\\ 5&\mbox{if }m=36,\dots,47,\\ 6&\mbox{if }m=48,\dots,59.\end{array}\right.$$

(6)

The design-based variance of $\mbox{logit}(\hat{p}_{it}^{\texttt{HT}})$ may then be estimated using the delta method or resampling methods such as the jackknife (J. Pedersen and Liu 2012). The smoothing of the direct estimates can then proceed using the model described in equations (2) – (3). When multiple surveys exist, one may choose to either model the survey-specific effects as fixed or random (Mercer et al. 2015) or first aggregate the direct estimates from multiple surveys to obtain a “meta-analysis” estimate in each area and time period (Li et al. 2019), i.e., at each time $t$, we combine the $K_{t}$ available direct estimates from multiple surveys to form the estimate

$$\hat{p}_{it}^{\texttt{meta}}=\mbox{expit}\Big{(}\sum_{k=1}^{K_{t}}\frac{(\hat{V}_{it,k}^{\texttt{HT}})^{-1}}{\sum_{k^{\prime}=1}^{K_{t}}(\hat{V}_{it,k^{\prime}}^{\texttt{HT}})^{-1}}\mbox{logit}(\hat{p}_{it}^{\texttt{HT}})\Big{)},$$

and the associated design-based variance on the logit scale is $\Big{(}\sum_{k^{\prime}=1}^{K_{t}}(\hat{V}_{it,k^{\prime}}^{\texttt{HT}})^{-1}\Big{)}^{-1}$. To mitigate the sparsity of available data in each year, Li et al. (2019) also considers a temporal model defined at the yearly level while the direct estimates are calculated at multi-year periods. All these variations can be fit using the SUMMER package.

The space-time Fay-Herriot estimates are useful when there are enough observations at the spatial and temporal unit of the analysis. When the target of inference is at finer resolution, e.g., on a yearly time scale with admin-2 areas and surveys stratified at admin-1 levels, the direct estimates may contain many $0$s or $1$s and the design-based variance cannot be calculated reliably. In this case, we can consider unit-level models where the individual survey responses are modeled. In the rest of this section, we describe a model for the cluster-level risk, where we account for the additional within-cluster variation by allowing overdispersion in the likelihood. More detailed comparisons of this modeling choice was examined in (T. Q. Dong and Wakefield 2021). In a two-stage cluster design, the clusters are referred to as primary sampling units (PSUs) and the households are referred to as secondary sampling units (SSUs). Thus we refer to such models as cluster-level model to avoid confusion. We describe the model for the mortality estimation problem below, while the same formulation applies to the case of any generic binary indicators as well.

In the most general setting, we consider multiple surveys over time, indexed by $k$. The sampling frame that was used for survey $k$, will be denoted by $r[k]$. We assume a discrete hazards model as before. We consider a beta-binomial model for the probability (hazard) of death from month $m$ to $m+1$ in survey $k$ and at cluster $c$ in year $t$. This model allows for overdispersion relative to the binomial model. Assuming constant hazards within age bands, we assume the number of deaths occurring within age band $a[m]$, in cluster $c$, time $t$, and survey $k$ follow the beta-binomial distribution,

$$Y_{a[m],k,c,t}~{}|~{}p_{a[m],k,c,t}\sim\mbox{BetaBinomial}\left(~{}n_{a[m],k,c,t}~{},~{}p_{m,k,c,t}~{},d~{}\right),$$

(7)

where $p_{m,k,c,t}$ is the monthly hazard at $m$-th month of age, in cluster $c$, time $t$, and survey $k$ and $d$ is the overdispersion parameter. The latent logistic model we use is,

	$\displaystyle p_{m,k,c,t}=$	$\displaystyle\mbox{expit}(\alpha_{m,c,k,t}+\epsilon_{t}+b_{k}),$		(8)
	$\displaystyle\alpha_{m,k,c,t}=$	$\displaystyle\beta_{a[m],r[k],t}I(s_{c}\in\mbox{ rural })+\gamma_{a[m],r[k],t}I(s_{c}\in\mbox{ urban })$
		$\displaystyle\;+S_{i[s_{c}]}+e_{i[s_{c}]}+\delta_{i[s_{c}],t}+\mbox{BIAS}_{k,t}.$		(9)

This form consists of a collection of terms that are used for prediction and a number that are not, as we now describe. We include a survey fixed effect $b_{k}$ with the constraint $\sum_{k}b_{k}\mathbf{1}_{r[k]=r}=0$ for each sampling frame $r$, so that the main temporal trends are identifiable for each sampling frame. The $b_{k}$ terms are not included in the prediction, i.e., they are set to zero. The $\epsilon_{t}$ are unstructured temporal effects that allow for perturbations over time. It is a contextual choice whether they are used in predictions. We include terms in (9) that are analogous to those in equations (2)–(3), in particular the spatial main effects $S_{i}$ and $e_{i}$ and the space-time interactions $\delta_{it}$.

For the temporal main effects $\beta_{a[m],r[k],t}$ and $\gamma_{a[m],r[k],t}$, we have stratum-specific distinct trends for each age group $a[m]$ in surveys from each sampling frame. We include separate urban and rural temporal terms to acknowledge the sampling design, often urban clusters are oversampled and have different risk to rural clusters, and so it is important to acknowledge this aspect in the model (Paige et al. 2020). The urban-rural stratification effects may also be parameterized as time-invariant fixed effects, i.e., restricting $\beta_{a[m],r[k],t}=\gamma_{a[m],r[k],t}+\Delta_{a[m],r[k]}$. For a detailed discussion of the parameterization of stratification effects, we refer readers to Wu et al. (2021). In addition, it is usually reasonable to assume shared temporal trends up to a constant shift across some age groups. For example, we may let

$$\beta_{a[m],r[k],t}=\beta_{a[m],r[k],0}+\beta^{\star}_{a^{\star}[m],r[k],t}$$

where $\beta^{\star}_{a^{\star}[m],r[k],t}$ is a collection of temporal random effects with sum-to-zero constraint $\sum_{t}\beta^{\star}_{a^{\star}[m],r[k],t}=0$, and $a^{\star}[m]$ is a reduced set of age bands. The default choice for U5MR in the package is

$$a^{\star}[m]=\left\{\begin{array}[]{ll}1&\mbox{if }m=0,\\ 2&\mbox{if }m=1,\dots,11,\\ 3&\mbox{if }m=12,\dots,59.\end{array}\right.$$

(10)

That is, we assume the temporal trends for logit hazards in the last four age groups are parallel and only differ by the intercept term $\beta_{a[m],r[k],0}$.

In situations where biases are known for particular surveys and/or years, we can adjust for bias following Wakefield et al. (2019) by including the bias ratio term, $\mbox{BIAS}_{k,t}=\mbox{U5MR}^{\star}_{t}/\widehat{\mbox{U5MR}}_{k,t}$, where $\mbox{U5MR}^{\star}_{t}$ is the expected U5MR in year $t$ and $\widehat{\mbox{U5MR}}_{k,t}$ is the biased version. This approach has been used to adjust for mothers who have died from AIDS (Walker, Hill, and Zhao 2012); such mothers cannot be surveyed, and their children are more likely to have died, so the missingness is informative.

The predicted U5MRs in urban and rural regions of area $i$ and at time $t$ according to sampling frame $r$ are,

	$\displaystyle\mbox{U5MR}_{i,t,U,r}$	$\displaystyle=$	$\displaystyle 1-\prod_{a=1}^{6}\left[\frac{1}{1+\exp(\beta_{a,r,t}+S_{i}+e_{i}+\delta_{i,t})}\right]^{z[a]}$		(11)
	$\displaystyle\mbox{U5MR}_{i,t,R,r}$	$\displaystyle=$	$\displaystyle 1-\prod_{a=1}^{6}\left[\frac{1}{1+\exp(\gamma_{a,r,t}+S_{i}+e_{i}+\delta_{i,t})}\right]^{z[a]},$		(12)

where $z[a]=1,11,12,12,12,12$, for the default choice of age bands. The aggregate risk in area $i$ and in year $t$ according to sampling frame $r$ is

$$p_{itr}=q_{itr}\times\mbox{U5MR}_{i,t,U,r}+(1-q_{itr})\times\mbox{U5MR}_{i,t,R,r},$$

(13)

where $q_{itr}$ and $1-q_{itr}$ are the proportions of the under-5 population in area $i$ that are urban and rural in year $t$ according to the classification of sampling frame $r$. The final aggregation over different sampling frames can be done using meta-analysis combination, so that,

$$\widehat{\mbox{U5MR}}_{it}=\mbox{expit}\left(\sum_{r}w_{itr}\times\mbox{logit}(p_{itr})\right),$$

where $w_{itr}=U_{itr}^{-1}/\sum_{r^{\prime}}U_{itr^{\prime}}^{-1}$ is the scaled inverse of $U_{itr}$, which is the posterior variance of $\mbox{logit}(\widehat{\mbox{U5MR}}_{it}^{(r)})$. Beyond point estimates, we obtain the full posterior of $\mbox{U5MR}_{it}$, and various summaries can be reported or mapped. The estimate constructed for U5MR is not relevant to any child, because that child would have to experience the hazards for each age group simultaneously in time period $t$, rather than moving through age groups over multiple time periods. Nevertheless, the resultant U5MR is a useful summary and the conventional measure that is used to inform on child mortality.

In this section, we present additional analysis for the simulated dataset in Example 2 of the main paper. We first load the package, data, and compute the direct estimates of U5MR following the example in the main paper.

library(SUMMER)library(stringr)library(dplyr)library(ggplot2)library(patchwork)data(DemoData)periods <- c("85-89", "90-94", "95-99", "00-04", "05-09", "10-14")directU5 <- getDirectList(births = DemoData, years = periods,                           regionVar = "region", timeVar = "time",                           clusterVar = "~clustid + id", ageVar = "age",                           weightsVar = "weights")

We now describe the fitting of the cluster level model. For simplicity, we first assume the survey was designed so that each of the four regions was a strata (thus no additional stratification within regions). The stratified analysis is left to the next section.

First, we calculate the number of person-months and number of deaths for each cluster, time period, and age group. Notice that we do not need to impute all the 0’s for combinations that do not exist in the data. We first create the data frame using the getCounts() function. The getCounts() function prepares the aggregated count dataset without modifying the original column names. For the model fitting functions to correctly identify the data columns, we rename the cluster ID and time period columns to be ‘cluster’ and ‘years’. The response variable is ‘Y’ and the binomial total is ‘total’.

counts.all <- NULLfor(i in 1:length(DemoData)){    vars <- c("clustid", "region", "time", "age")    counts <- getCounts(DemoData[[i]][, c(vars, "died")], variables = 'died',                         by = vars, drop=TRUE)    counts <- counts %>% mutate(cluster = clustid, years = time, Y=died)    counts$survey <- names(DemoData)[i]     counts.all <- rbind(counts.all, counts)}head(counts.all)

##   clustid  region  time age died total cluster years Y survey
## 1      36 central 85-89   0    0     1      36 85-89 0   1999
## 2      38 central 85-89   0    0     1      38 85-89 0   1999
## 3      91 central 85-89   0    0     1      91 85-89 0   1999
## 4     101 central 85-89   0    0     1     101 85-89 0   1999
## 5     128 central 85-89   0    0     1     128 85-89 0   1999
## 6     129 central 85-89   0    0     1     129 85-89 0   1999

With the created data frame, we fit the cluster-level model using the smoothCluster() function and obtain the estimates with the getSmoothed() function. Notice that here we need to specify the age groups (age.groups), the length of each age group (age.n) in months, and how the age groups are mapped to the temporal random effects (age.rw.group). In the default case, age.rw.group = c(1, 2, 3, 3, 3, 3) means the first two age groups each has its own temporal trend, the the following four age groups share the same temporal trend. We start with the default temporal model of random walk or order 2 on the 5-year periods in this dataset (with real data, we can use a finer temporal resolution). We add survey iid effects to the model as well using survey.effect = TRUE argument.

fit.bb.sim <- smoothCluster(data = counts.all, Amat = DemoMap$Amat,                 family = "betabinomial",                year.label = c(periods, "15-19"),                 age.group = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"),                age.n = c(1, 11, 12, 12, 12, 12),                age.time.group = c(1, 2, 3, 3, 3, 3),                time.model = "rw2",                st.time.model = "ar1",                pc.st.slope.u = 2, pc.st.slope.alpha = 0.1,                  survey.effect = TRUE)est.bb.sim <- getSmoothed(fit.bb.sim, nsim = 1000)

We visualize the results from the cluster level model in Figure 10. We also overlay the survey-specific direct estimates using the data.add argument.

plot(est.bb.sim$overall, plot.CI=TRUE,  data.add = directU5,             option.add = list(point = "mean", by = "surveyYears"),             color.add="steelblue") + facet_wrap(~region, ncol = 4) 

We now describe the fitting of the cluster-level model for U5MR taking into account the urban/rural stratification. For this simulated dataset, the strata variable is coded as region crossed by urban/rural status. For our analysis with urban/rural stratified model, we first construct a new strata variable that contains only the urban/rural status, i.e., the additional stratification within each region and computes the counts of person-months annd death similar to the previous case.

for(i in 1:length(DemoData)){  strata <- DemoData[[i]]$strata  DemoData[[i]]$strata[grep("urban", strata)] <- "urban"  DemoData[[i]]$strata[grep("rural", strata)] <- "rural"}counts.all <- NULLfor(i in 1:length(DemoData)){  vars <- c("clustid", "region", "strata", "time", "age")  counts <- getCounts(DemoData[[i]][, c(vars, "died")], variables = 'died',               by = vars, drop=TRUE)  counts <- counts %>% mutate(cluster = clustid, years = time, Y=died)  counts$survey <- names(DemoData)[i]   counts.all <- rbind(counts.all, counts)}

With the created data frame, we fit the cluster-level model using the smoothCluster function. The temporal main effects are defined for each stratum separately (specified by strata.time.effect = TRUE, so in total six random walks are used to model the main temporal effect.

fit.bb  <- smoothCluster(data = counts.all, Amat = DemoMap$Amat,             family = "betabinomial",            year.label = c(periods, "15-19"),             age.group = c("0", "1-11", "12-23", "24-35", "36-47", "48-59"),            age.n = c(1, 11, 12, 12, 12, 12),            age.time.group = c(1, 2, 3, 3, 3, 3),            time.model = "rw2",            st.time.model = "ar1",            pc.st.slope.u = 2, pc.st.slope.alpha = 0.1,              survey.effect = TRUE,             strata.time.effect = TRUE)

## ----------------------------------
## Cluster-level model
##   Main temporal model:        rw2
##   Number of time periods:     7
##   Spatial effect:             bym2
##   Number of regions:          4
##   Interaction temporal model: ar1
##   Interaction type:           4
##   Interaction random slopes:  yes
##   Number of age groups: 6
##   Stratification: yes
##   Number of age-specific fixed effect intercept per stratum: 6
##   Number of age-specific trends per stratum: 3
##   Strata-specific temporal trends: yes
##   Survey effect: yes
## ----------------------------------

est.bb <- getSmoothed(fit.bb, nsim = 1000, CI = 0.95, save.draws=TRUE)

## Starting posterior sampling...
## Cleaning up results...
## No strata weights has been supplied. Overall estimates are not calculated.

summary(fit.bb)

## ----------------------------------
## Cluster-level model
##   Main temporal model:        rw2
##   Number of time periods:     7
##   Spatial effect:             bym2
##   Number of regions:          4
##   Interaction temporal model: ar1
##   Interaction type:           4
##   Interaction random slopes:  yes
##   Number of age groups: 6
##   Stratification: yes
##   Number of age group fixed effect intercept per stratum: 6
##   Number of age-specific trends per stratum: 3
##   Strata-specific temporal trends: yes
##   Survey effect: yes
## ----------------------------------
## Fixed Effects
##                          mean   sd 0.025quant 0.5quant 0.97quant mode kld
## age.intercept0:rural     -3.0 0.14       -3.3     -3.0      -2.7 -3.0   0
## age.intercept1-11:rural  -4.7 0.11       -4.9     -4.7      -4.5 -4.7   0
## age.intercept12-23:rural -5.8 0.14       -6.1     -5.8      -5.5 -5.8   0
## age.intercept24-35:rural -6.6 0.20       -6.9     -6.6      -6.2 -6.6   0
## age.intercept36-47:rural -6.9 0.23       -7.3     -6.9      -6.4 -6.9   0
## age.intercept48-59:rural -7.3 0.29       -7.9     -7.3      -6.8 -7.3   0
## age.intercept0:urban     -2.7 0.13       -3.0     -2.7      -2.5 -2.7   0
## age.intercept1-11:urban  -5.0 0.13       -5.2     -5.0      -4.7 -5.0   0
## age.intercept12-23:urban -5.7 0.17       -6.1     -5.7      -5.4 -5.7   0
## age.intercept24-35:urban -7.1 0.29       -7.6     -7.1      -6.5 -7.1   0
## age.intercept36-47:urban -7.6 0.37       -8.3     -7.6      -6.9 -7.6   0
## age.intercept48-59:urban -8.0 0.46       -8.9     -8.0      -7.1 -8.0   0
##
## Slope fixed effect index:
## time.slope.group1: 0:rural
## time.slope.group2: 1-11:rural
## time.slope.group3: 12-23:rural, 24-35:rural, 36-47:rural, 48-59:rural
## time.slope.group4: 0:urban
## time.slope.group5: 1-11:urban
## time.slope.group6: 12-23:urban, 24-35:urban, 36-47:urban, 48-59:urban
## ----------------------------------
## Random Effects
##            Name             Model
## 1   time.struct         RW2 model
## 2 time.unstruct         IID model
## 3 region.struct        BYM2 model
## 4    region.int Besags ICAR model
## 5   st.slope.id         IID model
## 6     survey.id         IID model
## ----------------------------------
## Model hyperparameters
##                                                     mean       sd 0.025quant 0.5quant
## overdispersion for the betabinomial observations   0.002    0.001      0.001    0.002
## Precision for time.struct                         91.454  111.521     11.674   58.458
## Precision for time.unstruct                      831.467 3076.794     15.769  225.650
## Precision for region.struct                      252.997  671.939      5.264   88.691
## Phi for region.struct                              0.346    0.240      0.027    0.297
## Precision for region.int                         771.901 3485.913      7.931  168.610
## Group PACF1 for region.int                         0.893    0.203      0.253    0.970
## Precision for st.slope.id                         24.590   63.442      0.801    9.168
##                                                  0.97quant   mode
## overdispersion for the betabinomial observations     0.005  0.002
## Precision for time.struct                          349.819 27.909
## Precision for time.unstruct                       4750.233 35.610
## Precision for region.struct                       1367.096 11.338
## Phi for region.struct                                0.849  0.077
## Precision for region.int                          4531.658 15.399
## Group PACF1 for region.int                           0.999  1.000
## Precision for st.slope.id                          130.531  1.920
##                                        [,1]
## log marginal-likelihood (integration) -3455
## log marginal-likelihood (Gaussian)    -3449

The est.bb object above computes the U5MR estimates by urban/rural strata. In order to obtain the overall region-specific U5MR, we need additional information on the population fractions of each stratum within regions. For illustration purpose, here we simulate the population totals over the years and use these population totals to compute the population fractions later.

pop.base <- expand.grid(region = c("central", "eastern", "northern", "western"),                    strata = c("urban", "rural"))pop.base$population <- round(runif(dim(pop.base)[1], 1000, 20000))periods.all <- c(periods, "15-19")pop <- NULLfor(i in 1:length(periods.all)){  tmp <- pop.base  tmp$population <- pop.base$population + round(rnorm(dim(pop.base)[1], mean = 0, sd = 200))  tmp$years <- periods.all[i]  pop <- rbind(pop, tmp)}head(pop)

##     region strata population years
## 1  central  urban      16529 85-89
## 2  eastern  urban      10630 85-89
## 3 northern  urban      14664 85-89
## 4  western  urban      11470 85-89
## 5  central  rural       5836 85-89
## 6  eastern  rural       5479 85-89

In order to compute the aggregated estimates, we need the proportion of urban/rural populations within each region in each time period, as computed below in the weight.strata object.

weight.strata <- expand.grid(region = c("central", "eastern", "northern", "western"),                    years = periods.all)weight.strata$urban <- weight.strata$rural <- NAfor(i in 1:dim(weight.strata)[1]){  which.u <- which(pop$region == weight.strata$region[i] &                   pop$years == weight.strata$years[i] &                   pop$strata == "urban")  which.r <- which(pop$region == weight.strata$region[i] &                   pop$years == weight.strata$years[i] &                   pop$strata == "rural")  weight.strata[i, "urban"] <- pop$population[which.u] /                        (pop$population[which.u] + pop$population[which.r])  weight.strata[i, "rural"] <- 1 - weight.strata[i, "urban"]}head(weight.strata)

##     region years rural urban
## 1  central 85-89  0.26  0.74
## 2  eastern 85-89  0.34  0.66
## 3 northern 85-89  0.53  0.47
## 4  western 85-89  0.34  0.66
## 5  central 90-94  0.26  0.74
## 6  eastern 90-94  0.36  0.64

Now we can recompute the smoothed estimates with the population fractions.

est.bb <- getSmoothed(fit.bb, nsim = 1000, CI = 0.95, save.draws = TRUE,                       weight.strata = weight.strata)head(est.bb$overall)

##    region years time area variance median mean upper lower rural urban is.yearly
## 5 central 85-89    1    1  0.00149   0.22 0.22  0.30  0.16  0.26  0.74     FALSE
## 6 central 90-94    2    1  0.00060   0.20 0.20  0.25  0.15  0.26  0.74     FALSE
## 7 central 95-99    3    1  0.00033   0.18 0.18  0.22  0.15  0.25  0.75     FALSE
## 1 central 00-04    4    1  0.00026   0.17 0.17  0.21  0.14  0.25  0.75     FALSE
## 2 central 05-09    5    1  0.00027   0.16 0.16  0.20  0.13  0.26  0.74     FALSE
## 3 central 10-14    6    1  0.00043   0.15 0.15  0.19  0.11  0.26  0.74     FALSE
##   years.num
## 5        NA
## 6        NA
## 7        NA
## 1        NA
## 2        NA
## 3        NA

We can compare the stratum-specific and aggregated U5MR estimates now.

g1 <- plot(est.bb$stratified, plot.CI = TRUE) + facet_wrap(~strata) + ylim(0, 0.5)g2 <- plot(est.bb$overall, plot.CI = TRUE)  + ylim(0, 0.5) + ggtitle("Aggregated estimates")g1 + g2

In this example, we use the 2010 and 2015–2016 Malawi DHS surveys. The DHS website (https://dhsprogram.com/data/available-datasets.cfm?ctryid=24) provides links to download all the surveys with registration. Once access is approved, we can use the rdhs package to load data directly from the DHS API (Watson and Eaton, 2019). We first download both the birth records (BR) and GPS data (GE) for these two surveys. Notice the following codes can only be run after registration with DHS and set up the credentials using the rdhs package. Additional information on this step can be found in the documentation of the rdhs package.

library(rdhs)sv <- dhs_surveys(countryIds = "MW", surveyType = "DHS",                   surveyYearStart = 2010)BR <- dhs_datasets(surveyIds = sv$SurveyId, fileFormat = "Flat",                    fileType = "BR")BRfiles <- get_datasets(BR$FileName, reformat=TRUE)GPS <- dhs_datasets(surveyIds = sv$SurveyId, fileFormat = "Flat",                     fileType = "GE")GPSfiles <- get_datasets(GPS$FileName, reformat=TRUE)

The downloaded data can then be loaded by the returned file paths.

Surv2010 <- readRDS(BRfiles[[1]])Surv2015 <- readRDS(BRfiles[[2]])DHS2010.geo <- readRDS(GPSfiles[[1]])DHS2015.geo <- readRDS(GPSfiles[[2]])

load("data/downloaded.rda")

We then use the getBirths() function to process the raw birth history into person-month format. We label the person-month records with 6 age groups specified by month.cut. In this example, instead of using 5 year periods, we work directly with yearly estimates specified by year.cut.

DHS2010 <- getBirths(data = Surv2010,                  month.cut = c(1, 12, 24, 36, 48, 60),                  year.cut = seq(2000, 2020, by = 1), strata = "v022")DHS2015 <- getBirths(data = Surv2015,                  month.cut = c(1, 12, 24, 36, 48, 60),                  year.cut = seq(2000, 2020, by = 1), strata = "v022")

We perform similar processing steps for the 2015–2016 DHS survey birth records. Since only a small fraction of observations are available in 2016, we remove the partial year observations in this survey.

DHS2015 <- subset(DHS2015, time != 2016)

The DHS dataset does not usually have sufficient spatial resolution information in the birth records file, so we need to use the GPS datasets of cluster locations to assign records to the admin-2 areas. This process can be highly survey-specific and require more extensive data manipulation. In the 2010 DHS dataset, the GPS file contains the cluster ID (DHSCLUST), urbanicity indicator (URBAN_RURA), and DHSCLUST variable in the format of “admin-2_urbanrural”, with some spelling and capitalization differences as in the polygon file. From the GPS dataset, we processed a list of clusters and their corresponding admin-2 areas.

cluster.list <- data.frame(DHS2010.geo) %>%         distinct(DHSCLUST, DHSREGNA, URBAN_RURA) %>%         mutate(admin2 = gsub(" - rural", "", DHSREGNA)) %>%         mutate(admin2 = gsub(" - urban", "", admin2)) %>%         mutate(admin2 = recode(admin2, "nkhatabay" = "nkhata bay")) %>%         mutate(admin2 = recode(admin2, "nkhota kota" = "nkhotakota")) %>%         mutate(admin2 =  str_to_title(admin2)) %>%         mutate(urban = ifelse(URBAN_RURA=="U", 1, 0))  %>%         select(v001 = DHSCLUST, admin2, urban) head(cluster.list, n=3)     

##   v001  admin2 urban
## 1    1   Dedza     0
## 2    2  Balaka     0
## 3    3 Mchinji     0

We then merge this list to the main person-month file, which adds the admin2 and urban columns to the data frame.

DHS2010 <- DHS2010 %>% left_join(cluster.list)head(DHS2010)

##    dob survey_year died id.new          caseid v001 v002 v004    v005 v021 v022    v023
## 1 1316          NA    0      1         1  1  2    1    1    1 1902867    1   12 central
## 2 1316          NA    0      1         1  1  2    1    1    1 1902867    1   12 central
## 3 1316          NA    0      1         1  1  2    1    1    1 1902867    1   12 central
## 4 1316          NA    0      1         1  1  2    1    1    1 1902867    1   12 central
## 5 1316          NA    0      1         1  1  2    1    1    1 1902867    1   12 central
## 6 1316          NA    0      1         1  1  2    1    1    1 1902867    1   12 central
##      v024  v025    v139 bidx agemonth obsStart obsStop obsmonth year  age time strata
## 1 central rural central    1        0     1316    1316     1316  109    0 2009     12
## 2 central rural central    1        1     1317    1317     1317  109 1-11 2009     12
## 3 central rural central    1        2     1318    1318     1318  109 1-11 2009     12
## 4 central rural central    1        3     1319    1319     1319  109 1-11 2009     12
## 5 central rural central    1        4     1320    1320     1320  109 1-11 2009     12
## 6 central rural central    1        5     1321    1321     1321  110 1-11 2010     12
##   admin2 urban
## 1  Dedza     0
## 2  Dedza     0
## 3  Dedza     0
## 4  Dedza     0
## 5  Dedza     0
## 6  Dedza     0

We perform similar processing steps for the 2015–2016 DHS GPS data. To illustrate the idiosyncratic data processing steps required for each survey, the DHSCLUST variable in this dataset is in the form of “admin-1_urbanrural” and does not allow us to parse the admin-2 area names. The strata variable v022 in the birth records, however, is in the “admin-2_urbanrural” format. Therefore, we first merge the v022 variables to the GPS data and proceed in a similar fashion to the previous example.

cluster.list <- data.frame(DHS2015.geo) %>%         mutate(v001 = DHSCLUST) %>%         left_join(DHS2015, by = "v001") %>%         distinct(v001, v022, URBAN_RURA) %>%         mutate(admin2 = gsub(" - rural", "", v022)) %>%         mutate(admin2 = gsub(" - urban", "", admin2)) %>%         mutate(admin2 = recode(admin2, "nkhatabay" = "nkhata bay")) %>%         mutate(admin2 = recode(admin2, "nkhota kota" = "nkhotakota")) %>%         mutate(admin2 =  str_to_title(admin2)) %>%         mutate(urban = ifelse(URBAN_RURA=="U", 1, 0))  %>%         select(v001, admin2, urban) DHS2015 <- DHS2015 %>% left_join(cluster.list)head(DHS2015)

##    dob survey_year died id.new          caseid v001 v002 v004   v005 v021            v022
## 1 1334          NA    0      1        1  12  2    1   12    1 118748    1 ntchisi - urban
## 2 1334          NA    0      1        1  12  2    1   12    1 118748    1 ntchisi - urban
## 3 1334          NA    0      1        1  12  2    1   12    1 118748    1 ntchisi - urban
## 4 1334          NA    0      1        1  12  2    1   12    1 118748    1 ntchisi - urban
## 5 1334          NA    0      1        1  12  2    1   12    1 118748    1 ntchisi - urban
## 6 1334          NA    0      1        1  12  2    1   12    1 118748    1 ntchisi - urban
##              v023           v024  v025           v139 bidx agemonth obsStart obsStop
## 1 ntchisi - urban central region urban central region    1        0     1334    1334
## 2 ntchisi - urban central region urban central region    1        1     1335    1335
## 3 ntchisi - urban central region urban central region    1        2     1336    1336
## 4 ntchisi - urban central region urban central region    1        3     1337    1337
## 5 ntchisi - urban central region urban central region    1        4     1338    1338
## 6 ntchisi - urban central region urban central region    1        5     1339    1339
##   obsmonth year  age time          strata  admin2 urban
## 1     1334  111    0 2011 ntchisi - urban Ntchisi     1
## 2     1335  111 1-11 2011 ntchisi - urban Ntchisi     1
## 3     1336  111 1-11 2011 ntchisi - urban Ntchisi     1
## 4     1337  111 1-11 2011 ntchisi - urban Ntchisi     1
## 5     1338  111 1-11 2011 ntchisi - urban Ntchisi     1
## 6     1339  111 1-11 2011 ntchisi - urban Ntchisi     1

Then we use the getCounts() function to obtain the count data format input and stack the two DHS survey data into a single data frame.

vars <- c("v001", "v025", "admin2", "time", "age", "v005")dat1 <- getCounts(DHS2010[, c(vars, "died")], variables = 'died',                   by = vars, drop=TRUE)dat1$survey = "DHS2010"dat2 <- getCounts(DHS2015[, c(vars, "died")], variables = 'died',                   by = vars, drop=TRUE)dat2$survey = "DHS2015"DHS.counts <- rbind(dat1, dat2) %>%                mutate(cluster = v001, strata = v025, region = admin2,                        years = time, Y = died)

The analysis in the main paper is carried out only on the 2015–2016 dataset, which is obtained using the same procedure.

We use getDirectList to obtain direct estimates from both surveys and combine into the ‘meta-analysis’ estimator using aggregateSurvey. Notice that the survey variable in the returned data frame contains a numeric index of each survey in the list, and the surveyYears contains the survey names we assigned in the input list.

direct <- getDirectList(births = list(DHS2010=DHS2010, DHS2015=DHS2015),         years = 2000:2019, regionVar = "admin2", timeVar = "time",         clusterVar = "~v001 + v002", ageVar = "age", weightsVar = "v005")direct.comb <- aggregateSurvey(direct)

When additional information is available to adjust the direct estimates from the surveys, we use the methods described in Li et al. (2019). We can perform the ratio adjustment to the direct estimates using the getAdjusted() function. For the two surveys in Malawi, the calculated HIV adjustment ratios as described in Walker et al. (2012) are stored in MalawiData$HIV.yearly. In order to get the correct uncertainty bounds, we also need to specify the columns corresponding to the unadjusted uncertainty bounds, the lower and upper columns, which are on the probability scale in this case.

data(MalawiData)direct.2010 <- subset(direct, survey == 1)direct.2010.hiv <- getAdjusted(data = direct.2010,                 ratio = subset(MalawiData$HIV.yearly, survey == "DHS2010"),                 logit.lower = NA, logit.upper = NA,                 prob.lower = "lower", prob.upper = "upper")direct.2015 <- subset(direct, survey == 2)direct.2015.hiv <- getAdjusted(data = direct.2015,                 ratio = subset(MalawiData$HIV.yearly, survey == "DHS2015"),                 logit.lower = NA, logit.upper = NA,                 prob.lower = "lower", prob.upper = "upper")

Finally, we combine the direct estimates into direct.comb.hiv.

direct.hiv <- rbind(direct.2010.hiv, direct.2015.hiv)direct.comb.hiv <- aggregateSurvey(direct.hiv)

We now fit a national Fay-Harriot model with the calculated direct estimates using smoothDirect() and getSmoothed() functions.

fit.national.unadj <- smoothDirect(data = direct.comb, Amat = NULL,                         year.label = 2000:2019, year.range = c(2000, 2019),                         time.model = "rw2", m = 1)est.unadj <- getSmoothed(fit.national.unadj)

For comparison, we smooth both the unadjusted direct estimates and the direct estimates with HIV adjustments.

fit.national.hiv <- smoothDirect(data = direct.comb.hiv,  Amat = NULL,                         year.label = 2000:2019, year.range = c(2000, 2019),                         time.model = "rw2", m = 1)est.hiv <- getSmoothed(fit.national.hiv)

In addition, we also demonstrate the benchmarking procedure described in Li et al. (2019), where we first fit a smoothing model and then benchmark the results with the UN IGME estimates – the latter are based on more extensive data (Alkema and New, 2014) . We first calculate the adjustment ratio compared to the 2019 UN IGME estimates.

UN <- MalawiData$IGME2019UN.est <- UN$mean[match(2000:2019, UN$years)]Smooth.est <- est.hiv$median[match(2000:2019, est.hiv$years)]UN.adj <- data.frame(years = 2000:2019, ratio = Smooth.est / UN.est)head(UN.adj, n = 3)

##   years ratio
## 1  2000   1.1
## 2  2001   1.2
## 3  2002   1.1

We then fit the smoothing model on the benchmarked direct estimates.

direct.comb.benchmark <- getAdjusted(data = direct.comb.hiv, ratio = UN.adj,                                logit.lower = NA, logit.upper = NA,                                 prob.lower = "lower", prob.upper = "upper")fit.benchmark <- smoothDirect(data = direct.comb.benchmark, Amat = NULL,                          year.label = 2000:2019, year.range = c(2000, 2019),                          time.model = "rw2",  m = 1)est.benchmark <- getSmoothed(fit.benchmark)

We compare the different Fay-Harriot estimates in Figure 13. Compared to the raw estimates, HIV adjustments lead to higher estimates in the earlier years. The benchmarking step produces a national trend that follows the same trajectory as the UN IGME estimates.

g1 <- plot(est.unadj, is.subnational=FALSE, proj.year = 2016) +            ggtitle("Unadjusted") + ylim(c(0, 0.22))g2 <- plot(est.hiv, is.subnational=FALSE, proj.year = 2016) +            ggtitle("With HIV adjustment") + ylim(c(0, 0.22))g3 <- plot(est.benchmark, is.subnational=FALSE, proj.year = 2016, data.add = UN,     option.add = list(point = "mean"), label.add = "UN", color.add = "red") +    ggtitle("Benchmarked to UN IGME") + ylim(c(0, 0.22))g1 + g2 + g3 

Finally, the subnational models can be fit in the same manner. Since the direct estimates are already on the yearly scale, we do not need to perform the transformation from periods to years. So we set is.yearly to FALSE to fit the period model directly.

fit.smooth.direct <- smoothDirect(data = direct.comb.benchmark,                             Amat = MalawiGraph,                             year.label = 2000:2019, year.range = c(2000, 2019),                             time.model = "rw2", m = 1,                             type.st = 4, pc.alpha = 0.05, pc.u = 1)est.smooth.direct <- getSmoothed(fit.smooth.direct)

We now describe the fitting of the cluster level model using the two DHS surveys. With the created data frame, we first fit the national model with survey-year-specific HIV adjustment factors specified using bias.adj and bias.adj.by arguments. The adjustments are performed as offsets in the likelihood as described before. We also add a sum-to-zero survey effect term.

fit.bb.nat <- smoothCluster(data = DHS.counts, Amat = NULL,                     family = "betabinomial", year.label = 2000:2019,                     time.model = "rw2",                     bias.adj = MalawiData$HIV.yearly,                     bias.adj.by = c("years", "survey"),                    survey.effect = TRUE)

The getSmoothed function then produces estimates from the fitted cluster-level model. Similar to before, we take nsim draws of the posterior distribution to calculate the summaries of the U5MR estimates.

est.bb.nat <- getSmoothed(fit.bb.nat, nsim = 1000, save.draws = TRUE) 

In this example, no strata weights are provided and thus the overall estimates are empty. Given a data frame of strata proportions, we can rerun the getSmoothed function to re-aggregate the stratified estimates. The save.draws argument in the getSmoothed call allows the raw posterior draws to be returned as part of the output object. This can be helpful in such situations, as posterior draws already computed can be inserted into new getSmoothed calls using the draws argument to avoid resampling again.

Figure 14 shows the national estimates of U5MR in Malawi for urban and rural stratum respectively using the cluster-level model.

plot(est.bb.nat$stratified, is.subnational=FALSE, proj.year = 2016) + facet_wrap(~strata) 

We can also fit the subnational model in the same fashion. Additional analysis can be similarly carried out as in the previous example with simulated data. For the cluster-level model, benchmarking to national estimates is also implemented in the package using the procedure described in [@okonek2022computationally] with additional information on population fractions by region. We refer readers to the package vignette for more details.

fit.bb <- smoothCluster(data = DHS.counts, Amat = MalawiGraph,                     family = "betabinomial",                     year.label = 2000:2019,                     time.model = "rw2", st.time.model = "ar1",                    pc.st.slope.u = 2,                     pc.st.slope.alpha = 0.1,                    bias.adj = MalawiData$HIV.yearly,                     bias.adj.by = c("years", "survey"),                    survey.effect = TRUE,                     strata.time.effect = TRUE)est.bb <-  getSmoothed(fit.bb, nsim = 1000, save.draws = TRUE) 

The U5MR by strata can be visualized directly. Notice that in order to produce the overall estimates by region and time, additional information on population fractions in urban/rural is necessary, similar to the analysis conducted in the main paper using simulated data. For more details on obtaining such factions, we refer readers to [@fuglstad2021two] and [@wu2021spatial].

plot(est.bb$stratified) + facet_wrap(~strata)