# A spatial variance-smoothing area level model for small area estimation of
demographic rates

Peter A. Gao and Jon Wakefield

###### Abstract

Accurate estimates of subnational health and demographic indicators are
critical for informing health policy decisions. Many countries collect
relevant data using complex household surveys, but when data are limited,
direct survey weighted estimates of small area proportions may be unreliable.
Area level models treating these direct estimates as response data can improve
precision but often require known sampling variances of the direct estimators
for all areas. In practice, the sampling variances are typically estimated, so
standard approaches do not account for a key source of uncertainty. In order
to account for variability in the estimated sampling variances, we propose a
hierarchical Bayesian spatial area level model that smooths both the estimated
means and sampling variances to produce point and interval estimates of small
area proportions. Our model explicitly targets estimation of small area
proportions rather than means of continuous variables and we consider examples
of both moderate and low prevalence events. We demonstrate the performance of
our approach via simulation and application to vaccination coverage and HIV
prevalence data from the Demographic and Health Surveys.

Key words: Bayesian statistics, small area estimation, area level model,
spatial statistics, survey statistics.

##  1 Introduction

Subnational estimates of health and demographic indicators such as vaccination
rates or neonatal mortality rates are critical for highlighting disparities
and guiding policy interventions. Data on key health outcomes are often
collected using national household surveys designed to produce reliable direct
weighted estimators for national and regional demographic rates of interest.
However, often estimates are desired for smaller subregions for which direct
weighted estimation may be insufficiently precise. When data are limited,
model-based small area estimation methods can often improve estimates by
sharing information between areas. Among model-based methods, unit level
approaches are popular in demographic research as they explicitly model
individual survey responses and can incorporate available individual level
covariate information. However, if care is not taken to account for
informative sampling or other features of the survey design, estimates from
unit level models may be biased or improperly calibrated [1, 2]. On the other
hand, area level model-based approaches model direct survey weighted estimates
as noisy response data, using random effects to produce smoothed small area
estimates. The Fay-Herriot area level model [3] assumes that for each area,
the direct estimator $\widehat{p}_{a}$ is available and can be modeled using
the Gaussian distribution centered around the finite population mean $p_{a}$:

$\widehat{p}_{a}\mid p_{a},V_{a}\sim N(p_{a},V_{a})$ (1)

where $V_{a}$ denotes sampling variances, which are typically assumed to be
known. This sampling model is combined with a so-called linking model for the
finite population means:

$p_{a}\mid\boldsymbol{\beta},\sigma_{u}^{2}\sim
N(\mathbf{x}_{a}^{T}\boldsymbol{\beta},\sigma_{u}^{2})$ (2)

where $\mathbf{x}_{a}$ denotes a vector of area level covariate values for
area $a$ and $\boldsymbol{\beta}$ denotes the corresponding coefficients. The
linking model variance $\sigma_{u}^{2}$ controls the magnitude of deviations
from the mean model. Since the direct estimators $\widehat{p}_{a}$ account for
the survey design via the use of survey weights, area level methods are less
sensitive to the effects of informative sampling and other design features
than unit level methods. Under certain regularity conditions, the resulting
estimators are design consistent; for a review, see Rao and Molina (2015) [4].

In practice, the variances $V_{a}$ are usually estimated using sample-based
estimators $\widehat{V}_{a}$, but the standard Fay-Herriot model does not
account for uncertainty in $\widehat{V}_{a}$. This is a well known problem in
the small area estimation literature [5, 6, 7] and has motivated a number of
proposed extensions of Fay-Herriot that incorporate variance modeling.
However, existing approaches often rely on the availability of informative
area level covariates for modeling the variance of survey estimators and do
not account for uncertainty in the modeled variance estimates, simply treating
them as known and using them to replace the direct variance estimates in the
standard Fay-Herriot model. To account for this source of uncertainty, we
propose a fully Bayesian area level model for small area proportions that
jointly models the direct estimators and associated variance estimators. Our
approach assumes that area level random effects are spatially correlated and
induces spatial smoothing of both means and variances. While many existing
extensions to Fay-Herriot focus on estimation of means of continuous
variables, our method is designed for estimation of small area proportions. In
health and demographic applications, such proportions are the most common
targets of inference. In simulations, we find that our proposed method
produces interval estimates that may help correct for undercoverage observed
for the standard Fay-Herriot approach.

Below we expand upon our proposed model and draw connections to existing
approaches. In Section 2, we describe two examples of estimation problems
involving estimation of subnational demographic rates to provide motivation
for our new approach. Section 3 reviews existing area level models and
discusses recent efforts to incorporate variance smoothing for small area
estimation. In Section 4, we outline our new spatial variance-smoothing area
level model for estimation of small area proportions. We compare our approach
with other area level methods via simulation in Section 5 and by application
to data from the Demographic and Health Surveys in Section 6. Finally, we
compare our method with existing approaches and identify potential directions
for future research in Section 7.

##  2 Motivation

We consider two motivating examples of estimating subnational demographic
rates using data from the Demographic and Health Surveys (DHS) Program. In the
first, we use 2018 Nigeria DHS data to estimate regional vaccination coverage
rates for the first dose of measles-containing-vaccine (MCV1) among children
aged 12–23 months [8]. In the second, we use 2015-16 Malawi DHS data to
estimate HIV prevalence among women aged 15–49 [9]. Figure 1 provides maps of
direct survey-weighted estimators for both indicators. The measles vaccination
example represents an estimation problem where the estimated area level
proportions have a large spread and are generally located away from zero or
one; in the HIV prevalence example, the direct estimates exhibit less
variability and are on average closer to zero.

Figure 1: Direct weighted estimates of vaccination coverage rate for first
dose of measles-containing-vaccine (MCV1) among children aged 12–23 months in
Nigeria, 2018 (left) and HIV prevalence rate for women aged 15-49 in Malawi,
2015-2016 (right)

In Nigeria and Malawi, the DHS Program uses a stratified two-stage cluster
sampling design. Countries are divided into administrative regions which are
further partitioned into urban and rural areas. The sampling strata are
defined by crossing these regions with urban/rural status. In Nigeria, the
divisions used for defining strata are called Admin-1 regions; in Malawi, they
are called Admin-2 regions. Each stratum is divided into collections of
households called enumeration areas (EAs) or clusters. The first stage of
sampling selects a pre-specified number of EAs in each stratum with
probability proportional to the number of households in the EA. The second
stage of sampling selects a fixed number of households in each sampled EA.

The 2018 Nigeria DHS collected data on measles vaccination status for children
in sampled households based on vaccination cards or caregiver recall. In
Nigeria, the Admin-1 regions consist of 36 states and the Federal Capital
Territory of Abuja. For our analysis, we adhere to the Database of Global
Administrative Areas (GADM) boundaries
(https://gadm.org/download_country_v3.html). The sampling frame used for the
Nigeria DHS was based on a 2006 national census which identified 664,999 EAs.
Data were successfully collected in 1389 EAs, but due to security issues
during the survey, a number of EAs were dropped. In particular, estimates for
the Admin-1 area of Borno may have been affected (see Appendix A.3 of [8]).

The 2015-16 Malawi DHS used voluntary finger prick blood sampling to collect
data on HIV prevalence. We desire estimates of HIV prevalence for each of
Malawi’s 28 districts, also referred to as Admin-2 areas. For this survey, the
sampling frame was obtained from a 2008 census which identified 12,558 EAs
distributed between 56 strata. Ultimately, data were collected from 827 EAs,
from which a total of 8,497 women aged 15-49 were eligible for HIV testing.
Ultimately, 93% of eligible women were tested, but the HIV test results were
anonymized, with volunteers not informed of their results and instead
receiving access to educational materials and free counseling and testing [9].

For both Nigeria and Malawi, the DHS provides GPS coordinates for nearly all
EAs, but the locations have been adjusted to maintain privacy by adding small
distances at random. Figure 2 provides maps of the small area boundaries and
sampled EA locations in Nigeria and Malawi. In Malawi, since the island region
of Likoma is disconnected from the mainland, we omit its data from our
analysis.

Figure 2: Small area boundaries and sampled enumeration area locations for
Nigeria (left) and Malawi (right)

##  3 Existing Work

###  3.1 Variance Smoothing for Continuous Response

Standard area level models, including the Fay-Herriot model, do not account
for uncertainty in the estimates of sampling variance $\widehat{V}_{a}$. As
such, the basic Fay-Herriot model for continuous response data has been
extended to account for unknown $V_{a}$ in a number of ways as introduced by
Kleffe and Rao (1992) and Arora and Lahiri (1997)[10, 5] and reviewed by Bell
(2008) [7]. Below, we review extensions of area level models. Following Rivest
and Vandal (2002) [6] and Wang and Fuller (2003) [11], You and Chapman (2006)
[12] assume the following sampling model for the variance estimators
$\widehat{V}_{a}$:

$\widehat{V}_{a}\mid V_{a}\sim\frac{V_{a}}{d_{a}}\chi^{2}_{d_{a}}$ (3)

where $d_{a}$ denotes the degrees of freedom for area $a$. In addition,
$\widehat{V}_{a}$ are assumed to be independent of the mean estimators
$p_{a}$. If the response values for area $a$ were independently and
identically distributed Gaussian random variables, the above model (3) would
hold for the variance estimator $\widehat{V}_{a}=s_{a}^{2}/n_{a}$ with
$d_{a}=n_{a}-1$, where $n_{a}$ denotes the sample size for area $a$. When
responses are sampled at random with replacement within areas, such an
assumption may be appropriate, but for complicated sampling schemes, different
values of $d_{a}$ or even alternative models may be necessary. You and Chapman
(2006) [12] adopted a hierarchical Bayesian approach and placed inverse Gamma
priors on the variance parameters $\sigma_{u}^{2}\sim IG(r_{0},s_{0})$ and
$V_{a}\sim IG(r_{a},s_{a})$, with $r_{a},s_{a}$ chosen to be small for all
areas $a=1,\ldots A$. Notably, they allow the prior for $V_{a}$ to vary across
areas, which makes the $V_{a}$ values independent across areas. Maiti et al.
(2014) [13] assume the same variance sampling model (3), but adopt an
empirical Bayes approach, setting the prior $\sigma_{a}^{2}\sim IG(r,s)$ and
estimating $\\{r,s,\sigma_{u}^{2}\\}$ via maximum likelihood.

In addition to modeling $\widehat{V}_{a}$, Hwang et al. (2009) [14] and Dass
et al. (2012) [15] noted that assuming a common prior for $V_{a}$ for all
areas $a$ could induce shrinkage in the resulting variance estimates and
produce improved interval estimates for the parameters of interest. In this
vein, Sugasawa et al. (2017)[16] explore different priors for the sampling
variances $V_{a}$, adopting a fully Bayesian approach to estimation.
Alternatively, Polettini (2017) [17] induces shrinkage for the sampling
variance estimates using a semiparametric Dirichlet process model with random
variances.

###  3.2 Variance Smoothing for Binary Response

When response values are binary and the target of estimation is a small area
proportion, it may be helpful to account for the mean-variance relationship
observed in binary response data. Generalized variance functions (GVFs), which
model the functional relationship between the expectation and variance of a
survey estimator, can be used as an alternative to linearization-based
approximations or resampling methods for estimating $V_{a}$. If the model used
is appropriate, the resulting modeled variance estimates could improve upon
the direct variance estimates in terms of precision. An introduction to GVFs
is provided in Chapter 7 of Wolter (2007) [18].

For small area estimation of proportions, several GVF-like approaches to
variance estimation have been previously proposed based on treating the
responses like binomial data. Liu et al. (2014) [19] assume the following
model for $V_{a}$:

$V_{a}=\frac{p_{a}(1-p_{a})}{n_{a}}DEFF_{a}$ (4)

where $DEFF_{a}$ denotes the design effect, defined as the ratio of the
variance of $p_{a}$ under the implemented survey design to the the variance of
$p_{a}$ under simple random sampling. As described in their paper, Liu et al.
estimate design effects using available information on sample sizes and survey
weights and treat them as known. Model-based estimates of $V_{a}$ can be
produced by replacing the unknown $p_{a}$ values above with their direct
estimators. Hawala and Lahiri (2018) [20] propose a similar GVF for count
data. Maples (2016) [21] similarly proposes a GVF for producing variance
estimates based on estimating the design effect using additional information
about any unequal weighting or clustering in the sampling procedure. Franco
and Bell (2013) adopt a different strategy using a GVF to compute an effective
sample size for each area of interest, which they use to fit a binomial model
[22].

Mohadjer et al. (2012) [23] similarly use a GVF to produce variance estimates
for use in an area level model, assuming the following model

$\log(V_{a}/p_{a}^{2})=\eta_{0}+\eta_{1}\log(\widetilde{p}_{a})+\eta_{2}\log(1-\widetilde{p}_{a})+\eta_{3}\log(n_{a})+\varepsilon_{a}$
(5)

where $\epsilon_{a}\sim N(0,\sigma_{\epsilon}^{2})$ and $\widetilde{p}_{a}$
denotes a predictor of $p_{a}$ based on a model dependent solely on auxiliary
covariate information and not explicitly on any direct estimates.

The GVF approaches described thus far treat the resulting variance estimates
as known, so the resulting Fay-Herriot estimates do not account for
uncertainty in the variance model. Maples et al. (2009) [24] address this by
combining a GVF with a sampling model for the direct variance estimates. In
particular, they assume Model (3) holds for the direct variance estimates
$\widehat{V}_{a}$ and then propose the following linking model for $V_{a}$:

$V_{a}\mid\alpha,\boldsymbol{\gamma}\sim
IG(\alpha+1,\alpha\exp(\mathbf{z}_{a}^{T}\boldsymbol{\gamma}))$ (6)

where $\alpha$ controls the precision of the variance linking model,
$\mathbf{z}_{a}$ are area level covariates and $\boldsymbol{\gamma}$ are
corresponding coefficients estimated using an empirical Bayes approach. Maples
et al. outline a procedure for using bootstrap sampling to estimate effective
sample size for each area, which informs their choice for the degrees of
freedom $d_{a}$ in the variance sampling model. They show that this model
produces smoothed variance estimates that could help to correct
underestimation by the direct variance estimators.

###  3.3 Alternative sampling and linking models

The linking and sampling mean models in the Fay-Herriot approach assume
responses are continuous, but since $p_{a}$ are bounded between 0 and 1, it
may be inappropriate to treat $\widehat{p}_{a}$ as Gaussian, especially when
$p_{a}$ is close to 0 or 1 and when $V_{a}$ is large. In the health and
demography setting, Mercer et al. [25] apply a logit transformation to direct
estimates of mortality rates before fitting a Fay-Herriot-type model. You and
Rao (2002) [19] proposed to use unmatched sampling and linking models,
combining the sampling model given by Equation (1) with an alternative linking
model that transforms the finite population parameters of interest $p_{a}$ to
make a Gaussian approximation more appropriate. As an example, Liu et al.
(2014) [19] considered the following logit-normal linking model:

$\mathrm{logit}(p_{a})\mid\boldsymbol{\beta},\sigma_{u}^{2}\sim
N(\mathbf{x}_{a}^{T}\boldsymbol{\beta},\sigma_{u}^{2})$ (7)

Mohadjer et al. (2012) [23] apply this model to estimation of adult literacy
rates. Liu et al. also consider alternative models including a beta-logistic
model combining a beta sampling model with the above logistic linking model,
which accounts for the limited range of $\widehat{p}_{a}$ but will not reflect
its true sampling distribution. Franco and Bell (2013) [22] and Chen et al.
(2014) [26] consider binomial sampling models, treating observed area level
counts as being drawn from a binomial distribution with size parameter given
by some measure of effective sample size. As an alternative to unmatched
sampling and linking models, Mercer et al. (2015) [25] describe an approach
that uses Gaussian sampling and linking models to model
$\mathrm{logit}(\widehat{p}_{a})$ and $\mathrm{logit}(p_{a})$.

###  3.4 Spatial area level models

The Fay-Herriot model has been extended to incorporate spatial and
spatiotemporal random effects [27, 28]. These models leverage similarities
between areas that are close in space or time, producing smoothed estimates
and accounting for potential spatial or spatiotemporal structure in the
response data. Chung et al. (2020) [29] noted that when informative area level
covariate information is unavailable but responses are spatially correlated,
using spatial random effects models may be especially effective. In LMICs,
covariate information is often limited and spatial random effects are often
used in area level and unit level modeling. As an example, Mercer et al.
(2015) use an area level model with spatiotemporal random effects to estimate
child mortality rates in Tanzania [25].

##  4 Methods

We assume that for all $a=1,\ldots,A$, we have direct estimates of area level
proportions $\widehat{p}_{a}$ and corresponding variance estimates
$\widehat{V}_{a}$. We propose a Bayesian joint model for the full data
$(\widehat{\mathbf{p}},\widehat{\mathbf{V}})$ that induces spatial smoothing
for both the proportion and variance estimates. Our approach uses two sets of
unmatched models, one for the estimated proportions $\widehat{\mathbf{p}}$ and
one for the variance estimates $\widehat{\mathbf{V}}$, with these models being
linked through the use of a generalized variance function. We use a spatial
linking model for the proportions that induces spatial smoothing for both the
proportions and the estimated variances.

###  4.1 Mean model

For modeling the direct estimates $\widehat{p}_{a}$, we use unmatched sampling
and linking models, combining a Gaussian sampling model with a spatial logit-
Gaussian linking model:

$\widehat{p}_{a}\mid p_{a},V_{a}\stackrel{{\scriptstyle
ind}}{{\sim}}N(p_{a},V_{a}).$ (8)
$\mathrm{logit}(\mathbf{p})\mid\boldsymbol{\beta},\sigma_{u}^{2},\phi\sim\boldsymbol{\mathcal{N}}(\mathbf{X}\boldsymbol{\beta},\boldsymbol{\Sigma}_{BYM2}(\sigma_{u}^{2},\phi)).$
(9)

In the above, $p_{a}$ denotes the finite population area-specific proportion
and $V_{a}$ denotes the sampling variance of the direct estimator
$\widehat{p}_{a}$. We use the shorthand $\mathrm{logit}(\mathbf{p})$ to denote
the vector $(\mathrm{logit}(p_{1}),\ldots,\mathrm{logit}(p_{A}))^{T}$, which
we assume is drawn from a multivariate Gaussian distribution with mean
$\mathbf{X}\boldsymbol{\beta}$, where $\mathbf{X}$ is an $A\times(p+1)$ design
matrix containing covariate information and $\boldsymbol{\beta}$ is a
$(p+1)$-vector containing the intercept and corresponding coefficients.
Finally $\boldsymbol{\Sigma}_{BYM2}(\sigma_{u}^{2},\phi)$ denotes a spatial
covariance matrix dependent on marginal variance parameter $\sigma_{u}^{2}$
and spatial correlation parameter $\phi$. We use the BYM2 model, a
reparametrization of the Besag-York-Mollié [30] model proposed by Riebler et
al. (2016) [31] which determines the structure of
$\boldsymbol{\Sigma}_{BYM2}(\sigma_{u}^{2},\phi)$. Below, we review the BYM2
model, rewriting the mean linking model as follows for clarity:

$\mathrm{logit}(\mathbf{p})=\mathbf{X}\boldsymbol{\beta}+\mathbf{u}$ (10)
$\mathbf{u}\mid\boldsymbol{\sigma}_{u}^{2},\phi\sim\boldsymbol{\mathcal{N}}(\mathbf{0},\boldsymbol{\Sigma}_{BYM2}(\sigma_{u}^{2},\phi)).$
(11)

Under the BYM2 model, we assume $\mathbf{u}$ can be partitioned into an
unstructured component $\widetilde{\mathbf{u}}_{1}$ and a structured spatial
component $\widetilde{\mathbf{u}}_{2*}$:

$\mathbf{u}=\sigma_{u}\left(\sqrt{1-\phi}\widetilde{\mathbf{u}}_{1}+\sqrt{\phi}\widetilde{\mathbf{u}}_{2*}\right)$
(12)

We assume $\widetilde{\mathbf{u}}_{1}\sim N(\mathbf{0},\mathbf{I})$ is a
vector of iid Gaussian random area effects and assume an intrinsic conditional
autoregressive (ICAR) Gaussian prior for $\widetilde{\mathbf{u}}_{2*}$. The
ICAR prior, as proposed by Besag et al. (1991) [30], assumes that spatial
components $\widetilde{u}_{2a*}$ and $\widetilde{u}_{2b*}$, representing the
values of $\widetilde{\mathbf{u}}_{2*}$ for areas $a$ and $b$, are correlated
if areas $a$ and $b$ are defined to be neighbors. Under an ICAR prior, we
assume that for a particular area $a$, the mean of $\widetilde{u}_{2a*}$ is
equal to the mean of all neighboring effects and the precision of
$\widetilde{u}_{2a*}$ is proportional to the number of neighbors. Using this
parameterization, $\sigma_{u}$ denotes the marginal variance of $\mathbf{u}$
and $\phi$ represents the proportion of variation assigned to the spatial
component.

Under a BYM2 model, $\mathbf{u}$ has the covariance matrix

$\boldsymbol{\Sigma}_{BYM2}(\sigma_{u}^{2},\phi)=\sigma_{u}((1-\phi)\mathbf{I}+\phi\mathbf{Q}_{*}^{-})$
(13)

Here, $\mathbf{Q}_{*}$ denotes the precision matrix of
$\widetilde{\mathbf{u}}_{2*}$ and $\mathbf{Q}_{*}^{-}$ is its generalized
inverse. Note that the precision matrix implied by the ICAR prior,
$\mathbf{Q}_{*}$, is singular, yielding an improper prior. To ensure
identifiability, we must place a sum-to-zero constraint on $\mathbf{u}$. In
order to make the marginal variance parameter $\sigma_{u}$ interpretable, we
scale $\mathbf{Q}_{*}$ to make the geometric mean of the marginal variances
equal to one, as recommended by Riebler et al. (2016) [31]

###  4.2 Variance model

We similarly use unmatched models for the corresponding variance estimates
$\widehat{\mathbf{V}}$, using a chi-squared sampling model with a log-normal
linking model. We use the chi-squared sampling model described in Equation
(3), assuming that for all $a$, the variance estimate $\widehat{V}_{a}$ is an
unbiased estimator of $V_{a}$. The linking model assumes the true
$\log(V_{a})$ values are Gaussian distributed with expected values given by a
generalized variance function $f(p_{a},\mathbf{z}_{a};\boldsymbol{\gamma})$
whose inputs are the area proportion $p_{a}$, other area level predictors
$\mathbf{z}_{a}$, and parameters $\boldsymbol{\gamma}$. We can write down the
unmatched models as follows:

$\frac{d_{a}\widehat{V}_{a}}{V_{a}}\mid d_{a},V_{a}\stackrel{{\scriptstyle
ind}}{{\sim}}\chi^{2}_{d_{a}}$ (14) $\log(V_{a})\mid
p_{a},\mathbf{z}_{a},\boldsymbol{\gamma},\sigma_{\tau}^{2}\stackrel{{\scriptstyle
ind}}{{\sim}}N\left(f(p_{a},\mathbf{z}_{a};\boldsymbol{\gamma}),\sigma_{\tau}^{2}\right)$
(15)

Here, $d_{a}$ denotes the degrees of freedom parameter for area $a$, which we
determine based on the survey design and sample size as discussed below. We
use $\sigma_{\tau}^{2}$ to denote the variance of the linking model errors
which allow for area-specific deviations from the linking model.

We define the generalized variance function as follows:

$f(p_{a},\mathbf{z}_{a};\boldsymbol{\gamma})=\gamma_{0}+\gamma_{1}\log(p_{a}(1-p_{a}))+\gamma_{2}\log(n_{a})$
(16)

where $n_{a}$ denotes the sample size for area $a$. Note that if we set
$\gamma_{0}=0,\gamma_{1}=1,\gamma_{2}=-1$, the right hand side resembles the
logarithm of the binomial variance. As such, this GVF can be viewed as a
generalized version of the binomial variance. The GVF used here could also be
altered to introduce additional covariates or different functional
relationships between $p_{a}$ and $V_{a}$. We can view the variance linking
model (15) as a prior that shrinks the estimate $\widehat{V}_{a}$ towards a
model-based prediction dependent on the binomial mean-variance relationship.

As described above, the mean linking model induces spatial smoothing for
estimates $\widehat{p}_{a}$. By combining the mean and variance models and
incorporating the means $p_{a}$ into the GVF, we induce spatial correlation
into the resulting samples of $V_{a}$, potentially aiding estimation in areas
with fewer samples.

We treat the degrees of freedom parameter $d_{a}$ as known for all areas $a$.
The appropriate choice for $d_{a}$ depends on the sampling design. As
mentioned above, if the data for a given area were iid Gaussian (for example,
reflecting simple random sampling with replacement), the typical variance
estimator would follow a $\chi^{2}$ distribution with $d_{a}=n_{a}-1$ degrees
of freedom. However, for sampling without replacement and cluster sampling
designs, other choices of $d_{a}$ may be more appropriate depending on how
$\widehat{V}_{a}$ is computed for each area. Maples et al. (2009) [24] outline
a resampling procedure for estimating degrees of freedom for their variance
sampling model. When computing variance estimates from DHS data, we use a
simplified variance estimator based on the with-replacement variance estimator
for multistage designs presented in Equation (4.6.2) of Särndal et al.
(1993)[32], which is computed as a sum over clusters:

$\widehat{V}_{a}=\frac{1}{m_{a}(m_{a}-1)}\sum_{j\in
S_{1a}}\left(\frac{\widehat{t}_{aj}}{\pi_{j}}-\widehat{t}_{a}\right)^{2}$ (17)

where $S_{1a}$ denotes the set of indices of sampled clusters for area $a$,
$m_{a}$ denotes the number of sampled clusters, and $\pi_{j}$ denotes the
probability of sampling cluster $j$. Finally $\widehat{t}_{aj}$ denotes the
direct estimator for the total for cluster $j$ in area $a$ and
$\widehat{t}_{a}$ denotes the direct estimator for the total of area $a$.
Since this is a sum of squared error terms over $m_{a}$ clusters, we set
$d_{a}$ to be equal to $m_{a}-1$.

###  4.3 Estimation

We adopt a fully Bayesian approach to estimation by placing priors on the
following hyperparameters:

$\\{\boldsymbol{\beta},\sigma_{u}^{2},\phi,\boldsymbol{\gamma},\sigma_{\tau}^{2}\\}\sim\Pi(\boldsymbol{\theta})$
(18)

where $\boldsymbol{\theta}$ denotes any parameters used to specify the priors.
Details on the priors used in each example are provided in Appendix A. We
compute approximate posterior distributions for $p_{a}$ for all areas $a$
using Markov chain Monte Carlo sampling as implemented in the Stan programming
language [33]. Functions for fitting the models described above have been
collected in an R package called VSALM available at
https://github.com/peteragao/VSALM and the code for the below simulations and
analysis is available at the associated repository
https://github.com/peteragao/VSALM-paper. In this context, the Bayesian
approach offers a number of potential benefits. In particular, we are able to
sample from the joint posterior distributions for the proportions of interest
for all areas, giving a natural way to quantify uncertainty and also enabling
comparisons between areas. Moreover, the sampling approach implemented in Stan
is fast and flexible, enabling users to fit and compare potential models
quickly.

##  5 Simulations

###  5.1 Population generating model

We use simulations to evaluate our spatial variance smoothing estimator,
comparing its performance with that of the direct weighted Hájek estimator and
an estimator derived from a model without variance smoothing. For our
simulations, we generate an artificial population that mimics data from the
2018 Nigeria DHS. First, we generate synthetic cluster locations across
Nigeria using a pixel grid of estimated population counts for Nigeria in 2006
(mimicking the sampling frame used for the DHS survey) [34]. For each of the
73 strata used for the DHS, we sample 300 pixels without replacement with
probability proportional to population. These sampled pixels represent
enumeration areas or clusters. For each cluster location, we randomly generate
cluster sizes $N_{c}\sim\mathrm{Poisson(10)}$, yielding a population of
$N=\sum\limits_{c}N_{c}$ individuals.

For each individual $i$ in our population, we generate data using a population
generating model motivated by the models used by Corral et al. (2020) [35]
(see Section 7.2) and Gao and Wakefield (2022) [36]. For each cluster $c$, we
simulate cluster level covariate information as follows:

  1. 1.

The covariate $x_{1,c}$ is the realized value of a binary random variable
$X_{1,c}$ with $P(X_{1,c}=1)=0.5$;

  2. 2.

The covariate $x_{2,c}$ is the realized value of a binary random variable
$X_{2,c}$ with $P(X_{2,c}=1)=0.3+0.5\frac{a(c)}{37}$, where $a(c)$ is the
index of the area containing cluster $c$;

  3. 3.

The covariate $x_{3,c}=x_{3,a(c)}$ is obtained from a $37\times 1$ ICAR random
vector with marginal variance 1 for the Admin-1 areas.

  4. 4.

The covariate $x_{4,c}=x_{4,a(c)}$ is obtained from a $774\times 1$ ICAR
random vector with marginal variance 1 for the Admin-2 areas.

  5. 5.

The covariate $x_{5,c}$ is obtained from a random vector generated using a
stochastic partial differential equation (SPDE) -based approximation [37] to a
Gaussian process with Matérn covariance with smoothness parameter 1 and
marginal variance of 1.

Maps of these covariates are provided in Appendix B.1. Based on these
covariates, we simulate a cluster level risk parameter $q_{c}$ for each
cluster from the following model:

$\mathrm{logit}(q_{c})=\mathrm{logit}(\mu)+0.25x_{1,c}-0.25x_{2,c}+0.5x_{3,c}+0.25x_{4,c}+0.25x_{5,c}+u_{a(i)}+v_{c}$
(19)

where $u_{a(i)}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.25^{2})$ are
independent and identically distributed area level random effects, and
$v_{c}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.5^{2})$ represents random and
independent and identically distributed cluster level effects. In the above,
$\mu$ denotes the global superpopulation prevalence. The covariates are held
constant for all simulations, but the response variables and random area and
cluster effects are resampled for each new simulation. We repeatedly generate
$Y_{i}\mid q_{c(i)}\sim\mathrm{Bernoulli}(q_{c(i)})$, where $c(i)$ is the
cluster of individual $i$. For each simulation, we can thus compute true
population Admin-1 area level proportions $p_{a}$.

To obtain our simulated samples, we use a cluster sampling design. In each
simulation, we sample eight clusters from each stratum, keeping all
individuals in sampled clusters. We compute sampling weights $w_{i}$ for each
individual from the corresponding inverse inclusion probabilities.

We compare our unmatched joint smoothing model-based estimator with the direct
Hájek estimator and a number of alternative model-based estimates. First, we
consider spatial joint sampling (Spatial Unmatched JS) and non-spatial joint
sampling (Unmatched JS) models, where the non-spatial version is obtained by
replacing the BYM2 prior for the area effects $\mathbf{u}$ with an iid
multivariate Gaussian prior. We also consider an estimator produced using a
model that omits the variance smoothing model entirely, which we refer to as
the mean smoothing (MS) model-based estimator. This model is specified using
the unmatched models (8) and (9) but treating $\widehat{V}_{a}$ as known for
all $a$. We consider both spatial (Spatial Unmatched MS) and non-spatial
versions (Unmatched MS). For all model-based estimators, we adopt a fully
Bayesian approach to inference as described in Section 4. We obtain point
estimates $\widehat{p}_{a}$ and 90% credible interval estimates
$(l_{a},u_{a})$ by sampling from these approximate posterior distributions.
Further details on the estimation procedure are provided in the Appendix.

For each simulation, we calculate several performance metrics, including root
mean squared error (RMSE) and mean absolute error (MAE). We also compute the
average coverage of the 90% interval estimates and the mean interval lengths
(MIL) across all areas. For a single simulation, these metrics are defined as
follows.

$\displaystyle\mathrm{RMSE}(\widehat{p}_{a})$
$\displaystyle=\sqrt{\frac{1}{A}\sum_{a}(p_{a}-\widehat{p}_{a})^{2}}$ (20)
$\displaystyle\mathrm{MAE}(\widehat{p}_{a})$
$\displaystyle=\frac{1}{A}\sum_{a}|p_{a}-\widehat{p}_{a}|$ (21)
$\displaystyle\mathrm{Cov}_{90}(\widehat{p}_{a})$
$\displaystyle=\frac{1}{A}\sum_{a}\mathbf{1}\\{p_{a}\in(l_{a},u_{a})\\}$ (22)
$\displaystyle\mathrm{MIL}_{90}(\widehat{p}_{a})$
$\displaystyle=\frac{1}{A}\sum_{a}(u_{a}-l_{a})$ (23)

| $\mu=0.1$ | $\mu=0.5$
---|---|---
Method | RMSE | MAE | 90% Cov. | MIL | RMSE | MAE | 90% Cov. | MIL
Direct (Hájek) | 3.44 | 2.63 | 83 | 10.42 | 6.16 | 4.92 | 85 | 19.55
Unmatched MS | 3.53 | 2.55 | 82 | 9.10 | 5.80 | 4.62 | 85 | 17.48
Spatial Unmatched MS | 3.28 | 2.39 | 83 | 8.60 | 5.61 | 4.46 | 85 | 17.20
Unmatched JS | 3.24 | 2.46 | 89 | 9.78 | 6.44 | 5.06 | 89 | 20.20
Spatial Unmatched JS | 3.03 | 2.28 | 90 | 9.34 | 6.19 | 4.87 | 90 | 19.72
Table 1: RMSE ($\times 100$), MAE ($\times 100$), coverage rates, and mean
interval length ($\times 100$) of estimators of area level means across 1,000
simulated populations with spatially correlated binary responses based on
sample data obtained via informative sampling. The reduced model omits one of
the spatial covariates in the full model.

We consider two sets of simulations with differing global prevalence rates. We
let $\mu=0.1$ for the first set, which is similar to the overall HIV
positivity rate in the Malawi data. For the second set, we let $\mu=0.5$,
which is similar to the national MCV-1 vaccination rate in the Nigeria data.
Table 1 summarizes results for our two sets of simulations. In each setting,
the results represent the average values of the metrics (20)-(23) across 1,000
simulated populations. We observe that in the low prevalence examples, the
spatial unmatched joint model-based estimates perform best in terms of RMSE
and MAE. Moreover, prediction intervals constructed based on the direct
estimator and the mean-only smoothing model-based estimators tend to exhibit
undercoverage, whereas the joint model-based intervals achieve closer to
nominal coverage. In the moderate prevalence examples, the joint model-based
estimates perform slightly worse than the mean-only model-based estimates in
terms of the RMSE and MAE; however, the Hájek and mean-only model intervals
show slight undercoverage. The joint modeling approach thus yields slightly
more conservative prediction intervals which may be desirable for decision
making.

##  6 Applications

We apply our joint smoothing model-based estimator to two examples involving
DHS data, demonstrating its use for a low prevalence indicator (Malawi HIV
prevalence rates) and for a moderate prevalence indicator (Nigeria measles
vaccination rates). We show that our method induces spatial smoothing of
estimated variances and produces more conservative interval estimates than an
approach using an area level model that only smooths means.

For both examples, we compare direct weighted estimation with the model-based
smoothing methods described above. We first fit both the spatial mean
smoothing unmatched model (Spatial Unmatched MS), which omits the variance
model, as well as the full spatial joint smoothing unmatched model (Spatial
Unmatched JS) For all models, we use no covariates and we compute approximate
posterior distributions for all area level proportions $p_{a}$ and obtain
corresponding point and interval estimates by sampling from these posteriors.

Figure 3: Direct and model-based point estimates (top) and length of
corresponding 90% interval estimates (bottom) of vaccination coverage rate for
first dose of measles-containing-vaccine (MCV1) among children aged 12–23
months in Nigeria, 2018. Figure 4: Direct and model-based point estimates
(top) and length of corresponding 90% interval estimates (bottom) of HIV
prevalence rate for women aged 15-49 in Malawi, 2015-2016.

Figure 3 compares point estimates of MCV-1 coverage rates (top) and the length
of interval estimates (bottom) for Admin-1 areas among children aged 12-23
months in Nigeria in 2018. Figure 4 similarly provides point estimates of HIV
prevalence rates (top) and the length of interval estimates (bottom) for
Admin-1 areas for women aged 15-49 in Malawi, 2015-2016. For both examples,
the bottom set of maps illustrates the estimated uncertainty of the direct and
model-base estimates using the length of 90% credible intervals. In general,
we observe that the point estimates agree well for all three methods. However,
we observe some spatial smoothing of the interval lengths, suggesting that the
joint smoothing model induces spatial smoothing of the direct estimator
variances $V_{a}$.

Figure 5: Comparison of model-based 90% credible interval lengths with Hájek
90% confidence interval lengths for Malawi HIV example (left) and Nigeria MCV
example (right).

Figure 5 provides a scatter plot comparing the model-based 90% credible
interval lengths produced by the mean smoothing and joint smoothing models
with the design-based 90% confidence interval length associated with the Hájek
estimator. Here, we see that for the Malawi example (left), the credible
interval lengths are similar for the mean smoothing and joint smoothing
models. However, for the Nigeria example (right), the joint smoothing
intervals are more conservative when the Hájek intervals are short and
narrower when the Hájek intervals are long. In addition to the spatial
smoothing of interval lengths seen in Figures 3 and 4, this suggests the joint
model can help smooth variance estimates globally.

##  7 Discussion

Our proposed model-based estimator for small area proportions combines several
proposed extensions of the basic Fay-Herriot area level model, including using
unmatched linking and sampling models to address the non-Gaussian response
data [38], incorporating spatial smoothing via correlated random effects in
the mean linking model [28], and introducing a variance smoothing model so
that the resulting estimators exhibit both smoothed means and variances [12,
13, 16]. We propose a spatial joint smoothing model and adopt a fully Bayesian
approach to estimation, which facilitates quick computation of point and
interval estimates. Through simulation and application, we have shown that
inferences based on our model can improve upon those based on a model that
only incorporates smoothing of means. Interval estimates obtained from our
model can correct for the undercoverage seen in models that only smooth means,
suggesting our model may more accurately account for uncertainty in estimated
variances of direct weighted estimators.

For our clustered binary response data, the variance smoothing model we have
adopted may help address undercoverage of interval estimates caused by
treating variances of direct weighted estimators as known. However, for other
designs and contexts, such a model may be inappropriate. In general, the
choice of variance sampling and linking models should depend on a number of
factors including any clustering and stratification in the design as well as
the distribution and presumed mean-variance relationship of the response
variables.

Moreover, we acknowledge that our variance smoothing model is a simplification
of the true distribution of design-based variances $V_{a}$. in particular, the
use of a chi-squared distribution for the variance sampling model relies upon
the assumption that the direct estimator of variance $\widehat{V}_{a}$ for a
particular area $a$ is computed as the sum of several squared Gaussian terms.
Since our data is non-Gaussian, this assumption may be violated and other
sampling models for $\widehat{V}_{a}$, such as a Gaussian model, could be
explored i future work. Moreover, within each area, we have assumed that the
variance for each stratum is equal, but this assumption may be inappropriate
for our data since each area of interest is divided into urban and rural
subregions, which may be qualitatively different from one another. Finally,
the appropriate number of degrees of freedom $d_{a}$ depends on the specific
design used; using resampling methods like those explored by Maples et al.
(2009) [24] to choose $d_{a}$ may help improve the fit of the variance
sampling model.

Although we have presented one approach applied to two different types of
problems, in practice, for decision making purposes, different estimation
problems may have varying priorities. For example, when implementing targeted
vaccination programs, it is important to identify communities with especially
low vaccination rates, whereas designing policy for providing resources
associated with HIV involves identifying communities with high rates of
positivity. Given that the variance of a direct prevalence rate estimator may
depend on its expected value, various modeling decisions such as choosing to
apply a transformation for $\widehat{p}_{a}$ may lead to different results
depending on the expected value of $\widehat{p}_{a}$. As such, it is crucial
to carefully consider the distribution of direct estimators before selecting a
model. Above, we have used unmatched sampling and linking models for the area
level proportions, but we also considered first computing the logit-
transformed direct estimators and then applying matched sampling and linking
models treating both $\widehat{p}_{a}$ and $p_{a}$ as Gaussian random
variables. In our simulations and application, this approach did not
outperform the unmatched models we adopted, but future research could help
illustrate when such an approach could be useful.

When mapping subnational health and demographic indicators in LMICs, unit
level models, and in particular geostatistical models using spatial Gaussian
processes, are often used as they allow estimates to be generated at arbitrary
resolutions and can incorporate unit level covariate information. However,
such approaches may often struggle to account for design effects such as those
caused by clustering and informative sampling. While unit level models may be
able to generate prevalence estimates at the individual cluster level,
aggregating those cluster level estimates upwards to produce area level
estimates may introduce additional errors and lead to improperly calibrated
interval estimates [2, 39, 40]. Area level models are specified to generate
estimates for a preselected set of regions. Moreover, area level models are
often simpler and faster to implement than unit level models. For these
reasons, we have explored the feasibility of using area level models to
generate maps of health indicators such as vaccination rates and disease
prevalence rates in LMICs. Our method, like many area level methods, directly
accounts for survey design by incorporating available sampling weight
information. By incorporating a spatial variance smoothing model and using
unmatched sampling and linking models, we are able to address some of the
difficulties related to applying area level models for use in this specific
context.

## References

  * [1] Parker PA, Holan SH, Janicki R. Computationally efficient Bayesian unit-level models for non-Gaussian data under informative sampling with application to estimation of health insurance coverage. The Annals of Applied Statistics 2022; 16(2): 887–904.  doi: 10.1214/21-AOAS1524
  * [2] Fuglstad GA, Li ZR, Wakefield J. The Two Cultures for Prevalence Mapping: Small Area Estimation and Spatial Statistics. arXiv:2110.09576 [stat] 2022\.
  * [3] Fay RE, Herriot RA. Estimates of Income for Small Places: An Application of James-Stein Procedures to Census Data. Journal of the American Statistical Association 1979; 74(366a): 269–277.  doi: 10.1080/01621459.1979.10482505
  * [4] Rao JNK, Molina I. Small Area Estimation.  Wiley Series in Survey MethodologyJohn Wiley & Sons.  2nd ed. 2015.
  * [5] Arora V, Lahiri P. On the Superiority of the Bayesian Method Over the BLUP in Small Area Estimation Problems. Statistica Sinica 1997; 7(4): 1053–1063.
  * [6] Rivest LP, Vandal N. Mean squared error estimation for small areas when the small area variances are estimated. In: ; 2002; Ottawa, Canada.
  * [7] Bell WR. Examining Sensitivity of Small Area Inferences to Uncertainty About Sampling Error Variances. working Paper, United States Census Bureau; 2008.
  * [8] National Population Commission - NPC/Nigeria and ICF. . Nigeria Demographic and Health Survey 2018. 2019.
  * [9] ICF NSONa. Malawi Demographic and Health Survey 2015-16. 2017.
  * [10] Kleffe J, Rao J. Estimation of mean square error of empirical best linear unbiased predictors under a random error variance linear model. Journal of Multivariate Analysis 1992; 43(1): 1–15.  doi: 10.1016/0047-259X(92)90107-Q
  * [11] Wang J, Fuller WA. The Mean Squared Error of Small Area Predictors Constructed with Estimated Area Variances. Journal of the American Statistical Association 2003; 98(463): 716–723.
  * [12] You Y, Chapman B. Small Area Estimation Using Area Level Models and Estimated Sampling Variances. Survey Methodology 2006; 32(1): 97–103.
  * [13] Maiti T, Ren H, Sinha S. Prediction Error of Small Area Predictors Shrinking Both Means and Variances. Scandinavian Journal of Statistics 2014; 41(3): 775–790.
  * [14] Hwang JTG, Qiu J, Zhao Z. Empirical Bayes Confidence Intervals Shrinking Both Means and Variances. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 2009; 71(1): 265–285.
  * [15] Dass SC, Maiti T, Ren H, Sinha S. Confidence interval estimation of small area parameters shrinking both means and variances. Survey Methodology 2012; 38(2): 173–187.
  * [16] Sugasawa S, Tamae H, Kubokawa T. Bayesian Estimators for Small Area Models Shrinking Both Means and Variances. Scandinavian Journal of Statistics 2017; 44(1): 150–167.  doi: 10.1111/sjos.12246
  * [17] Polettini S. A Generalised Semiparametric Bayesian Fay–Herriot Model for Small Area Estimation Shrinking Both Means and Variances. Bayesian Analysis 2017; 12(3): 729–752.  doi: 10.1214/16-BA1019
  * [18] Wolter K. Introduction to Variance Estimation.  Springer Science & Business Media .  2007\.
  * [19] Liu B, Lahiri P, Kalton G. Hierarchical Bayes Modeling of Survey-Weighted Small Area Proportions. Survey Methodology 2014; 40(1): 1–13.
  * [20] Hawala S, Lahiri P. Variance Modeling for Domains. Statistics and Applications 2018; 16(1): 399–409.
  * [21] Maples JJ. Estimating Design Effects in Small Areas/Domains through Aggregation. In: ; 2016: 671–681.
  * [22] Franco C, Bell WR. Applying Bivariate Binomial/Logit Normal Models to Small Area Estimation. 2013: 690–702.
  * [23] Mohadjer L, Rao JNK, Liu B, Krenzke T, Kerckhove dWV. Hierarchical Bayes Small Area Estimates of Adult Literacy using Unmatched Sampling and Linking Models. Journal of the Indian Society of Agricultural Statistics 2012: 55–63.
  * [24] Maples JJ, Bell WR, Huang ET. Small-Area Variance Modeling with Application to County Poverty Estimates from the American Community Survey. In: ; 2009: 5056–5067.
  * [25] Mercer LD, Wakefield J, Pantazis A, Lutambi AM, Masanja H, Clark S. Space-Time Smoothing of Complex Survey Data: Small Area Estimation for Child Mortality. The annals of applied statistics 2015; 9(4): 1889–1905.  doi: 10.1214/15-AOAS872
  * [26] Chen C, Wakefield J, Lumley T. The use of sampling weights in Bayesian hierarchical models for small area estimation. Spatial and Spatio-temporal Epidemiology 2014; 11: 33–43.  doi: 10.1016/j.sste.2014.07.002
  * [27] Ghosh M, Natarajan K, Stroud TWF, Carlin BP. Generalized Linear Models for Small-Area Estimation. Journal of the American Statistical Association 1998; 93(441): 273–282.  doi: 10.2307/2669623
  * [28] Pratesi M, Salvati N. Small area estimation: the EBLUP estimator based on spatially correlated random area effects. Statistical Methods and Applications 2008; 17(1): 113–141.  doi: 10.1007/s10260-007-0061-9
  * [29] Chung HC, Datta GS. Bayesian Hierarchical Spatial Models for Small Area Estimation. Tech. Rep. 2020-07, Center for Statistical Research & Methodology, Research and Methodology Directorate, U.S. Census Bureau; 2020.
  * [30] Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 1991; 43(1): 1–20.  doi: 10.1007/BF00116466
  * [31] Riebler A, Sørbye SH, Simpson D, Rue H. An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research 2016; 25(4): 1145–1165.  doi: 10.1177/0962280216660421
  * [32] Särndal CE, Swensson B, Wretman J. Model Assisted Survey Sampling.  Springer Science & Business Media .  2003\.
  * [33] Carpenter B, Gelman A, Hoffman MD, et al. Stan: A Probabilistic Programming Language. Journal of Statistical Software 2017; 76: 1–32.  doi: 10.18637/jss.v076.i01
  * [34] WorldPop , Center for International Earth Science Information Network (CIESIN), Columbia University . Global High Resolution Population Denominators Project - Funded by The Bill and Melinda Gates Foundation (OPP1134076). 2006.
  * [35] Corral P, Molina I, Nguyen M. Pull Your Small Area Estimates Up by the Bootstraps.  Policy Research Working PapersThe World Bank .  2020\.
  * [36] Gao PA, Wakefield J. Smoothed Model-Assisted Small Area Estimation. 2022\.  Number: arXiv:2201.08775 arXiv:2201.08775 [stat]
  * [37] Lindgren F, Rue H, Lindström J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2011; 73(4): 423–498.  doi: 10.1111/j.1467-9868.2011.00777.x
  * [38] You Y, Rao JNK. Small Area Estimation Using Unmatched Sampling and Linking Models. The Canadian Journal of Statistics / La Revue Canadienne de Statistique 2002; 30(1): 3–15.  doi: 10.2307/3315862
  * [39] Paige J, Fuglstad GA, Riebler A, Wakefield J. Design- and Model-Based Approaches to Small-Area Estimation in a Low- and Middle-Income Country Context: Comparisons and Recommendations. Journal of Survey Statistics and Methodology 2022; 10(1): 50–80.  doi: 10.1093/jssam/smaa011
  * [40] Paige J, Fuglstad GA, Riebler A, Wakefield J. Spatial Aggregation with Respect to a Population Distribution. 2022.  Number: arXiv:2207.06700 arXiv:2207.06700 [stat]
  * [41] Simpson D, Rue H, Riebler A, Martins TG, Sørbye SH. Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors. Statistical Science 2017; 32(1): 1–28.  doi: 10.1214/16-STS576

##  Appendix A Parameter estimation

Below we provide further details on the estimation process and priors used for
each method described in the simulations and applications.

###  A.1 Direct estimation

We use the R package survey for computing direct weighted Hájek estimators
(and corresponding variance estimates) for all areas.

###  A.2 Mean-smoothing model-based estimation

We adopt a fully Bayesian approach to estimating the mean-smoothing unmatched
model described above, assuming priors for model parameters and then using
MCMC as implemented in the R package STAN to sample from the posterior
distributions of the area level proportions $p_{a}$ for all $a=1,\ldots,A$. We
place a $N(0,1000)$ prior on the area level model intercept and fixed effects.
We use penalized complexity priors for the variance parameter $\sigma_{u}$, as
described by Simpson et al. (2017) [41] We specify these priors such that
$P(\sigma_{u}>1)=0.01$. For the spatial models, we place a
$\mathrm{Beta}(1/2,1/2)$ prior on the spatial correlation prior $\phi$.

###  A.3 Joint-smoothing model-based estimation

We keep the same priors as for the mean-smoothing model on area level model
intercept and fixed effects as well as the variance parameters $\sigma_{u}$
and for the spatial models, $\phi$. We then place a penalized complexity prior
for the variance parameter $\sigma_{\tau}$ such that
$P(\sigma_{\tau}>1)=0.01$. We place $N(0,1)$ prior on $\gamma_{0}$ and
$N(1,.5)$ prior on $\gamma_{1}$ and $N(-1,.5)$ prior on $\gamma_{2}$ to shrink
the resulting variances estimates towards that of a binomial random variable.

##  Appendix B Additional results

###  B.1 Covariate maps

Figure 6 provides maps of the simulated cluster locations and covariate values
used in the simulations described in Section 5.

Figure 6: Simulated cluster locations and covariate values used to generate
population data.

###  B.2 Large sample simulations

Table 2 provides results for an additional set of simulations that were
identical to the simulations with $\mu=0.5$ described in Section 5, except
with a larger sample size, with twenty-five clusters sampled per stratum
rather than the eight used in the main text. The results illustrate that for
large sample sizes, the joint smoothing and mean smoothing model-based
estimators perform similarly to the direct weighted estimators with 90%
prediction interval coverage rates that are close to nominal.

Method | RMSE | MAE | 90% Cov. | MIL
---|---|---|---|---
Direct (Hájek) | 2.40 | 1.91 | 90 | 7.98
Unmatched MS | 2.35 | 1.87 | 90 | 7.87
Spatial Unmatched MS | 2.33 | 1.86 | 90 | 7.80
Unmatched JS | 2.37 | 1.89 | 91 | 7.99
Spatial Unmatched JS | 2.35 | 1.87 | 91 | 7.92
Table 2: RMSE ($\times 100$), MAE ($\times 100$), coverage rates, and mean
interval length ($\times 100$) of estimators of area level means across 1,000
simulated populations with spatially correlated binary responses based on
sample data obtained via informative sampling. The reduced model omits one of
the spatial covariates in the full model.

###  B.3 Applications

Tables 3 and 4 provide point estimates and corresponding 90% prediction
intervals for model hyperparameters.

Tables 5 and 6 provide full estimates with prediction intervals for all areas
of interest for the methods described in the manuscript applied to both the
measles vaccination rate and HIV prevalence rate applications.

Parameter | Spatial Unmatched MS | Spatial Unmatched JS
---|---|---
$\mathrm{logit}(\mu)$ | 0.36 | (0.22, 0.5) | 0.36 | (0.23, 0.49)
$\sigma_{u}$ | 0.71 | (0.54, 0.96) | 0.71 | (0.54, 0.94)
$\phi$ | 0.84 | (0.38, 1) | 0.86 | (0.39, 1)
$\gamma_{0}$ |  |  | -0.73 | (-2.22, 0.74)
$\gamma_{1}$ |  |  | 0.47 | (-0.14, 1.01)
$\gamma_{2}$ |  |  | -0.9 | (-1.18, -0.63)
$\sigma_{\tau}$ |  |  | 0.42 | (0.29, 0.59)
Table 3: Point estimates and 90% interval estimates for model parameters for Nigeria measles vaccination example Parameter | Spatial Unmatched MS | Spatial Unmatched JS
---|---|---
$\mathrm{logit}(\mu)$ | -2.03 | (-2.14, -1.93) | -2.02 | (-2.14, -1.92)
$\sigma_{u}$ | 0.41 | (0.28, 0.62) | 0.43 | (0.3, 0.64)
$\phi$ | 0.89 | (0.36, 1) | 0.88 | (0.36, 1)
$\gamma_{0}$ |  |  | -0.22 | (-1.9, 1.5)
$\gamma_{1}$ |  |  | 0.91 | (0.51, 1.27)
$\gamma_{2}$ |  |  | -0.96 | (-1.29, -0.64)
$\sigma_{\tau}$ |  |  | 0.19 | (0.07, 0.37)
Table 4: Point estimates and 90% interval estimates for model parameters for Malawi HIV prevalence example State | Hájek | Spatial Unmatched MS | Spatial Unmatched JS
---|---|---|---
Lagos | 0.89 | (0.84, 0.94) | 0.88 | (0.83, 0.93) | 0.89 | (0.83, 0.93)
Ekiti | 0.87 | (0.79, 0.94) | 0.83 | (0.76, 0.89) | 0.83 | (0.75, 0.89)
Anambra | 0.80 | (0.73, 0.88) | 0.79 | (0.73, 0.85) | 0.79 | (0.72, 0.85)
Enugu | 0.80 | (0.71, 0.88) | 0.77 | (0.7, 0.85) | 0.77 | (0.69, 0.85)
Edo | 0.79 | (0.7, 0.88) | 0.78 | (0.7, 0.85) | 0.78 | (0.69, 0.85)
Osun | 0.77 | (0.69, 0.85) | 0.76 | (0.69, 0.83) | 0.76 | (0.68, 0.83)
Abia | 0.75 | (0.69, 0.82) | 0.75 | (0.69, 0.81) | 0.75 | (0.68, 0.82)
Delta | 0.75 | (0.69, 0.8) | 0.75 | (0.7, 0.8) | 0.76 | (0.69, 0.82)
Abuja | 0.73 | (0.68, 0.79) | 0.72 | (0.67, 0.77) | 0.72 | (0.65, 0.78)
Imo | 0.73 | (0.63, 0.84) | 0.74 | (0.65, 0.83) | 0.73 | (0.64, 0.82)
Bayelsa | 0.73 | (0.65, 0.8) | 0.73 | (0.66, 0.81) | 0.73 | (0.65, 0.81)
Ondo | 0.69 | (0.58, 0.8) | 0.70 | (0.6, 0.79) | 0.69 | (0.6, 0.78)
Rivers | 0.68 | (0.59, 0.77) | 0.69 | (0.61, 0.77) | 0.69 | (0.61, 0.77)
Cross River | 0.65 | (0.52, 0.78) | 0.66 | (0.55, 0.77) | 0.66 | (0.55, 0.76)
Adamawa | 0.65 | (0.57, 0.73) | 0.63 | (0.56, 0.71) | 0.63 | (0.55, 0.71)
Nassarawa | 0.65 | (0.54, 0.75) | 0.63 | (0.54, 0.72) | 0.63 | (0.54, 0.72)
Ebonyi | 0.63 | (0.57, 0.7) | 0.64 | (0.58, 0.7) | 0.64 | (0.58, 0.71)
Akwa Ibom | 0.63 | (0.55, 0.71) | 0.64 | (0.56, 0.72) | 0.64 | (0.56, 0.73)
Benue | 0.63 | (0.54, 0.71) | 0.63 | (0.55, 0.71) | 0.63 | (0.55, 0.7)
Oyo | 0.60 | (0.51, 0.7) | 0.61 | (0.52, 0.7) | 0.61 | (0.52, 0.7)
Plateau | 0.59 | (0.52, 0.65) | 0.58 | (0.52, 0.65) | 0.58 | (0.51, 0.65)
Kano | 0.56 | (0.5, 0.62) | 0.56 | (0.5, 0.62) | 0.56 | (0.5, 0.62)
Jigawa | 0.54 | (0.48, 0.6) | 0.53 | (0.48, 0.59) | 0.53 | (0.47, 0.59)
Kwara | 0.51 | (0.35, 0.67) | 0.55 | (0.42, 0.68) | 0.54 | (0.43, 0.66)
Ogun | 0.51 | (0.4, 0.62) | 0.55 | (0.45, 0.65) | 0.55 | (0.45, 0.66)
Borno | 0.49 | (0.42, 0.57) | 0.49 | (0.42, 0.56) | 0.49 | (0.41, 0.57)
Yobe | 0.45 | (0.4, 0.5) | 0.45 | (0.4, 0.5) | 0.45 | (0.39, 0.51)
Kaduna | 0.43 | (0.35, 0.5) | 0.43 | (0.36, 0.5) | 0.43 | (0.36, 0.5)
Kogi | 0.42 | (0.32, 0.53) | 0.49 | (0.39, 0.59) | 0.50 | (0.39, 0.6)
Taraba | 0.42 | (0.36, 0.48) | 0.43 | (0.37, 0.49) | 0.43 | (0.36, 0.5)
Niger | 0.39 | (0.27, 0.51) | 0.40 | (0.3, 0.51) | 0.41 | (0.32, 0.5)
Bauchi | 0.36 | (0.3, 0.43) | 0.37 | (0.31, 0.43) | 0.37 | (0.31, 0.44)
Katsina | 0.34 | (0.26, 0.41) | 0.34 | (0.27, 0.41) | 0.34 | (0.27, 0.41)
Kebbi | 0.31 | (0.25, 0.38) | 0.31 | (0.24, 0.37) | 0.31 | (0.24, 0.37)
Gombe | 0.28 | (0.2, 0.36) | 0.31 | (0.23, 0.38) | 0.31 | (0.24, 0.39)
Sokoto | 0.18 | (0.13, 0.23) | 0.18 | (0.13, 0.23) | 0.17 | (0.12, 0.23)
Zamfara | 0.14 | (0.07, 0.21) | 0.17 | (0.12, 0.23) | 0.18 | (0.12, 0.24)
Table 5: Point estimates of measles vaccination rates and 90% interval estimates for Admin-1 areas among children aged 12-23 months in Nigeria in 2018 State | Hájek | Spatial Unmatched MS | Spatial Unmatched JS
---|---|---|---
Mulanje | 0.27 | (0.23, 0.32) | 0.25 | (0.21, 0.29) | 0.25 | (0.21, 0.29)
Blantyre | 0.24 | (0.2, 0.28) | 0.21 | (0.17, 0.26) | 0.22 | (0.18, 0.26)
Phalombe | 0.24 | (0.19, 0.29) | 0.23 | (0.19, 0.27) | 0.23 | (0.19, 0.28)
Zomba | 0.20 | (0.16, 0.24) | 0.19 | (0.16, 0.22) | 0.19 | (0.16, 0.22)
Ntcheu | 0.16 | (0.12, 0.21) | 0.15 | (0.11, 0.19) | 0.16 | (0.13, 0.2)
Mangochi | 0.16 | (0.12, 0.19) | 0.15 | (0.12, 0.18) | 0.15 | (0.12, 0.18)
Nsanje | 0.16 | (0.11, 0.2) | 0.15 | (0.11, 0.19) | 0.16 | (0.12, 0.2)
Thyolo | 0.15 | (0.11, 0.2) | 0.16 | (0.13, 0.19) | 0.16 | (0.13, 0.2)
Balaka | 0.15 | (0.11, 0.19) | 0.15 | (0.12, 0.18) | 0.15 | (0.12, 0.19)
Chiradzulu | 0.15 | (0.1, 0.19) | 0.16 | (0.12, 0.19) | 0.16 | (0.12, 0.2)
Neno | 0.14 | (0.1, 0.18) | 0.15 | (0.11, 0.18) | 0.15 | (0.12, 0.18)
Mwanza | 0.13 | (0.09, 0.17) | 0.14 | (0.1, 0.17) | 0.14 | (0.1, 0.17)
Karonga | 0.12 | (0.09, 0.15) | 0.11 | (0.08, 0.13) | 0.10 | (0.07, 0.13)
Chikwawa | 0.11 | (0.07, 0.15) | 0.13 | (0.1, 0.16) | 0.13 | (0.1, 0.16)
Nkhata Bay | 0.10 | (0.07, 0.14) | 0.09 | (0.06, 0.12) | 0.09 | (0.07, 0.12)
Nkhotakota | 0.10 | (0.07, 0.13) | 0.09 | (0.07, 0.11) | 0.09 | (0.07, 0.11)
Machinga | 0.09 | (0.07, 0.12) | 0.11 | (0.08, 0.13) | 0.10 | (0.08, 0.13)
Kasungu | 0.09 | (0.06, 0.12) | 0.08 | (0.06, 0.11) | 0.09 | (0.07, 0.11)
Rumphi | 0.09 | (0.05, 0.12) | 0.08 | (0.06, 0.11) | 0.09 | (0.06, 0.11)
Lilongwe | 0.09 | (0.06, 0.11) | 0.08 | (0.06, 0.11) | 0.08 | (0.06, 0.11)
Dedza | 0.08 | (0.04, 0.12) | 0.09 | (0.07, 0.12) | 0.10 | (0.08, 0.13)
Mchinji | 0.08 | (0.06, 0.1) | 0.08 | (0.06, 0.1) | 0.07 | (0.05, 0.09)
Salima | 0.07 | (0.05, 0.1) | 0.08 | (0.06, 0.1) | 0.07 | (0.06, 0.1)
Ntchisi | 0.07 | (0.05, 0.1) | 0.07 | (0.06, 0.09) | 0.07 | (0.05, 0.09)
Dowa | 0.07 | (0.04, 0.1) | 0.07 | (0.05, 0.09) | 0.07 | (0.05, 0.09)
Chitipa | 0.06 | (0.03, 0.09) | 0.07 | (0.05, 0.09) | 0.06 | (0.04, 0.09)
Mzimba | 0.06 | (0.04, 0.08) | 0.06 | (0.05, 0.08) | 0.06 | (0.05, 0.08)
Table 6: Point estimates of HIV prevalence and 90% interval estimates for
Admin-1 areas among women aged 15-29 in Nigeria in 2018