# Accounting for survey design in Bayesian disaggregation of survey-based
areal estimates of proportions: an application to the American Community
Survey

Marco H<EMAIL_ADDRESS>[ Veronica
J<EMAIL_ADDRESS>[ Roderick J.
<EMAIL_ADDRESS>[ Nationwide Children’s Hospital
Center for Injury Research and Policy
575 Children’s Crossroad
Columbus, OH 43205
Department of Statistics
School of Information and Computer Sciences
Donald Bren Hall
University of California, Irvine
Irvine, CA 92697
Department of Biostatistics
School of Public Health
1415 Washington Heights
University of Michigan
Ann Arbor, MI 48109

###### Abstract

Understanding the effects of social determinants of health on health outcomes
requires data on characteristics of the neighborhoods in which subjects live.
However, estimates of these characteristics are often aggregated over space
and time in a fashion that diminishes their utility. Take, for example,
estimates from the American Community Survey (ACS), a multi-year nationwide
survey administered by the U.S. Census Bureau: estimates for small municipal
areas are aggregated over 5-year periods, whereas 1-year estimates are only
available for municipal areas with populations $>$65,000. Researchers may wish
to use ACS estimates in studies of population health to characterize
neighborhood-level exposures. However, 5-year estimates may not properly
characterize temporal changes or align temporally with other data in the
study, while the coarse spatial resolution of the 1-year estimates diminishes
their utility in characterizing neighborhood exposure. To circumvent this
issue, in this paper we propose a modeling framework to disaggregate estimates
of proportions derived from sampling surveys which explicitly accounts for the
survey design effect. We illustrate the utility of our model by applying it to
the ACS data, generating estimates of poverty for the state of Michigan at
fine spatio-temporal resolution.

Spatio-temporal change of support problem,

Bayesian hierarchical model,

Multi-resolution approximation,

Latent spatio-temporal process,

American Community Survey,

Survey-based estimates,

###### keywords:

, and

t1Post-doctoral Scientist t2Associate Professor t3Richard D. Remington
Distinguished University Professor

##  1 Introduction

Interest and attention in the social determinants of health, that is, the
social and economic factors that characterize where and how people live, have
soared in the last 20 years (Braverman, Egerter and Williams, 2011; Marmot et
al., 2012) as awareness of health disparities within countries’ populations
has become more prevalent. Discussions on social determinants of health have
also been at the forefront of national news (TV, newspapers, magazines, etc.)
during the first months of the current COVID-19 pandemic, as various social
determinants of health – poverty, homelessness, smoke exposure, etc. – are
suspected to worsen COVID-19 outcomes (Abrams and Szefler, 2020; Rollston and
Galea, 2020; Singu et al., 2020).

A public source of information on social determinants of health is the
American Community Survey (ACS), a multi-year national survey administered by
the United States Census Bureau. Sampling annually approximately 3.5 million
Americans, including those residing in unincorporated territories (U.S. Census
Bureau, 2008), the ACS releases every year up-to-date, timely, and accurate
population and housing information to the general public and to data-users.
Due to privacy concerns and sample size limitations, often these estimates are
aggregated over space and/or time. Currently, ACS estimates for small
municipal sub-divisions, such as census tracts, are aggregated over 5-year
time periods, whereas estimates of neighborhood characteristics corresponding
to 1-year time periods are only available for municipal sub-divisions with
populations greater than 65,000.

This aggregation, if justified by privacy and statistical considerations, can
result in estimates whose spatial and/or temporal resolution is misaligned
with the target spatial and/or temporal resolution of a research study. As an
example, a researcher who wishes to incorporate an ACS estimate of poverty
(e.g. proportion of households living in poverty) in an epidemiological
analysis is typically faced with a choice: (a) utilize estimates with fine
spatial resolution whose 5-year temporal resolution is unlikely to conform to
other data sources and, in the case of longitudinal studies, fail to properly
characterize yearly changes; or (b) utilize 1-year estimates, whose
aggregation over large areal units diminishes their ability to characterize
neighborhoods in a meaningful way. Having access to estimates at fine spatial
and temporal resolution would eliminate these problems.

This is the goal of our paper. Taking the ACS as a case study, we propose a
Bayesian hierarchical spatio-temporal model that aims to generate estimates of
certain social indicators at fine spatial and temporal resolution starting
from estimates - the ACS estimates - that are either available at fine
resolution in space but not in time, or are temporally resolved but not in
space. Hence, we offer a model that solves the so-called spatio-temporal
_change of support problem_ (COSP), that is, the problem of performing
inference about a spatial or spatio-temporal process at a resolution (or
support) that differs from that of the data, in the case of multi-year
estimates of proportions derived from a complex survey. The COSP is one of the
most common problem in spatial statistics, and reviews of methods to address
it can be found in Banerjee, Carlin and Gelfand (2004) and Gotway and Young
(2002). Gelfand, Zhu and Carlin (2001) offer an extension to the space-time
setting. Our model is not the first attempt at solving the COSP for ACS data,
nor it is the first paper that models these data spatially: Bradley, Wikle and
Holan (2015), Bradley, Wikle and Holan (2016), Bradley, Holan and Wikle
(2016), Savitsky (2016), and Simpson et al. (2019) have all contributed to
this literature. In particular, Bradley, Wikle and Holan (2015) were the first
to present a statistical model that derives estimates of a socio-economic
indicator at a different spatial support than that of the ACS data, namely
over three different Native American reservations.

Our model differs from previous work in several ways. First, it deals with
spatio-temporal estimates of proportions: previous efforts considered either
variables that could be modeled using a Gaussian distribution or dealt with
estimates that referred to counts, and thus could be modeled as Poisson random
variables. As we show in Section 5.3, application and adaptation of the
aforementioned models to handle proportions, while they yield point-level
estimates that are for the most part in agreement with the ACS estimates, tend
to underrepresent the estimates’ uncertainty, leading to credible intervals
that, when validated with hold-out data, do not achieve the correct nominal
coverage.

Another key distinction of our modeling approach is that it explicitly
accounts for the survey design effect, thus merging survey methodology with
spatial statistical modeling frameworks. Although Bradley, Wikle and Holan
(2016) did account for the sampling design when modeling ACS estimates for a
given year, they only did so when specifying a COSP model for estimates of
count data: in that case, they provided a model for both the ACS estimates and
the ACS sampling-based variance, leveraging the known relationship between the
mean and the variance of a Poisson distribution. No model for the ACS
sampling-based variance was formulated in the case of Gaussian-distributed
indicators (see Bradley, Wikle and Holan (2015)): rather, the ACS variance was
taken as known and used as the variance of the normal likelihood.

Our model proposes to account for the sampling design in two ways: first, by
including the design effect (Kish, 1965), building upon the work of Korn and
Graubard (1998), Ghitza and Gelman (2013), Mercer et al. (2014), and Chen,
Wakefield and Lumley (2014), secondly by introducing random effects specified
at the spatial resolution of the sampling frame. Specifically, using both the
ACS estimates of proportions and their sampling based variance, we create two
working variables - the _effective number of cases_ (ENC) and the _effective
sample size_ (ESS) - which we use in a Binomial likelihood. Furthermore, since
the ACS sampling design uses counties as sampling frames, to account for the
fact that estimates relative to administrative areal units within the same
county might be more strongly correlated than estimates relative to areal
units that are spatially close but within different counties, our model
introduces county-level random effects.

As in Bradley, Wikle and Holan (2015) and Bradley, Wikle and Holan (2016), we
handle the COSP by assuming that the true area-level proportions result from
the aggregation of an underlying, point-referenced spatio-temporal process
over the specified area. As in Bradley, Wikle and Holan (2015) and Bradley,
Wikle and Holan (2016), such specification allows us to derive estimates over
spatio-temporal resolutions that are equal or larger than the smallest spatial
and temporal resolution for which we have data. In our application, we focus
on generating estimates at the 1-year time scale and at census tract level,
but our modeling framework could be applied to generate estimates over any
type of areal unit. To handle the large number of areal units for which we
have data, another contribution of the paper is to introduce an approximation
that alleviates computation when trying to infer upon a point-referenced,
spatio-temporal process. The approximation, called the Spatio-Temporal Multi-
Resolution Approximation (ST-MRA), is achieved through a novel basis function
expansion, which builds upon the Multi-Resolution Approximation (MRA) of a
Gaussian process presented by Katzfuss (2017).

In its focus on yielding estimates of socio-economic indicators over areal
units, our model shares similarities with other efforts within the rich small-
area estimation (SAE) literature. In particular, of the two broad classes of
methods within SAE (Pfeffermann, 2013), our model fits within the class of
model-based methods. The latter comprises statistical approaches where a
stochastic formulation is offered for the sample data, and optimal predictors,
or approximately optimal predictors, are used to derive estimates of the
quantity of interest. In using the ACS estimates as data and in specifying a
hierarchical model, we follow the same approach as Fay and Herriot (1979),
however, differently from the latter, we account directly for the spatial
dependence in the estimates. Including spatial random effects into SAE model-
based methods is not unheard of: Singh, Shukla and Kundu (2005), Pratesi and
Salvati (2008), Pereira and Coelho (2010), and Porter et al. (2014), to name a
few, have all explicitly accounted for spatial correlation in the estimates.
However, differently from us, these models do not adjust for the sampling
design, nor do they explicitly address the change of support problem in multi-
year survey estimates, which is the raison d’être of our modeling effort.

We apply our model to ACS multi-year estimates of the proportion of families
in Michigan living in poverty, and we show the ability of our model to
generate estimates with high precision, highlighting the potential for this
model to become a tool that can be used by epidemiological researchers to
derive reliable, fine-scale estimates of socioeconomic indicators. These
estimates can be subsequently incorporated into health studies examining the
role of social determinants of health on various health outcomes.

The remainder of this paper is organized as follows. Section 2 provides more
detailed background information on the ACS. Sections 3.1 to 3.5 describe our
modeling framework whereas Section 3.6 provides a succinct description of
alternative models that we apply to survey-based estimates of proportions.
Section 4 illustrates the capabilities of our model in simulation experiments,
while Section 5 presents results for the proportion of families living in
poverty in Michigan from 2006 through 2016. In both cases, the predictive
performance of our model is compared to that of alternative models. The paper
concludes with a discussion in Section 6.

##  2 Data

In this section, we provide general information on the American Community
Survey (ACS) and we present results of an exploratory data analysis performed
on the ACS estimates of the proportion of families living in poverty in
Michigan between 2006 and 2016.

###  2.1 The American Community Survey

The American Community Survey is an ongoing survey conducted by the U.S.
Census Bureau (U.S. Census Bureau, 2008, 2014). It replaced the Census long
form in the 2000 Census. It samples approximately 3.5 million households
annually, collecting data on social, housing, economic, and other community
characteristics. In contrast to the Census long form, for which data were
gathered every 10 years, the ACS surveys are administered continuously,
allowing for the timely dissemination of up-to-date community information that
are statistically representative of the time period during which the surveys
were administered.

A comprehensive report on the ACS sampling methodology is available in (U.S.
Census Bureau, 2014). Here we provide a brief overview and focus on the
sampling of housing units rather than group quarters (e.g. college dormitories
or correctional facilities). The ACS sampling procedure is broken up into two
phases: the first phase consists of the initial sample selection while the
second phase deals with follow-up surveys being sent to unmailable and non-
responding addresses. Housing units are sampled into the ACS independently for
each county in the US. To ensure that no household is selected for the ACS
more than once in a 5-year period, the sampling frame within each county is
subdivided into five disjoint sub-frames, which are rotated through every five
years. For example, ACS surveys from 2006, 2011, and 2016 are all selected
from the same sub-frame. Each year, the first phase of the ACS sampling begins
by sorting any new housing units into one of the five sub-frames.

The ACS sampling rate varies depending on the characteristics of the
neighborhood in which a housing unit resides. Housing units belong to several
municipal sub-divisions of varying sizes, or sampling entities, for example,
the unit’s city, census tract, or school district. Each of these sampling
entities is provided with a measure of size (MOS), which is approximated based
on the number of addresses contained within the entity. Blocks of housing
units are stratified based on the MOS of the smallest sampling entity that
contains that unit, which is referred to as the units’ smallest entity’s
measure of size (SEMOS). The sampling rates for the ACS are inversely
proportional to the housing units’ SEMOS. Tables 4-1 and 4-2 in (U.S. Census
Bureau, 2014) provide details on the ACS sampling rates. Once the initial
sample is selected, each address is assigned a month in which it will receive
the survey.

In the second phase of sampling, follow-up surveys are sent to a set of
randomly selected, non-responding households with higher sampling fractions
for populations with high rates of non-response.

Much like sample selection, computation of the ACS sampling weights takes
place in several stages. The first stage provides a housing unit with a so-
called basic sampling weight, which is inversely related to the unit’s
probability of selection. A series of additional calibrations then occurs,
including adjustments to ensure that the weighted estimates derived from the
ACS conform to the Census Bureau’s Population Estimates Program (PEP).
Weighted estimates of neighborhood characteristics are weighted functions of
survey responses within a neighborhood and time period. Margins of error of
are computed using successive differences replication (U.S. Census Bureau,
2014).

Given statistical accuracy, precision and privacy concerns, ACS estimates are
released with varying spatial and temporal resolution. Specifically, estimates
for small municipal subdivisions, such as census tracts, are aggregated and
provided in the form of averages over a 5-year time period, whereas yearly
estimates are provided for administrative regions that have over 65,000
inhabitants. While certain counties meet this criterion, a sizeable number of
counties in the US have less than 65,000 residents and are therefore excluded
from a dataset with a 1-year temporal resolution ACS estimates. An alternative
to using county-level estimates is to use 1-year estimates at the Public Use
Microdata Areas (PUMA) level, that is, collections of contiguous counties
and/or census tracts whose total population exceeds 100,000 people.

###  2.2 Families in poverty in Michigan

The proportion of families in poverty in an area is one of the indicators that
the US Census Bureau employs to measure poverty in the population. A family is
deemed to live in poverty if the total income of all the family members living
together is lower than a predetermined threshold. There are multiple poverty
thresholds (now, a total of 48) depending on the size of the family and the
age of the family members. Thresholds do not vary geographically across the
U.S. but are updated annually for inflation.

In this paper, we consider data on the proportion of families living in
poverty in Michigan in the period 2006-2016. Specifically, we will utilize
1-year PUMA level estimates (for a total of 68 PUMAs in Michigan) and 5-year
census tract estimates. We will use both sets of estimates to derive census
tract-level, 1-year estimates of the proportion of families living in poverty
in Michigan for every year from 2006 to 2016. Of the 2,813 census tracts in
Michigan, 84 (3%) did not have enough data to provide estimates due to a low
number of residential buildings. Hence, these census tracts were not
considered in the analysis.

Our exploratory data analysis started with an inspection of the 1-year
estimates at the PUMA level, which showed considerable spatial variability in
poverty across Michigan. While in some PUMA’s only 1.4% of the families lived
in poverty, in others that percentage raised to about 45%. However, when
averaged across Michigan, the average proportion of families living in poverty
in Michigan’s PUMAs varied between 12.1% and 12.8% in the period 2006-2016.
Investigating whether the level of poverty changed over time, we fitted a
linear mixed model to the entire times series of estimates. We considered both
a model with a linear time trend and a model with linear and quadratic terms
of time. In both models, the temporal correlation was accounted for through
the inclusion of PUMA-specific intercepts which were assumed to be
independent, identically distributed and following a common normal
distribution. The model with the quadratic trend fitted the data better and
indicated, on average, a growing level of poverty among families in Michigan
from 2006 until 2011 followed by a gradual decline.

As the ACS estimates of the proportion of families in poverty are multi-year
estimates, another goal of our exploratory data analysis was to investigate
the type of spatio-temporal dependence in the data. As we discuss in Section
3.2, our model assumes an underlying, continuous in space, discrete in time
spatio-temporal process driving the true areal proportions. Thus, to examine
the nature of the spatio-temporal dependence, we took the centroids of the
Michigan PUMA’s as observation sites, treated the data as geostatistical data,
and we used two approaches: (i) we conducted a formal Likelihood Ratio Test
(LRT) to assess separability of space and time; and (ii) we performed a more
exploratory investigation based on comparing yearly variograms. As the
boundaries of the PUMAs in Michigan changed following the 2010 Census, in
assessing space-time separability, we split the 2006-2016 data into two sets:
one consisting of data relative to the 2006-2011 pre-boundary-changes period,
and one made of estimates relative to the 2012-2016 time period. Working on
the log scale, and performing the test of separability proposed by Mitchell,
Genton and Gumpertz (2005) on the two sets of data individually, we obtained
LRT values of $7.4\times 10^{-8}$ and $2.4\times 10^{-5}$, respectively,
suggesting time and space separability. We reached a similar conclusion when
comparing the empirical semi-variograms and the associated parameters, derived
using the log of the ACS estimates of the proportion of families in poverty
for each year. Despite some annual variation, the estimates of the marginal
variance and decay parameter were generally similar over time.

In light of these results, available in the Supplementary Material, when
modeling the underlying process driving the true, areal proportion of families
living in poverty in Michigan, we decided to adopt a separable space-time
covariance function.

##  3 Modeling Approach

Our model uses both the 5-year ACS estimates at census tract level, and the
1-year ACS estimates at the PUMA level. Following Bradley, Wikle and Holan
(2015), we denote by $z_{t}^{(l)}(A)$ the estimate of a proportion
corresponding to areal unit $A$ for the $l$-year time period ending in time
$t$. Thus, $z_{t}^{(5)}(A_{ig})$ indicates the ACS estimate for the 5-year
time period ending at year $t$ for census tract $g$, $g=1,\ldots,G_{i}$,
within PUMA $i$, $i=1,\ldots,N$, whereas $z_{t}^{(1)}(A_{i})$ refers to the
1-year ACS estimate for PUMA $i$ at year $t$. We denote by
$\tau^{2(5)}_{t}(A_{ig})$ and $\tau^{2(1)}_{t}(A_{i})$ the design-based
variance of $z_{t}^{(5)}(A_{ig})$ and $z_{t}^{(1)}(A_{i})$, respectively,
derived from the margins of error provided in the ACS dataset.

###  3.1 Modeling survey-based estimates of areal proportions accounting for
the design effect

Following Bradley, Wikle and Holan (2015), we assume that the survey-based
estimate of the proportion corresponding to areal unit $A$ over the $l$-unit
time period ($l=1$ or $5$), ending at time $t$, $z_{t}^{(l)}(A)$, is related
to the true proportion, $\pi_{t}^{(l)}(A)$, through some distribution
function.

A first idea would be to model the ACS estimate $z_{t}^{(l)}(A)$ as following
a normal distribution with mean equal to the true proportion,
$\pi_{t}^{(l)}(A)$, and variance equal to the design-based variance
$\tau^{2(l)}_{t}(A)$. However, as also noted in another context by Chen,
Wakefield and Lumley (2014), such modeling choice would be inaccurate for
small samples and will not ensure that the estimated $\pi_{t}^{(l)}(A)$
belongs to the interval $[0,1]$. For this reason, building upon the work of
Korn and Graubard (1998), and following Mercer et al. (2014) and Chen,
Wakefield and Lumley (2014), we introduce a working likelihood for a random
variable $q^{*(l)}_{t}(A)$ that we construct from the ACS estimate
$z_{t}^{(l)}(A)$ and from the effective sample size $m^{*(l)}_{t}(A)$. The
latter represents the sample size that a simple random sample (SRS) should
have to yield an estimator for the proportion that has a variance that matches
the design-based variance of the ACS estimate. To derive the effective sample
size, we use the notion of design effect $d$ introduced by Kish (1995), who
calls a survey’s design effect the ratio of the variance of an estimator under
SRS to the sampling-based variance of a survey-based estimator. By setting the
design effect equal to 1 and solving the equation for the SRS sample size
$m_{t}^{(l)}(A)$, we obtain the sample size of the SRS that will yield an
estimator with variance corresponding to the ACS design-based variance. We
call this sample size the effective sample size (ESS). Including the ESS in
the distribution function that relates $q^{*(l)}_{t}(A)$ to the true
proportion $\pi_{t}^{(l)}(A)$ allows us to account for the survey’s design
effect in our modeling framework.

More specifically, in the case of $z_{t}^{(l)}(A)$, for a SRS of size
$m_{t}^{(l)}(A)$, the estimated variance of $z_{t}^{(l)}(A)$ would be equal to
$\frac{z_{t}^{(l)}(A)(1-z_{t}^{(l)}(A))}{m_{t}^{(l)}(A)}$. Setting the
survey’s design effect $d$ equal to 1, yields the following equation

$\tau^{2(l)}_{t}(A)=\frac{z_{t}^{(l)}(A)(1-z_{t}^{(l)}(A))}{m_{t}^{(l)}(A)},$

leading to the following expression for the effective sample size,
$m_{t}^{*(l)}(A)$:

$m_{t}^{*(l)}(A)=\left[\frac{z_{t}^{(l)}(A)(1-z_{t}^{(l)}(A))}{\tau^{2(l)}_{t}(A)}\right],$
(3.1)

with $\left[\cdot\right]$ denoting rounding to the nearest integer. Although
not necessarily needed in (3.1), we introduce rounding to ensure that the
effective sample size is an integer.

We then use the effective sample size $m_{t}^{*(l)}(A)$ and the ACS estimate
$z_{t}^{(l)}(A)$ for areal unit $A$ and for the $l$-year time period ending in
year $t$ to derive the effective number of cases, $q_{t}^{*(l)}(A)$, for the
same areal unit and for the same time period, that is:

$q_{t}^{*(l)}(A):=\left[m_{t}^{*(l)}(A)\cdot z_{t}^{(l)}(A)\right],$ (3.2)

again rounded to the nearest integer. This quantity represents the number of
cases that we would have observed in a SRS of size $m_{t}^{*(l)}(A)$ to obtain
an estimate of the proportion $\pi_{t}^{(l)}(A)$ that is equal to the ACS
estimate $z_{t}^{(l)}(A)$ and with the same variance as the design-based
variance $\tau^{2(l)}_{t}(A)$.

Using now the effective number of cases in our working likelihood, our
Bayesian hierarchical model specifies at the first stage a Binomial likelihood
for $q_{t}^{*(l)}(A)$ with number of trials equal to $m_{t}^{*(l)}(A)$ and
success probability equal to the true proportion, $\pi_{t}^{(l)}(A)$, our
parameter of interest, i.e.:

$q_{t}^{*(l)}(A)|\;\pi_{t}^{(l)}(A)\sim\text{Binomial}\left(m_{t}^{*(l)}(A),\pi_{t}^{(l)}(A)\right).$

As we fit our model to 5-year and 1-year ACS estimates of areal proportions,
from (3.1) and (3.2), we derive the corresponding number of cases and
effective sample sizes - $q_{t}^{(5)}(A_{ig}),m_{t}^{*(5)}(A_{ig})$ and
$q_{t}^{(1)}(A_{i}),m_{t}^{*(1)}(A_{i})$ \- which we employ in our working
likelihood, made of the following two components

$\displaystyle q_{t}^{*(5)}(A_{ig})|\;\pi_{t}^{(5)}(A_{ig})$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$
$\displaystyle\text{Binomial}\left(m_{t}^{*(5)}(A_{ig}),\pi_{t}^{(5)}(A_{ig})\right)$
$\displaystyle q_{t}^{*(1)}(A_{i})|\;\pi_{t}^{(1)}(A_{i})$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$
$\displaystyle\text{Binomial}\left(m_{t}^{*(1)}(A_{i}),\pi_{t}^{(1)}(A_{i})\right)$
(3.3)

with $g=1,\ldots,G_{i}$ and $i=1,\ldots,N$.

We note that in (3.3) we are following the tradition of spatial generalized
linear models (see Diggle, Tawn and Moyeed (1998) for details) where spatial
dependence in the data is accounted for by assuming that the model parameters
are spatially correlated.

To disaggregate the ACS estimates, we link $\pi_{t}^{(5)}(A_{ig})$ and
$\pi_{t}^{(1)}(A_{i})$ to $\pi_{t}^{(1)}(A_{ig})$,
$g=1,\ldots,G_{i};i=1,\ldots,N$, the true proportions at our desired spatial
and temporal resolution, via:

$\displaystyle\pi_{t}^{(5)}(A_{ig})$ $\displaystyle=$
$\displaystyle\frac{1}{5}\sum_{k=t-4}^{t}\pi_{k}^{(1)}(A_{ig})$
$\displaystyle\pi_{t}^{(1)}(A_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{N_{t}(A_{i})}\sum_{h=1}^{G_{i}}N_{t}(A_{ih})\pi_{t}^{(1)}(A_{ih})$
(3.4)

where $N_{t}(A)$ generally denotes the number of households in areal unit $A$
at time $t$.

###  3.2 Addressing the Change of Support Problem (COSP)

In practice, we may want to infer about proportions over areal units that are
not conveniently comprised of combinations of $A_{ig}$ and/or $A_{i}$. To this
end, we further decompose $\pi_{t}^{(1)}(A_{ig})$, i.e. the true proportion at
one-year and census tract resolution. Following in the tradition of models
handling the spatial and spatio-temporal COSP, we assume that a random
variable for an areal unit and a time period can be expressed as the
aggregation over the areal unit and the time period of a point-referenced
spatio-temporal process, continuous in space and discrete in time. Because the
process is discrete in time:

$\pi_{t}^{(5)}(A_{ig})=\frac{1}{l}\sum_{k=t-l+1}^{t}\pi_{k}^{(1)}(A_{ig})\qquad
g=1,\ldots,G_{i};\;i=1,\ldots,N.$

To allow the flexibility to work over any areal unit, we link
$\pi_{t}^{(1)}(A_{ig})$ to an underlying point-referenced spatio-temporal
process $\zeta_{t}(\mathbf{s})$, $\mathbf{s}\in\mathcal{S}$, via the probit
link function, $\Phi^{-1}(\cdot)$, thus yielding

$\Phi^{-1}\left(\pi_{t}^{(1)}(A_{ig})\right)=\frac{1}{|A_{ig}|}\int_{\mathbf{s}\in
A_{ig}}\zeta_{t}(\mathbf{s})d\mathbf{s}+\xi(C_{A_{ig}})+\epsilon_{t}(A_{ig}),$
(3.5)

with $\epsilon_{t}(A_{ig})\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\tau^{2}_{\epsilon})$ and
$\xi(C_{A_{ig}})\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\tau^{2}_{C})$. In
(3.5), $\epsilon_{t}(A_{ig})$ denote i.i.d. error terms that account for model
specification error in linking $\pi^{(1)}_{t}(A_{ig})$ to the latent process
$\zeta_{t}(\mathbf{s})$, whereas $\xi(C_{A_{ig}})$ denotes a random effect
defined at the same areal unit level as the clustering units of the sampling
survey. In (3.5), $C_{A_{ig}}$ indicates the cluster that contains areal unit
$A_{ig}$. The cluster-level random effect, $\xi(C_{A_{ig}})$, is introduced to
enforce stronger dependence among certain estimates in a way that is
reflective of the survey sampling design.

To provide an interpretation of the spatio-temporal process
$\zeta_{t}(\mathbf{s})$, $\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, in (3.5), we
consider the application to the proportion of families in poverty in any areal
unit $A$. In this case, $\zeta_{t}(\mathbf{s})$ represents a function of the
likelihood that a family living at location $\mathbf{s}\in A$ is in poverty in
year $t$. Decomposing the point-referenced spatio-temporal process
$\zeta_{t}(\mathbf{s})$, $\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, into a
large-scale spatio-temporal trend, $\mu_{t}(\mathbf{s})$, representing the
mean of the process, and a spatio-temporal random effect, $w_{t}(\mathbf{s})$,
$\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, (3.5) becomes:

$\Phi^{-1}\left(\pi_{t}^{(1)}(A_{ig})\right)=\frac{1}{|A_{ig}|}\int_{\mathbf{s}\in
A_{ig}}\left(\mu_{t}(\mathbf{s})+w_{t}(\mathbf{s})\right)d\mathbf{s}+\xi(C_{A_{ig}})+\epsilon_{t}(A_{ig}).$
(3.6)

In light of the results of our exploratory data analysis, discussed in Section
2.2, we model the spatio-temporal random effect, $w_{t}(\mathbf{s})$,
$\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, as a Gaussian spatio-temporal process
with a separable space-time covariance function with an AR(1) structure in
time and a spatial dependence encoded through the covariance function
$C(\mathbf{s},\mathbf{s}^{\prime};\boldsymbol{\theta})$,
$\mathbf{s},\mathbf{s}^{\prime}\in\mathcal{S}$.

For computations involving a large number of areal units, we approximate the
spatio-temporal process
$w_{t}(\mathbf{s}),\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, with a linear
combination of spatial basis functions with appropriate spatio-temporal basis
function weights. Given the nested geographies of the ACS, we elect to use the
basis functions implied by the Multi-Resolution Approximation (MRA; Katzfuss
(2017)), also characterized by a nested structure.

###  3.3 The Spatio-Temporal Multi-resolution Approximation (ST-MRA)

As the MRA is defined only in a spatial context, here we extend it to the
spatio-temporal setting. Let
$w_{t}(\mathbf{s}),\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, denote a mean-zero
spatio-temporal Gaussian process defined on a spatial domain $\mathcal{S}$
with a separable space-time covariance function that invokes a first-order
autoregressive structure in time and spatial covariance function
$C(\mathbf{s},\mathbf{s}^{\prime};\boldsymbol{\theta})$,
$\mathbf{s},\mathbf{s}^{\prime}\in\mathcal{S}$. As in the MRA, we start by
introducing a first set of $r$ knots on the spatial domain $\mathcal{S}$
(level 0). Then, at each level $m$ ($m=1,\ldots,M$), we recursively partition
the spatial domain $\mathcal{S}$ in $J^{m}$ non-overlapping subregions in
which we introduce $r$ knots. Let $S^{*}_{m,j}$ denote the set of $r$ knots
defined on partition $j$ of level $m$. We define the basis functions
$\mathbf{b}_{m,j}(\mathbf{s})$, for $j=1,\ldots,J^{m};m=0,\ldots,M$
recursively as:

$\displaystyle v_{0}(\mathbf{s}_{1},\mathbf{s}_{2})$ $\displaystyle=$
$\displaystyle C(\mathbf{s}_{1},\mathbf{s}_{2};\boldsymbol{\theta})$
$\displaystyle\mathbf{b}_{m,j}(\mathbf{s})$ $\displaystyle:=$ $\displaystyle
v_{m}(\mathbf{s},S^{*}_{m,j})$ $\displaystyle\mathbf{K}^{-1}_{m,j}$
$\displaystyle:=$ $\displaystyle v_{m}(S^{*}_{m,j},S^{*}_{m,j})$
$\displaystyle v_{m+1}(\mathbf{s}_{1},\mathbf{s}_{2})$ $\displaystyle=$
$\displaystyle 0\qquad\text{ if }\mathbf{s}_{1}\text{ and
}\mathbf{s}_{2}\text{ are in different regions at resolution $m$}$ (3.7)
$\displaystyle v_{m+1}(\mathbf{s}_{1},\mathbf{s}_{2})$ $\displaystyle:=$
$\displaystyle
v_{m}(\mathbf{s}_{1},\mathbf{s}_{2})-\mathbf{b}_{m,j}(\mathbf{s}_{1})^{\prime}\mathbf{K}_{m,j}\mathbf{b}_{m,j}(\mathbf{s}_{2})\qquad\text{
otherwise.}$

In the MRA construction (Katzfuss, 2017), the basis functions weights
$\boldsymbol{\eta}_{m,j}$ are specified to follow a multivariate normal
distribution $\boldsymbol{\eta}_{m,j}\sim N_{r}(\mathbf{0},\mathbf{K}_{m,j})$.
With this specification for the basis functions and the basis functions
weights, the linear combination
$\sum_{m=0}^{M}\sum_{j=1}^{J^{m}}\mathbf{b}_{m,j}(\mathbf{s})\boldsymbol{\eta}_{m,j}$
yields an $M$-level approximation to a mean-zero Gaussian process with
covariance function $C(\mathbf{s},\mathbf{s}^{\prime};\boldsymbol{\theta})$.

For our spatio-temporal process $w_{t}(\mathbf{s})$,
$\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, we let the basis function weights
$\boldsymbol{\eta}_{t,m,j}$ vary in time, modeling them with a stationary,
first-order autoregressive structure (Gelfand, Banerjee and Gamerman, 2005).
Hence, at time $t=1$ we assume that $\boldsymbol{\eta}_{1,m,j}\sim
N_{r}(\mathbf{0},\mathbf{K}_{m,j})$, while for $t=2,\dots,T$:

$\displaystyle\boldsymbol{\eta}_{t,m,j}|\boldsymbol{\eta}_{t-1,m,j},\boldsymbol{\eta}_{t-2,m,j},\ldots,\boldsymbol{\eta}_{1,m,j}$
$\displaystyle\sim$ $\displaystyle
N_{r}(\alpha\boldsymbol{\eta}_{t-1,m,j},\mathbf{U}_{m,j}),$ (3.8)
$\displaystyle\mathbf{U}_{m,j}$ $\displaystyle=$
$\displaystyle(1-\alpha^{2})\mathbf{K}_{m,j}.$

We call

$w_{t,M}(\mathbf{s}):=\sum_{m=0}^{M}\sum_{j=1}^{J^{m}}\mathbf{b}_{m,j}(\mathbf{s})\boldsymbol{\eta}_{t,m,j}$
(3.9)

with basis functions $\mathbf{b}_{m,j}(\mathbf{s})$ defined as in (3.7), the
M-level ST-MRA approximation of the separable, spatio-temporal process
$w_{t}(\mathbf{s})$, $\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, with AR(1)
dependence in time and spatial covariance function
$C(\mathbf{s},\mathbf{s}^{\prime};\boldsymbol{\theta})$. Section 1 of the
Supplementary Material (Benedetti, Berrocal and Little, 2021) shows that the
above expression does indeed provide an approximation of the desired spatio-
temporal dependence structure.

###  3.4 The Bayesian spatio-temporal disaggregation model

Combining the formulations in Sections 3.2 and 3.3, we obtain:

$\displaystyle\Phi^{-1}\left(\pi_{t}^{(1)}(A_{ig})\right)$ $\displaystyle=$
$\displaystyle\frac{1}{|A_{ig}|}\int_{\mathbf{s}\in
A_{ig}}\left(\mu_{t}(\mathbf{s})+w_{t}(\mathbf{s})\right)d\mathbf{s}+\xi(C_{A_{ig}})+\epsilon_{t}(A_{ig})$
(3.10) $\displaystyle\approx$
$\displaystyle\frac{1}{|A_{ig}|}\int_{\mathbf{s}\in
A_{ig}}\left(\mu_{t}(\mathbf{s})+w_{t,M}(\mathbf{s})\right)d\mathbf{s}+\xi(C_{A_{ig}})+\epsilon_{t}(A_{ig})$
$\displaystyle=$ $\displaystyle\frac{1}{|A_{ig}|}\int_{\mathbf{s}\in
A_{ig}}\left(\mu_{t}(\mathbf{s})+\sum_{m=0}^{M}\sum_{j=1}^{J^{m}}\mathbf{b}_{m,j}(\mathbf{s})\boldsymbol{\eta}_{t,m,j}\right)d\mathbf{s}+\xi(C_{A_{ig}})+\epsilon_{t}(A_{ig}),$

where the spatio-temporal random effect $w_{t}(\mathbf{s})$,
$\mathbf{s}\in\mathcal{S}$, has been replaced by its $M$-level ST-MRA
approximation $w_{t,M}(\mathbf{s})$ defined in (3.9).

As the sampling frames in the ACS survey consist of counties, denoting by
$C_{A_{ig}}$ the county containing census tract $A_{ig}$, our model for
disaggregating spatially and temporally the ACS estimates of proportions, has
the following hierarchical specification:

$\displaystyle q_{t}^{*(5)}(A_{ig})|\;\pi_{t}^{(5)}(A_{ig})$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$
$\displaystyle\text{Binomial}\left(m_{t}^{*(5)}(A_{ig}),\pi_{t}^{(5)}(A_{ig})\right)$
$\displaystyle q_{t}^{*(1)}(A_{i})|\;\pi_{t}^{(1)}(A_{i})$
$\displaystyle\stackrel{{\scriptstyle ind}}{{\sim}}$
$\displaystyle\text{Binomial}\left(m_{t}^{*(1)}(A_{i}),\pi_{t}^{(1)}(A_{i})\right)$
$\displaystyle\pi_{t}^{(5)}(A_{ig})$ $\displaystyle=$
$\displaystyle\frac{1}{5}\sum_{k=t-4}^{t}\pi_{k}^{(1)}(A_{ig})$ (3.11)
$\displaystyle\pi_{t}^{(1)}(A_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{N_{t}(A_{i})}\sum_{h=1}^{G_{i}}N_{t}(A_{ih})\pi_{t}^{(1)}(A_{ih})$
$\displaystyle\Phi^{-1}\left(\pi_{t}^{(1)}(A_{ig})\right)$
$\displaystyle\approx$ $\displaystyle\frac{1}{|A_{ig}|}\int_{\mathbf{s}\in
A_{ig}}\left(\mu_{t}(\mathbf{s})+\sum_{m=0}^{M}\sum_{j=1}^{J^{m}}\mathbf{b}_{m,j}(\mathbf{s})\boldsymbol{\eta}_{t,m,j}\right)d\mathbf{s}$
$\displaystyle+$ $\displaystyle\xi(C_{A_{ig}})+\epsilon_{t}(A_{ig})$
$\displaystyle\xi(C_{A_{ig}})$ $\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}$ $\displaystyle N(0,\tau_{C}^{2})$
$\displaystyle\epsilon_{t}(A_{ig})$ $\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}$ $\displaystyle N(0,\tau^{2}_{\epsilon})$

with $q_{t}^{*(5)}(A_{ig})$, $q_{t}^{*(1)}(A_{i})$, $m_{t}^{*(5)}(A_{ig})$,
and $m_{t}^{*(1)}(A_{i})$ defined as in (3.2) and (3.1), respectively. The
county-level random effects in the expression of
$\Phi^{-1}\left(\pi_{t}^{(1)}(A_{ig})\right)$ allow the ACS estimates for
census tracts within the same county to exhibit greater dependence with one
another than with ACS estimates for census tracts in different counties, even
when the distances between those tracts are the same. We speculate that this
will account for the fact that factors such as sampling procedure or response
rate within a county-wide sampling frame might systematically affect ACS
estimates corresponding to most or all of the census tracts within that
county.

The integral in (3.10) can be re-expressed as:

$\displaystyle\Phi^{-1}\left(\pi_{t}^{(1)}(A_{ig})\right)$
$\displaystyle\approx$
$\displaystyle\mu_{t}(A_{ig})+\sum_{m=0}^{M}\sum_{j=1}^{J^{m}}\mathbf{b}_{m,j}(A_{ig})\boldsymbol{\eta}_{t,m,j}+\xi(C_{A_{ig}})+\tilde{\epsilon}_{t}(A_{ig})$
$\displaystyle\tilde{\epsilon}_{t}(A_{ig})$
$\displaystyle\stackrel{{\scriptstyle iid}}{{\sim}}$ $\displaystyle
N(0,\tau^{2}),$ (3.12)

where $\mu_{t}(A_{ig})$ and $\mathbf{b}_{m,j}(A_{ig})$ denote the integrals of
$\mu_{t}(\mathbf{s})$ and of the basis functions
$\mathbf{b}_{m,j}(\mathbf{s})$, $m=0,\ldots,M$; $j=1,\ldots,J^{m}$, as
$\mathbf{s}$ varies in areal unit $A_{ig}$, with
$g=1,\ldots,G_{i};i=1,\ldots,N$. The term $\tilde{\epsilon}_{t}(A_{ig})$ in
(3.12) accounts for errors due to model misspecification, aggregation as well
as any error that occurs as a result of the multi-resolution space-time
approximation.

###  3.5 Prior distributions

Our Bayesian model includes Inverse Gamma prior distributions for the error
variance parameter $\tau^{2}$ in (3.12), and for the variance of the county-
level random effects, $\tau^{2}_{C}$. We assume
$\mu_{t}(\mathbf{s})\equiv\mu_{t},t=1,2,...,T$, and model these spatially-
constant temporal trend terms as independent a priori, with an improper prior
$p(\mu_{t})\propto 1,\forall t$. This modeling choice implies that $\forall
t=1,2,\ldots,T$, the spatio-temporal random effect
$w_{t}(\mathbf{s}),\mathbf{s}\in\mathcal{S}$, accounts for all the spatial
variation in the ACS estimates. We investigated whether allowing the mean
terms $\mu_{t}$ vary in space would lead to significantly different results in
terms of model fit, but we did not observe any meaningful change.

We assign a Uniform $([0,1])$ prior to the autoregressive parameter $\alpha$
of the basis functions weights in (3.8), while the definition of the ST-MRA
basis functions is determined once we choose the spatial covariance function
$C(\mathbf{s},\mathbf{s}^{\prime};\boldsymbol{\theta})$. Here we take it to be
the stationary Matérn covariance function with parameters $\sigma^{2}$, $\phi$
and $\nu$

$C(\mathbf{s},\mathbf{s}^{\prime};\boldsymbol{\theta})=\frac{\sigma^{2}}{2^{\nu-1}\Gamma(\nu)}\left(\frac{||\mathbf{s}-\mathbf{s}^{\prime}||}{\phi}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{||\mathbf{s}-\mathbf{s}^{\prime}||}{\phi}\right)$
(3.13)

where $\mathbf{s},\mathbf{s}^{\prime}\in\mathcal{S}$ and
$\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function of the second kind.

We specify a non-informative Inverse Gamma prior on the marginal variance
parameter $\sigma^{2}$, while we place a Gamma$(1,1)$ prior on the range
parameter $\phi$ and a Uniform$((0,2))$ prior on the smoothness parameter
$\nu$, as suggested by Finley, Banerjee and Carlin (2007). As the latter is
notoriously difficult to estimate, an alternative specification would entail
the use of penalized complexity priors as described in Simpson et al. (2017).

###  3.6 Other models

We describe succinctly alternative models that we compare with our model in
Sections 4 and 5. More details are available in Section 3 of the Supplementary
Material (Benedetti, Berrocal and Little, 2021). To evaluate the utility of
the _effective sample size_ and _effective number of cases_ , a first
competing model specifies a “standard” Binomial likelihood for the number of
cases, obtained by multiplying the ACS estimate, $z^{(l)}_{t}(A)$, by the
number of survey responses, $m_{t}^{(l)}(A)$, obtained in areal unit $A$ over
the $l$-unit time period ending in year $t$. Calling this product
$q^{(l)}_{t}(A)$, the _standard Binomial model_ for disaggregation applied to
the 1-year and 5-year ACS data assumes that

$\begin{array}[]{rcl}q_{t}^{(1)}(A_{i})|\pi^{(1)}_{t}(A_{i})&\sim&\mbox{Binomial}\left(m^{(1)}_{t}(A_{i}),\pi^{(1)}_{t}(A_{i})\right)\\\
q^{(5)}_{t}(A_{ig})|\pi^{(5)}_{t}(A_{ig})&\sim&\mbox{Binomial}\left(m^{(5)}_{t}(A_{ig}),\pi^{(5)}_{t}(A_{ig})\right).\end{array}$

We keep the other levels of this model exactly as in the Bayesian hierarchical
model in (3.11).

The second and third model we consider are extensions and adaptations of
models proposed by Bradley, Wikle and Holan (2016) and Bradley, Wikle and
Holan (2015), respectively, when analyzing Poisson spatial-only and Gaussian
space-time ACS data.

The _BWH Poisson space-time model_ extends the model for count data proposed
by Bradley, Wikle and Holan (2016) to the space-time setting. Interpreting the
counts $q^{(1)}_{t}(A_{i})$ and $q^{(5)}_{t}(A_{ig})$ as Poisson random
variables, we assume a latent process $Y_{t}(A_{ig})$, $t=1,\ldots,T$, defined
at the census tract level, such that

$\begin{array}[]{rcl}q^{(1)}_{t}(A_{i})|\left\\{Y_{t}(A_{ig});g=1,\ldots,G_{i},i=1,\ldots,N\right\\}&\sim&\mbox{Poisson}\left(\sum_{h=1}^{G_{i}}\exp\left(Y_{t}(A_{ih})\right)\right)\\\
q^{(5)}_{t}(A_{ig})|\left\\{Y_{k}(A_{ig});k=1,\ldots,T\right\\}&\sim&\mbox{Poisson}\left(\frac{1}{5}\sum_{k=t-4}^{t}\exp\left(Y_{k}(A_{ig})\right)\right)\end{array}$

for $g=1,\ldots,G_{i};i=1,\ldots,N$ and $t=1,\ldots,T$. For each
$t=1,\ldots,T$, following Bradley, Wikle and Holan (2016), $Y_{t}(A_{ig})$ is
decomposed as:

$Y_{t}(A_{ig})=\beta_{t}+\boldsymbol{\psi}\boldsymbol{\vartheta}+\varsigma_{t}(A_{ig})$
(3.14)

with $\boldsymbol{\psi}$ Moran’s basis functions, $\boldsymbol{\vartheta}$
basis functions weights defined as in Bradley, Wikle and Holan (2016), and
$\varsigma_{t}(A_{ig})$ error terms that account for aggregation and other
types of errors. Differently from Bradley, Wikle and Holan (2016), here we are
dealing with estimates over multiple years: to accommodate this added
dimension, in (3.14) we allow the intercept terms $\beta_{t}$ to vary in time,
hence representing a temporal trend. Similarly, we allow the error terms
$\varsigma_{t}(A_{ig})$ to change in time. Finally, as in Bradley, Wikle and
Holan (2016), the _BWH Poisson space-time model_ provides a stochastic
formulation for the ACS design based variances $\tau^{2(1)}_{t}(A_{i})$ and
$\tau^{2(5)}_{t}(A_{ig})$, assumed respectively to follow a lognormal
distribution:

$\begin{array}[]{rcl}\log\left(\tau^{2(1)}_{t}(A_{i})\right)&\sim&N\left(\log\left(\sum_{h=1}^{G_{i}}\exp\left(Y_{t}(A_{ih})\right)\right),\sigma^{2(1)}(A_{i})\right)\\\
\log\left(\tau^{2(5)}_{t}(A_{ig})\right)&\sim&N\left(\log\left(\frac{1}{5}\sum_{k=t-4}^{t}\exp\left(Y_{k}(A_{ig})\right)\right),\sigma^{2(5)}(A_{ig})\right).\end{array}$

The specification of the _BWH Poisson space-time model_ is completed by the
following prior distributions, which we take directly from Bradley, Wikle and
Holan (2016): $\beta_{t}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,10^{15}),\forall t=1,\ldots,T$;
$\varsigma_{t}(A_{ig})\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\sigma^{2}_{\varsigma}),\forall t=1,\ldots,T$,
$g=1,\ldots,G_{i},i=1,\ldots,N$; $\sigma^{2(1)}(A_{i})\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Gamma}(1,1),\forall i=1,\ldots,N$;
$\sigma^{2(5)}(A_{ig})\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Gamma}(1,1),\forall g=1,\ldots,G_{i},i=1,\ldots,N$; and
$\sigma^{2}_{\varsigma}\sim\mbox{Gamma}(1,1)$.

The last model we consider is an adaption of the spatio-temporal model
proposed by Bradley, Wikle and Holan (2015) for Gaussian-distributed ACS
variables to ACS estimates of proportions. To frame the ACS estimates of
proportions, $z^{(l)}_{t}(A)$, within a Gaussian likelihood, we apply a
logistic transformation to them, thus obtaining variables defined in
$\mathbf{R}$. As $\tau^{2(l)}_{t}(A)$ is the design-based variance of
$z^{(l)}_{t}(A)$,we employ the delta method to derive the expression of the
variance of $\log\left(\frac{z^{(l)}_{t}(A)}{(1-z^{(l)}_{t}(A))}\right)$ for
every areal unit $A$ and $l$-unit time period. Thus, the first stage of this
new Bayesian hierarchical model, which we call the _BWH Gaussian Delta method
model_ , is given by:

$\begin{array}[]{rcl}\log\left(\frac{z^{(1)}_{t}(A_{i})}{(1-z^{(1)}_{t}(A_{i}))}\right)|\pi^{(1)}_{t}(A_{i})&\sim&N\left(\log\left[\frac{\pi^{(1)}_{t}(A_{i})}{(1-\pi^{(1)}_{t}(A_{i}))}\right],\frac{\tau^{2(1)}_{t}(A_{i})}{z^{(1)}_{t}(1-z^{(1)}_{t}(A_{i}))}\right)\\\
\\\
\log\left(\frac{z^{(5)}_{t}(A_{ig})}{(1-z^{(5)}_{t}(A_{ig}))}\right)|\pi^{(5)}_{t}(A_{ig})&\sim&N\left(\log\left[\frac{\pi^{(5)}_{t}(A_{ig})}{(1-\pi^{(5)}_{t}(A_{ig}))}\right],\frac{\tau^{2(5)}_{t}(A_{ig})}{z^{(5)}_{t}(1-z^{(5)}_{t}(A_{ig}))}\right)\end{array}$

for $i=1,\ldots,N$, $g=1,\ldots,G_{i}$, $t=1,\ldots,T$.

Calling
$\tilde{y}^{(1)}_{t}(A_{i}):=\log\left[\frac{\pi^{(1)}_{t}(A_{i})}{(1-\pi^{(1)}_{t}(A_{i}))}\right]$
and
$\tilde{y}^{(5)}_{t}(A_{ig}):=\log\left[\frac{\pi^{(1)}_{t}(A_{ig})}{(1-\pi^{(1)}_{t}(A_{ig}))}\right]$,
we achieve their disaggregation in time through the following equality

$\tilde{y}^{(5)}_{t}(A_{ig})=\frac{1}{5}\sum_{k=t-4}^{t}\tilde{y}^{(1)}_{t}(A_{ig})\qquad
i=1,\ldots,N;g=1,\ldots,G_{i}$

whereas their disaggregation in space is handled, for any areal unit $A$,
through

$\tilde{y}^{(1)}_{t}(A)=\frac{1}{|A|}\int_{s\in
A}\zeta_{t}(\mathbf{s})d\mathbf{s}\qquad\forall t=1,\ldots,T$

with $\zeta_{t}(\mathbf{s})$ spatio-temporal Gaussian process. Following a
similar approach as discussed in Section 3.2, for all $t=1,\ldots,T$, we
decompose $\zeta_{t}(\mathbf{s})$ in the sum of a spatio-temporal trend term
$\mu_{t}(\mathbf{s})$ and spatio-temporal random effects $w_{t}(\mathbf{s})$,
$\mathbf{s}\in\mathcal{S};t=1,\ldots,T$, and we replace
$\zeta_{t}(\mathbf{s})$ with $\mu_{t}(\mathbf{s})+w_{t}(\mathbf{s})$ under the
integral. As in Section 3.2, $w_{t}(\mathbf{s})$ assumed to be equipped with a
separable space-time covariance function.

In fitting this model to data we use a dimension reduction approach, and we
approximate the spatio-temporal random effects $w_{t}(\mathbf{s})$,
$\mathbf{s}\in\mathcal{S};t=1,\ldots,T$ using the ST-MRA approximation
discussed in Section 3.3. Also Bradley, Wikle and Holan (2015) handled the
large dimensionality of the data through an approximation that involved a
basis function expansion. However, they used bisquare basis functions rather
than the basis functions we employ here. We believe that the difference in
basis functions employed in the approximation should not result in drastic
changes in terms of model performance.

###  3.7 Computation

We fit our Bayesian hierarchical model and all the other competing models
using Markov Chain Monte Carlo (MCMC) algorithms, with Gibbs sampling and
Metropolis-Hastings steps. For our model, posterior sampling exploits the data
augmentation method of Albert and Chib (1993) to sample the MRA basis function
coefficients $\boldsymbol{\eta}_{t,m,j}$ via Gibbs sampling, whereas posterior
samples of the Matérn covariance function parameters - $\phi$ and $\nu$ \- are
generated using a Metropolis-Hastings algorithm. We assess convergence of the
MCMC algorithms both visually, by inspecting trace plots, and numerically
using Geweke’s diagnostic for Markov chains (Geweke, 1992). We run each MCMC
algorithm for a number of iterations large enough that the effective sample
size post burn-in for each model parameter exceeds 1,000. The proposal
distributions used in the Metropolis-Hastings steps are tuned during burn-in
to achieve desirable acceptance rates (Roberts, Gelman and Gilks, 1997).

##  4 Simulation Studies

We now report results for two simulation studies: Simulation study 1 evaluates
the ability of the proposed model to disaggregate spatially and temporally
areal-level estimates of proportions, even when the data are not generated
according to our model specifications. On the other hand, Simulation study 2
gauges the need to account for the design effect.

###  4.1 Generating the true proportions

In both simulation studies, we use very similar data generating mechanisms,
all very different from our modeling framework. As a spatial domain we
envision a geographical configuration that is analogous to that of census
tracts and PUMAs. Specifically, we consider a $10\times 10$ square grid made
of 100 areal units, all assumed to have the same population size. The 100
areal units are in turn grouped into 4 distinct regions, each with the same
population size and each containing 25 areal units (see Figure 1). To simulate
data, we proceed as follows: we generate the true 1-year proportions
$\pi^{(1)}_{t}(A_{ig})$ for each time $t=1,\ldots,10$ and for each subregion
$g$, $g=1,\ldots,25$, nested within region $i$, $i=1,\ldots,4$. We repeat the
procedure thirty times, yielding a total of 30 simulated datasets per
simulation setting, and we consider 4 different data generating mechanisms.
This allows us to assess the performance of our model in settings that differ
from that of our model.

Figure 1:  Areal units utilized in the simulation studies.

Under each simulation setting we assume that in each subregion $A_{ig}$, there
is a latent covariate $x(A_{ig})$, not varying in time, distributed according
to a standard normal distribution. This latent covariate drives the true
proportion $\pi^{(1)}_{t}(A_{ig})$, $g=1,2,\ldots,25$; $i=1,\ldots,4$. In the
first simulation setting, for each subregion $A_{ig}$ and at each time point
$t$, the true proportion $\pi^{(1)}_{t}(A_{ig})$ is obtained by applying the
_expit_ (inverse logistic) function to the sum of the latent covariate
$x(A_{ig})$ and the randomly generated white noise, thus allowing for temporal
and spatial variability in the true proportions. Although the true proportions
$\pi^{(1)}_{t}(A_{ig})$ will not be the same across space and time, they are
independent in space and time.

To induce spatial correlation in the true proportions, in the second
simulation setting we introduce a point-referenced spatial process,
$\lambda(\mathbf{s})$, with a Matérn covariance function with unit marginal
variance (e.g. $\sigma^{2}_{\lambda}=1$), and range and smoothness parameters
($\phi_{\lambda}$ and $\nu_{\lambda}$, respectively) equal to 0.5 and 1. This
implies that the effective range of the spatial process $\lambda(\mathbf{s})$
is between 1 and 2 resulting in true proportions for neighboring subregions
that are spatially dependent. The true proportion, $\pi^{(1)}_{t}(A_{ig})$,
for areal unit $A_{ig}$ at time $t$ is obtained by applying the expit function
to the sum of the latent covariate $x(A_{ig})$, the spatial process
$\lambda(\mathbf{s}_{ig})$ evaluated at the centroid $\mathbf{s}_{ig}$ of
areal unit $A_{ig}$, and the white noise term $e_{t}(A_{ig})$. Although this
second data generating mechanism yields spatially correlated true proportions,
they are independent over time.

The third data generating mechanism allows for a temporal trend in the true
proportions by introducing a linear time trend $\alpha_{0}+\alpha_{1}t$. Thus,
the true proportion for areal unit $A_{ig}$ at time $t$ is now obtained by
applying the expit function to the sum of the latent covariate $x(A_{ig})$,
the linear trend, $\alpha_{0}+\alpha_{1}t$, and the white noise
$e_{t}(A_{ig})$. We employ $\alpha_{0}=-1.0$ and $\alpha_{1}=0.2$ in the
linear temporal trend, resulting in a noticeable increase over time of the
true proportions.

Despite the temporal dependence in the true proportions generated under the
third simulation setting, the $\pi^{(1)}_{t}(A_{ig})$’s are not spatially
correlated. To address this shortcoming, the fourth data generating mechanism
combines the second and third data generating mechanism together yielding true
proportions $\pi^{(1)}_{t}(A_{ig})$ that display a temporal trend and are
correlated in space. Thus, in short:

$\pi^{(1)}_{t}(A_{ig})=\frac{\exp\\{x(A_{ig})+\lambda(\mathbf{s}_{ig})+\alpha_{0}+\alpha_{1}t+e_{t}(A_{ig})\\}}{1+\exp\\{x(A_{ig})+\lambda(\mathbf{s}_{ig})+\alpha_{0}+\alpha_{1}t+e_{t}(A_{ig})\\}}\qquad
i=1,...,4;g=1,...,25\\\ $
$\begin{array}[]{lrllrrl}\mathbf{X}&=&\\{x(A_{ig})\\}_{i=1,...,4;g=1,...,25};&&x(A_{ig})&\stackrel{{\scriptstyle
iid}}{{\sim}}&N(0,1)\\\ \\\
\boldsymbol{\lambda}&=&\\{\lambda(\mathbf{s}_{ig})\\}_{i=1,...,4;g=1,...,25};&&\boldsymbol{\lambda}&\sim&\text{MVN}\left(0,\Sigma(\boldsymbol{\theta}_{\lambda})\right)\end{array}$

with $\Sigma(\boldsymbol{\theta}_{\lambda})$, 100$\times$100 covariance matrix
induced by a Matérn covariance function, e.g. by (3.13), with
$\boldsymbol{\theta}_{\lambda}=\left(\sigma^{2}_{\lambda},\phi_{\lambda},\nu_{\lambda}\right)=\left(1.0,\;0.5,\;1.0\right)^{\prime}$.

Table 1:  Data generating mechanism used in each of the four simulation settings of both simulation studies. Setting | Equation
---|---
1 | $\displaystyle\pi^{(1)}_{t}(A_{ig})=\frac{\exp\\{x(A_{ig})+e_{t}(A_{ig})\\}}{1+\exp\\{x(A_{ig})+e_{t}(A_{ig})\\}},\qquad e_{t}(A_{ig})\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.2^{2})$
2 | $\displaystyle\pi^{(1)}_{t}(A_{ig})=\frac{\exp\\{x(A_{ig})+\lambda(\mathbf{s}_{ig})+e_{t}(A_{ig})\\}}{1+\exp\\{x(A_{ig})+\lambda(\mathbf{s}_{ig})+e_{t}(A_{ig})\\}},\qquad e_{t}(A_{ig})\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.2^{2})$
3 | $\displaystyle\pi^{(1)}_{t}(A_{ig})=\frac{\exp\\{x(A_{ig})+\alpha_{0}+\alpha_{1}t+e_{t}(A_{ig})\\}}{1+\exp\\{x(A_{ig}))+\alpha_{0}+\alpha_{1}t+e_{t}(A_{ig})\\}},\qquad e_{t}(A_{ig})\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.2^{2})$
4 | $\displaystyle\pi^{(1)}_{t}(A_{ig})=\frac{\exp\\{x(A_{ig})+\lambda(\mathbf{s}_{ig})+\alpha_{0}+\alpha_{1}t+e_{t}(A_{ig})\\}}{1+\exp\\{x(A_{ig})+\lambda(\mathbf{s}_{ig})+\alpha_{0}+\alpha_{1}t+e_{t}(A_{ig})\\}},\qquad e_{t}(A_{ig})\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.2^{2})$

###  4.2 Generating the observed estimates

Having generated the true proportions $\pi^{(1)}_{t}(A_{ig})$ under the 4 data
generating mechanisms, we proceed to simulate the corresponding “ _observed_ ”
5-year and 1-year estimates, $z^{(5)}_{t}(A_{ig})$ and $z_{t}^{(1)}(A_{i})$,
respectively. These estimates play the equivalent role to the ACS estimates,
in that they represent the data to which our model is fit. They are obtained
by adding to the true proportions additional random error, which represents
the survey-based error. Specifically, we first generate the 1-year subregional
estimates $z_{t}^{(1)}(A_{ig})$ by adding error to the true proportions on the
logit scale:

$\begin{array}[]{rcl}\log\left(\frac{z^{(1)}_{t}(A_{ig})}{1-z^{(1)}_{t}(A_{ig})}\right)&=&\text{logit}\left(\pi^{(1)}_{t}(A_{ig})+\tilde{e}_{t}(A_{ig}\right)\\\
\\\ \tilde{e}_{t}(A_{ig})&\stackrel{{\scriptstyle
iid}}{{\sim}}&N(0,v_{t}(A_{ig}))\end{array}$

From these we then derive the 1-year regional and the 5-year subregional
observed estimates, $z^{(1)}_{t}(A_{i})$ and $z^{(5)}_{t}(A_{ig})$, as
follows:

$\displaystyle z_{t}^{(5)}(A_{ig})$ $\displaystyle=$
$\displaystyle\frac{1}{5}\sum_{k=t-4}^{t}z^{(1)}_{k}(A_{ig})$ $\displaystyle
z_{t}^{(1)}(A_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{25}\sum_{h=1}^{25}z^{(1)}_{t}(A_{ih}).$

We use two different strategies to determine the magnitude of the variances
$v_{t}(A_{ig})$’s. In simulation study 1, the $v_{t}(A_{ig})$’s are fixed
across the four simulation settings and are chosen so that the variation in
the simulated observed estimates at adjacent time periods resembles the year-
to-year variation in the ACS estimates of the proportion of families in
poverty. This is achieved when $v_{t}(A_{ig})=0.15^{2}$ for $t=1,2,\ldots,10$;
$i=1,\ldots,4$; $g=1,2,\ldots,25$.
In simulation study 2, we derive the magnitude of the $v_{t}(A_{ig})$’s as a
function of the design effect $d$. Since the goal of this simulation study is
to evaluate the inferential gain obtained by working with the _effective
sample size_ and _effective number of cases_ , rather than with the _observed
number of cases_ and the _observed sample size_ , we let $d$ vary. This
results in different values of the $v_{t}(A_{ig})$’s. The relationship between
the $v_{t}(A_{ig})$’s and the design effect $d$ can be determined based on the
following consideration: the $v_{t}(A_{ig})$’s ought to be such that,
conditional on the true proportions $\pi^{(1)}_{t}(A_{ig})$,

$\displaystyle\text{Var}\left(z^{(1)}_{t}(A_{ig})|\pi^{(1)}_{t}(A_{ig})\right)$
$\displaystyle=$ $\displaystyle
d\cdot\mbox{Var}\left[\text{expit}\left\\{\text{logit}\left(\pi^{(1)}_{t}(A_{ig})\right)+\tilde{e}_{t}(A_{ig})\right\\}\right]$
(4.1) $\displaystyle=$ $\displaystyle d\cdot\mbox{Var}^{(1)}_{SRS,t}(A_{ig}),$

with $\mbox{Var}^{(1)}_{SRS,t}(A_{ig})$ variance of the estimator
$\hat{\pi}^{(1)}_{SRS,t}(A_{ig})$ of the true proportion
$\pi^{(1)}_{t}(A_{ig})$ based on a simple random sample (SRS) . This leads to
the following expression for $v_{t}(A_{ig})$:

$\displaystyle v^{(1)}_{t}(A_{ig})$ $\displaystyle=$ $\displaystyle
d\times\frac{\left(\exp\left\\{\text{logit}\left(\pi^{(1)}_{t}(A_{ig})\right)\right\\}+1\right)^{4}\pi^{(1)}_{t}(A_{ig})\left(1-\pi^{(1)}_{t}(A_{ig})\right)}{m^{(1)}_{t}(A_{ig})\exp\left\\{2\times\text{logit}(\pi^{(1)}_{t}(A_{ig}))\right\\}}.$
(4.2)

Letting the sample sizes $m^{(1)}_{t}(A_{ig})$ for each subregion $A_{ig}$,
$g=1,\ldots,25;i=1,\ldots,4$, be equal to 100 for each time $t$, we derive
from (4.2) the values of the $v_{t}(A_{ig})$’s.

In simulation study 1, we fit to our Bayesian hierarchical model to the
“observed” estimates $z_{t}^{(5)}(A_{ig})$ and $z_{t}^{(1)}(A_{i})$, whereas
in simulation study 2, we fit to them both our Bayesian hierarchical model and
the standard Binomial model for disaggregation. In each case, we run the MCMC
algorithms for 10,000 iterations, discarding the first 2,000 for burn-in.

###  4.3 Simulation Results

Simulation study 1. Taking the posterior means of the
$\pi^{(1)}_{t}(A_{ig})$’s as estimates of the true proportions, and denoting
them by $\hat{\pi}^{(1)}_{t}(A_{ig})$, $i=1,\ldots,4$; $g=1,2,\ldots,25$, we
summarize the performance of our model by evaluating, for each
$\pi^{(1)}_{t}(A_{ig})$, the magnitude of the errors and the empirical
coverage of both the 50% and the 95% pointwise and joint credible intervals,
respectively.

Tables 2 and 3 present results from our simulation studies, including the Mean
Squared Error (MSE) and the Mean Absolute Error (MAE). The latter are defined
as the mean squared difference and the mean absolute difference between the
$\hat{\pi}^{(1)}_{t}(A_{ig})$’s and the true values. In addition, Tables 2 and
3 present the mean squared relative error (MSRE) and the mean absolute
relative error (MARE) defined respectively as:

$\displaystyle MSRE$ $\displaystyle=$
$\displaystyle\frac{1}{100}\sum_{g=1}^{4}\sum_{i=1}^{25}\frac{(\hat{\pi}^{(1)}_{t}(A_{ig})-\pi^{(1)}_{t}(A_{ig}))^{2}}{\pi^{(1)}_{t}(A_{ig})},$
$\displaystyle MARE$ $\displaystyle=$
$\displaystyle\frac{1}{100}\sum_{g=1}^{4}\sum_{i=1}^{25}\frac{|\hat{\pi}^{(1)}_{t}(A_{ig})-\pi^{(1)}_{t}(A_{ig})|}{\pi^{(1)}_{t}(A_{ig})}.$
(4.3)

The 50% and 95% pointwise credible intervals, computed by taking the
appropriate percentiles of the posterior samples for each
$\pi^{(1)}_{t}(A_{ig})$, contain the true values between 46.1% and 53.4% of
the time, and between 90.7% and 94.8% of the time, respectively. Similarly,
the 50% and the 95% joint credible intervals, constructed using the method of
Sørbye and Rue (2011), yield nearly nominal coverage.

We observe that in both cases, the credible intervals corresponding to the
middle of the time-series ($t=3,4,5,6,7$) have the highest coverage
probabilities.

The low values for the squared and absolute errors indicate successful
recovery of the true $\pi^{(1)}_{t}(A_{ig})$’s. Figure 2 presents scatterplots
of the true proportions against the $\hat{\pi}^{(1)}_{t}(A_{ig})$’s: all plots
illustrate our model’s ability to disaggregate survey-based estimates of areal
proportions regardless of the data generating mechanism.

Table 2:  Simulation study 1. Results corresponding to 30 simulated datasets
generated under the first two of the four settings described in Table 1. For
each time $t$, $t=1,\ldots,10$, the table reports: (i) the average empirical
coverage of the 95% pointwise credible interval for $\pi^{(1)}_{t}(A_{ig})$,
$i=1,\ldots,4$, $g=1,\ldots,25$ averaged across the 100 subregions $A_{ig}$;
(ii) the average empirical coverage of the 50% pointwise credible interval for
$\pi^{(1)}_{t}(A_{ig})$; (iii) the average empirical coverage of the 95%
simultaneous credible interval for
$\boldsymbol{\pi}^{(1)}(\mathcal{S})=\left\\{\pi^{(1)}_{t}(A_{ig}):i=1,\ldots,4;g=1,\ldots,25\right\\}$
averaged across the 30 simulated datasets; (iv) the average empirical coverage
of the 50% simultaneous credible interval for
$\boldsymbol{\pi}^{(1)}(\mathcal{S})$; (v) the mean squared error (MSE); (vi)
the mean absolute error (MAE); (vii) the mean squared relative error (MSRE);
and (viii) the mean absolute relative error (MARE) as defined in (4.3).

| Average | Average | Average | Average |  |  |  |
---|---|---|---|---|---|---|---|---
| Coverage | Coverage | Coverage | Coverage | MSE | MAE | MSRE | MARE
$t$ | 95% CI pointwise | 50% CI pointwise | 95% CI joint | 50% CI joint | $\times 10^{3}$ | $\times 10^{2}$ | $\times 10^{2}$ | $\times 10^{2}$
$1$ | $91.9\%$ | 47.7% | $90.0\%$ | 46.7% | $10.3$ | $8.3$ | $4.2$ | $25.2$
$2$ | $92.9\%$ | 48.0% | $90.0\%$ | 46.7% | $8.5$ | $7.3$ | $2.6$ | $22.0$
$3$ | $92.5\%$ | 48.9% | $93.3\%$ | 50.0% | $7.3$ | $6.8$ | $2.1$ | $18.8$
$4$ | $92.6\%$ | 49.4% | $96.7\%$ | 50.0% | $6.7$ | $6.5$ | $1.9$ | $16.9$
$5$ | $93.5\%$ | 49.6% | $96.7\%$ | 50.0% | $6.4$ | $6.3$ | $1.7$ | $16.1$
$6$ | $93.5\%$ | 50.3% | $96.7\%$ | 50.0% | $6.2$ | $6.2$ | $1.7$ | $16.1$
$7$ | $92.4\%$ | 49.8% | $93.3\%$ | 50.0% | $6.6$ | $6.4$ | $1.8$ | $17.0$
$8$ | $91.8\%$ | 50.3% | $93.3\%$ | 50.0% | $7.0$ | $6.5$ | $2.0$ | $19.3$
$9$ | $91.6\%$ | 48.9% | $90.0\%$ | 46.7% | $8.1$ | $7.2$ | $2.4$ | $21.4$
$10$ | $90.7\%$ | 47.4% | $90.0\%$ | 46.7% | $11.0$ | $8.5$ | $3.8$ | $26.1$
(a) Simulation study 1, setting 1

| Average | Average | Average | Average |  |  |  |
---|---|---|---|---|---|---|---|---
| Coverage | Coverage | Coverage | Coverage | MSE | MAE | MSRE | MARE
$t$ | 95% CI pointwise | 50% CI pointwise | 95% CI joint | 50% CI joint | $\times 10^{3}$ | $\times 10^{2}$ | $\times 10^{2}$ | $\times 10^{2}$
$1$ | $91.3\%$ | 47.7% | $90.0\%$ | 46.7% | $9.7$ | $8.1$ | $4.0$ | $24.8$
$2$ | $92.0\%$ | 49.7% | $93.3\%$ | 46.7% | $6.3$ | $6.4$ | $2.3$ | $22.3$
$3$ | $92.3\%$ | 51.9% | $93.3\%$ | 50.0% | $5.4$ | $5.8$ | $1.9$ | $21.4$
$4$ | $94.8\%$ | 52.9% | $93.3\%$ | 50.0% | $4.8$ | $5.3$ | $1.5$ | $19.0$
$5$ | $94.4\%$ | 53.4% | $96.7\%$ | 50.0% | $5.0$ | $5.4$ | $1.5$ | $17.4$
$6$ | $94.0\%$ | 53.0% | $96.7\%$ | 50.0% | $5.0$ | $5.4$ | $1.4$ | $16.2$
$7$ | $93.9\%$ | 52.8% | $93.3\%$ | 50.0% | $5.0$ | $5.5$ | $1.6$ | $17.9$
$8$ | $94.0\%$ | 51.0% | $93.3\%$ | 50.0% | $5.3$ | $5.8$ | $1.8$ | $19.4$
$9$ | $92.8\%$ | 50.4% | $93.3\%$ | 50.0% | $6.1$ | $6.3$ | $2.2$ | $20.6$
$10$ | $91.2\%$ | 47.2% | $90.0\%$ | 46.7% | $9.4$ | $8.0$ | $4.0$ | $25.6$
(b) Simulation study 1, setting 2

Table 3:  Simulation study 1. Results corresponding to 30 simulated datasets
generated under the last two of the four settings described in Table 1. For
each time $t$, $t=1,\ldots,10$, the table reports: (i) the average empirical
coverage of the 95% pointwise credible interval for $\pi^{(1)}_{t}(A_{ig})$,
$i=1,\ldots,4$, $g=1,\ldots,25$ averaged across the 100 subregions $A_{ig}$;
(ii) the average empirical coverage of the 50% pointwise credible interval for
$\pi^{(1)}_{t}(A_{ig})$; (iii) the average empirical coverage of the 95%
simultaneous credible interval for
$\boldsymbol{\pi}^{(1)}(\mathcal{S})=\left\\{\pi^{(1)}_{t}(A_{ig}):i=1,\ldots,4;g=1,\ldots,25\right\\}$
averaged across the 30 simulated datasets; (iv) the average empirical coverage
of the 50% simultaneous credible interval for
$\boldsymbol{\pi}^{(1)}(\mathcal{S})$; (v) the mean squared error (MSE); (vi)
the mean absolute error (MAE); (vii) the mean squared relative error (MSRE);
and (viii) the mean absolute relative error (MARE) as defined in (4.3).

| Average | Average | Average | Average |  |  |  |
---|---|---|---|---|---|---|---|---
| Coverage | Coverage | Coverage | Coverage | MSE | MAE | MSRE | MARE
$t$ | 95% CI pointwise | 50% CI pointwise | 95% CI joint | 50% CI joint | $\times 10^{3}$ | $\times 10^{2}$ | $\times 10^{2}$ | $\times 10^{2}$
$1$ | $91.9\%$ | 46.8% | $90.0\%$ | 46.7% | $9.2$ | $7.8$ | $5.8$ | $36.1$
$2$ | $92.5\%$ | 48.1% | $90.0\%$ | 46.7% | $7.1$ | $6.7$ | $3.0$ | $31.0$
$3$ | $91.7\%$ | 48.9% | $90.0\%$ | 46.7% | $6.4$ | $6.4$ | $2.4$ | $25.8$
$4$ | $92.9\%$ | 51.1% | $93.3\%$ | 53.3% | $6.2$ | $6.1$ | $1.9$ | $20.6$
$5$ | $92.9\%$ | 50.9% | $93.3\%$ | 53.3% | $6.0$ | $6.1$ | $1.6$ | $17.9$
$6$ | $93.0\%$ | 50.6% | $93.3\%$ | 50.0% | $6.0$ | $6.1$ | $1.6$ | $14.3$
$7$ | $93.0\%$ | 50.8% | $93.3\%$ | 50.0% | $6.2$ | $6.2$ | $1.5$ | $15.2$
$8$ | $92.6\%$ | 48.3% | $93.3\%$ | 46.7% | $6.7$ | $6.5$ | $1.6$ | $15.9$
$9$ | $91.2\%$ | 48.6% | $90.0\%$ | 46.7% | $7.2$ | $6.8$ | $1.6$ | $16.1$
$10$ | $91.3\%$ | 47.9% | $90.0\%$ | 46.7% | $8.9$ | $7.7$ | $1.8$ | $17.6$
(c) Simulation study 1, setting 3

| Average | Average | Average | Average |  |  |  |
---|---|---|---|---|---|---|---|---
| Coverage | Coverage | Coverage | Coverage | MSE | MAE | MSRE | MARE
$t$ | 95% CI pointwise | 50% CI pointwise | 95% CI joint | 50% CI joint | $\times 10^{3}$ | $\times 10^{2}$ | $\times 10^{2}$ | $\times 10^{2}$
$1$ | $91.8\%$ | 45.3% | $90.0\%$ | 46.7% | $9.1$ | $7.7$ | $5.2$ | $34.8$
$2$ | $91.3\%$ | 49.3% | $90.0\%$ | 46.7% | $6.5$ | $6.4$ | $3.3$ | $30.9$
$3$ | $91.1\%$ | 49.8% | $90.0\%$ | 50.0% | $5.7$ | $6.0$ | $2.5$ | $25.1$
$4$ | $92.2\%$ | 51.8% | $93.3\%$ | 53.3% | $5.3$ | $5.6$ | $1.8$ | $20.1$
$5$ | $92.6\%$ | 52.0% | $93.3\%$ | 53.3% | $5.1$ | $5.5$ | $1.5$ | $16.3$
$6$ | $93.0\%$ | 51.9% | $93.3\%$ | 53.3% | $5.2$ | $5.5$ | $1.5$ | $14.7$
$7$ | $92.7\%$ | 51.7% | $93.3\%$ | 50.0% | $5.3$ | $5.6$ | $1.5$ | $15.0$
$8$ | $92.2\%$ | 49.8% | $93.3\%$ | 46.7% | $5.9$ | $6.1$ | $1.6$ | $15.3$
$9$ | $92.0\%$ | 48.7% | $93.3\%$ | 46.7% | $6.7$ | $6.5$ | $1.7$ | $16.8$
$10$ | $91.9\%$ | 46.1% | $90.0\%$ | 46.7% | $9.2$ | $7.8$ | $1.8$ | $18.1$
(d) Simulation study 1, setting 4

(e) Simulation study 1, setting 1; $t$ = 2 (f) Simulation study 1, setting 1;
$t$ = 5

(g) Simulation study 1, setting 1; $t$ = 9 (h) Simulation study 1, setting 2;
$t$ = 2

(i) Simulation study 1, setting 2; $t$ = 5 (j) Simulation study 1, setting 2;
$t$ = 9

(k) Simulation study 1, setting 3; $t$ = 2 (l) Simulation study 1, setting 3;
$t$ = 5

(m) Simulation study 1, setting 3; $t$ = 9 (n) Simulation study 1, setting 4;
$t$ = 2

(o) Simulation study 1, setting 4; $t$ = 5 (p) Simulation study 1, setting 4;
$t$ = 9

Figure 2:  Simulation study 1. Scatterplots of the true
$\pi^{(1)}_{t}(A_{ig})$’s against their corresponding estimates,
$\hat{\pi}^{(1)}_{t}(A_{ig})$, at selected times, $t=2,5,9$. The simulated
values and their estimates refer to all the 30 simulated datasets generated
under one of the four different simulation settings described in Table 1.

Simulation study 2. Here we compare the performance of our proposed model to
that of the standard Binomial model for disaggregation. Taking again the
posterior means of the $\pi^{(1)}_{t}(A_{ig})$’s as our estimates,
$\hat{\pi}^{(1)}_{t}(A_{ig})$’s, $g=1,\ldots,25;i=1,\ldots,4;t=1,\ldots,10$,
Table 4 presents, for each model, the average mean squared error and the
average mean absolute error, averaged over areal units, time points, and
simulated datasets. Conversely, Table 5 provides the empirical coverage
probabilities of 50% and 95% pointwise and joint credible intervals.

We inspect the difference in inference provided by the two models as the
design effect varies. When $d=2$, there is little difference between the
standard Binomial model and our model. However, when $d=4,6$ or $8$, the
standard Binomial model has an inferior performance with respect to each of
the metrics considered, suggesting that by ignoring the design effect, the
standard Binomial model places too much certainty in the pseudo-survey-
estimates that we have generated. On the other hand, by correctly propagating
the uncertainty of the pseudo-survey-estimates through the use of the
effective sample size and the effective number of cases, our proposed model
achieves lower mean squared and lower mean absolute error, as well as near
nominal coverage probability.

Table 4: Simulation study 2. Average probability that a 95%, respectively, a 50% pointwise, respectively, joint credible interval covers the true value. Averages are taken over areal units, time points, datasets, and simulation settings for the pointwise credible intervals, whereas they are taken over time points, datasets and simulation settings for the joint credible intervals. | Coverage | Coverage | Coverage | Coverage | Coverage | Coverage | Coverage | Coverage
---|---|---|---|---|---|---|---|---
| 95% CI | 95% CI | 50% CI | 50% CI | 95% CI | 95% CI | 50% CI | 50% CI
| pointwise - | pointwise - | pointwise - | pointwise - | joint - | joint - | joint - | joint -
| Proposed | Standard | Proposed | Standard | Proposed | Standard | Proposed | Standard
$d$ | model | Binomial | model | Binomial | model | Binomial | model | Binomial
2 | 93.2% | 89.7% | 53.3% | 46.7% | 93.1% | 93.1% | 53.8% | 53.8%
4 | 95.4% | 87.4% | 53.7% | 44.7% | 92.1% | 88.7% | 51.4% | 45.2%
6 | 93.2% | 83.6% | 52.8% | 42.6% | 91.2% | 83.6% | 50.1% | 41.0%
8 | 93.9% | 80.4% | 52.3% | 40.1% | 93.1% | 79.3% | 48.8% | 39.8%
Table 5: Simulation study 2. Average mean squared error (MSE) $\times 10^{3}$ and average mean absolute error (MAE) $\times 10^{2}$ computed by taking, respectively, the squared difference between the estimated $\hat{\pi}^{(1)}_{t}(A_{ig})$ and the true value $\pi^{(1)}_{t}(A_{ig})$, and the absolute value of the difference between the estimated $\hat{\pi}^{(1)}_{t}(A_{ig})$ and the true value $\pi^{(1)}_{t}(A_{ig})$ for $i=1,\ldots,4;g=1,\ldots,25;t=1,\ldots,10$. For each modeling framework, averages are taken over areal units, time points, datasets, and simulation settings. | MSE $\times 10^{3}$ | MSE $\times 10^{3}$ | MAE $\times 10^{2}$ | MAE $\times 10^{2}$
---|---|---|---|---
$d$ | Proposed model | Standard Binomial | Proposed model | Standard Binomial
2 | 5.7 | 5.9 | 5.2 | 5.2
4 | 6.0 | 12.3 | 5.7 | 7.3
6 | 8.6 | 14.7 | 6.2 | 9.0
8 | 9.0 | 21.9 | 7.0 | 11.7

##  5 Data Analysis

###  5.1 Families in poverty

We apply the model in Sections 3.1 to 3.5 to the ACS estimates of the
proportion of families in Michigan living in poverty from 2006 to 2016, with
the goal of obtaining annual estimates at the census tract level.

We present results from our model in a variety of ways, including a comparison
of the mean and variance of our model-based estimates to those provided in the
ACS dataset. We also compare the out-of-sample predictive performance of our
model to that of the three competing models described in Section 3.6, which we
also fit to the 2006-2016 ACS data. However, our primary focus is on the
disaggregated estimates for selected neighborhoods in Detroit. Here, we chose
a set of census tracts in Midtown, a mixed-use area in Detroit located north
of downtown and comprising several business districts, Wayne State University,
and some residential neighborhoods. Some of the census tracts in Midtown have
very high poverty, while others host various sporting arenas and other
downtown attractions, and thus exhibit considerably lower poverty rate. Some
of the high-poverty tracts have been subject to gentrification and development
in recent years (Moehlman and Robins-Somerville, 2016; Aguilar, 2015). Due to
these spatial inhomogeneities and temporal changes, we believe that these
tracts illustrate the need for fine scale spatio-temporal estimates in order
to properly characterize neighborhood surroundings.

####  5.1.1 Comparison of 5-year model-based estimates as estimated by our
model to ACS 5-year estimates

(a) Model-based vs. ACS estimates (b) Posterior SD vs. ACS SE’s

Symbol | Tract ID | ACS Estimate | ACS SE
---|---|---|---
$\blacktriangle$ | 26163517200 | 0.00 | 0.19
$\bullet$ | 26161400100 | 0.57 | 0.20
$\blacktriangle$ | 26163550800 | 0.00 | 0 .04
$\blacktriangle$ | 26073000700 | 0.63 | 0.11
$\blacktriangle$ | 26057000400 | 0.00 | 0.61
$\blacktriangle$ | 26037011200 | 0.00 | 0.17
$\blacktriangle$ | 26163500400 | 0.64 | 0.10
$\blacktriangle$ | 26163512900 | 0.72 | 0.11
(c) ACS estimates and SE’s for tracts highlighted in (a)

Figure 3:  (a) Model-based estimates of 5-year average proportion of families
in poverty in Michigan at census tract level for years 2009-2013 against
corresponding ACS estimates. Census tracts deviating greatly from the identity
line are denoted by triangles. (b) Posterior standard deviation for the 5-year
average proportion of families in poverty in Michigan at census tract level as
yielded by our model against the ACS standard error of the corresponding
estimates. (c) Tabulation of Tract ID’s, ACS estimates, and ACS standard
errors for census tracts for which our model-based estimate deviates greatly
from the ACS estimate.

Our model is intended for spatial and temporal disaggregation, but we can also
generate 5-year census tract estimates, which should resemble the
corresponding estimates from ACS. Figure 3(a) shows a scatter plot comparing
these estimates. The points tend to fall around the identity line, indicating
good agreement. Figure 3(b) compares the standard errors of the ACS 5-year
census tract estimates with the posterior standard deviations from our model.
As our model borrows information from neighboring census tracts and from the
1-year PUMA-level estimates, it yields estimates with smaller posterior
standard deviation than the ACS standard errors. Many of the points that
deviate from the identity line in Figure 3(a) correspond to census tracts with
zero-valued ACS estimates, which have large ACS standard errors compared to
the average of 0.041 over all tracts (Figure 3(c)).

An example of such a census tract is displayed in panels (a) and (b) of Figure
4 along with its neighboring tracts. Census tract 26161400100 is located in
downtown Ann Arbor and, according to the ACS estimate, has an average poverty
rate of 0.59 for the 5-year time period from 2009 to 2013. This estimate
deviates greatly from that of the neighboring tracts. In addition, it has a
design-based standard error around 5 times the average standard error for ACS
estimates of poverty in Michigan. As our model borrows information from
neighboring census tracts, the model-based estimate for this tract is pushed
towards the average of the neighboring tracts more than towards the raw ACS
estimates. We observe regression of a model-based estimate towards the average
of its neighbors in situations where the design-based standard error is quite
large.

(a) ACS estimates, Ann Arbor (b) Model-based estimates, Ann Arbor

(c) ACS estimates, Romulus (d) Model-based estimates, Romulus

(e) Posterior density

Figure 4:  (a) ACS estimate for the 5-year average proportion of families in
poverty in Ann Arbor census tract 26161400100 and (b) our corresponding model-
based estimate. (c) ACS estimate for the 5-year average proportion of families
in poverty in a census tract in Romulus and (d) our model-based estimate.
Here, despite the lower-poverty level in the neighboring census tracts, our
model-based estimate is not smoothed towards the poverty-level of the
neighboring tracts. (e) Posterior densities of: (i) 5-year average proportion
and (ii) 1-year proportion of families living in poverty in census tract
26163563500 according to our model, as well as the truncated normal density
function obtained using as mean and standard deviation, respectively, the
5-year ACS estimate and its corresponding standard error.

To illustrate this phenomenon, Figure 4 shows in panels (c) and (d) a census
tract in Romulus also characterized by high poverty despite being surrounded
by census tracts with lower poverty level, according to the ACS. In this case,
since the ACS standard error is much lower, our model-based estimate of
poverty still reflects the spatial heterogeneity in the ACS estimates and it
is not smoothed towards the average of the neighboring tracts.

####  5.1.2 Posterior density of true population proportions

ACS estimates are provided with margins of error based on standard errors
derived from an iterative estimation procedure (see U.S. Census Bureau (2014)
for details) multiplied by standard normal quantiles. These may in turn be
utilized to construct confidence intervals for the estimates by adding and
subtracting their margins of error, with users being cautioned to use “logical
boundaries when creating confidence bounds from the margins of error” (U.S.
Census Bureau, 2008) (i.e. zero and one for proportions). This implies a
truncated normal distribution centered at the ACS estimate with variance
depending on the standard error.

Through our Bayesian modeling framework, we obtain the posterior distribution
of the true proportions at any spatial and temporal scale without imposing
assumptions of symmetry or truncation. For example, Figure 3(e) plots the
posterior density of the average proportion of families living in poverty in
census tract 26163563500 located in Midtown Detroit for the 5-year period
2009-2013. For this census tract, the confidence interval for the 5-year ACS
estimate for 2009-2013 would be subject to truncation at zero. Figure 3(e)
shows the truncated normal density with mean and standard deviation equal,
respectively, to the ACS estimate and its standard error. To facilitate direct
comparison to the ACS estimates, Figure 3(e) also presents the posterior
density of the 5-year average proportion of families living in poverty as
provided by our model, as well as the posterior density for the 1-year
proportions for years 2009, 2010, 2011, 2012 and 2013. As the figure shows,
thanks to the borrowing of information from neighboring census tracts, the
posterior density of the 5-year average proportion is characterized by smaller
uncertainty than the truncated normal density centered at the ACS estimate.

####  5.1.3 Disaggregated estimates of poverty for Detroit

Disaggregating the ACS estimates allows us to examine yearly trends in poverty
for individual census tracts, as well as for combinations of census tracts
which do not form a PUMA or are not part of a highly populated county. For
both of these cases we cannot assess temporal trends using the ACS estimates
alone. Figure 5 presents various maps of the disaggregated estimates of the
proportion of families in poverty in Michigan from our model. Panel (a)
displays census tract estimates for all of Michigan for year 2010 while panel
(b) presents the same results for Wayne County, which contains areas of
extreme poverty, as well as some of the wealthiest neighborhoods in Michigan.

Panels (c)-(k) of Figure 5 present yearly results over time for a subset of
census tracts in Wayne County located in Midtown Detroit. This set of census
tracts was selected because the poverty rates exhibit spatial heterogeneity,
with certain pairs of neighboring census tracts differing by over 20%, so a
single poverty estimate for this area would not properly characterize the
neighborhood conditions of its residents. Recent changes in these census
tracts have been well-documented (Moehlman and Robins-Somerville, 2016),
particularly in the Cass Corridor, an area of downtown Detroit that has faced
high crime and poverty, but has recently experienced sudden gentrification
(Aguilar, 2015). The 5-year ACS estimates at the census tract level may not
properly characterize yearly changes in these census tracts.

(a) Model-based estimates: Michigan 2010 (b) Model-based estimates: Wayne
County 2010

(c) 2007 (d) 2008 (e) 2009 (f) 2010 (g) 2011 (h) 2012 (i) 2013 (j) 2014 (k)
2015

Figure 5: Disaggregated estimates of the proportion of families in poverty in
Midtown Detroit.

Figure 6(a) shows the changes over time of poverty rates for a set of census
tracts in the Midtown area of Detroit. Consistent with national trends, Figure
6(b) indicates that, on average, the area experienced an increase in poverty
following the 2008 financial crisis in the US, with an eventual improvement in
later years. Figure 5, panels (c)-(k), highlights a census tract of particular
interest, indicated in green in Figure 6(a), which did not experience a
decrease in poverty until 2014. Year to year, it has had among the highest
poverty rates in Detroit. However, recent developments such as the
groundbreaking of Little Caesars Arena in 2014 and an influx of newly built
restaurants and bars, might have contributed to the drop in poverty in that
census tract from 2014 to 2015 (Moehlman and Robins-Somerville, 2016).

(a) Spaghetti plot of poverty for Midtown census tracts. (b) Mean poverty for
Midtown census tracts

Figure 6:  (a) Spaghetti plot displaying the estimated proportions of families
in poverty over time with highlighted the census tract shown in Figure 5 and
(b) the estimated average poverty rate across census tracts in Midtown Detroit
between 2006 and 2016.

###  5.2 Out-of-Sample Prediction

To assess our model’s out-of-sample predictive performance, we consider the
county-level proportion of families in poverty during the 3-year time periods
2010-2012 and 2011-2013. We generate 3-year county-level predictions as a
weighted average of the disaggregated estimates within each county’s census
tracts for the appropriate 3-year periods, with weights proportional to the
number of families living in each tract. Then, those yearly county estimates
are averaged over the 3-year time periods. Our “true values” are the 3-year
ACS estimates, which we did not use for model fitting. This allows us to
assess our model’s ability to predict over time periods and areal units that
are not utilized in model fitting.

Figure 7(a) compares the estimated 3-year proportion yielded by our model with
the ACS estimates. As the figure shows, the two sets of estimates tend to be
very similar, indicating the strong predictive performance of our model. The
mean squared and mean absolute prediction errors are, respectively,
4.16$\times 10^{-5}$ and 4.83$\times 10^{-3}$, whereas the mean squared
relative prediction error is 3.13$\times 10^{-4}$ and the mean absolute
relative prediction error is 4.18$\times 10^{-2}$. These values demonstrate
reduced predictive error in our modeling framework.

###  5.3 Comparison to other models

We also generate out-of-sample predictions of the proportions of families in
poverty at the 3-year county resolution for the three competing models
presented in Section 3.6. Details on how these predictions are derived are
provided in Section 3 of the Supplementary Material (Benedetti, Berrocal and
Little, 2021). We evaluate the quality of these out-of-sample predictions by
validating them against the ACS estimates: Figure 7(b), (c) and (d) show
scatter plots of the predicted proportion of families in poverty as yielded by
each of the three competing models against the ACS estimate. The standard
Binomial model and the BWH Poisson model produce estimates that are fairly in
line with the ACS values, while there is a larger discrepancy between the
estimates yielded by the BWH Gaussian Delta Method model and the ACS 3-year
estimates. Numerically, we compare the predictive performance of our proposed
model to that of the other 3 models in terms of average predictive bias, mean
squared predictive error, mean absolute predictive error, and coverage of the
50% and 95% prediction intervals. These statistics are reported in Table 6,
with the first three summary statistics all functions of the difference
between the predicted proportions and the ACS 3-year estimates.

As the table shows, our model performs almost equivalently to the standard
Binomial model in terms of predictive accuracy with our model yielding a
coverage slightly closer to the nominal level than the standard Binomial
model. This occurs for both the 50% and the 95% prediction intervals. However,
the difference is minimal: we attribute this to the large sample size and
careful sampling design of the ACS, which limits the impact of the design
effect on the model’s performance.

Moving onto the BWH Poisson space-time model, we can see that even though this
model exhibits accurate predictions, our model is slightly more accurate. The
main differences between the two models are with respect to the posterior
predictive standard deviations: the BWH Poisson model has much smaller
posterior predictive standard deviations and thus much lower coverage
probabilities than our model.
Finally, the BWH Gaussian Delta Method model offers a poorer predictive
performance than our model with respect to all metrics. We acknowledge that
the models that we have attributed to Bradley, Wikle, and Holan are not
necessarily the approaches that the authors would have taken to model the ACS
spatio-temporal estimates of proportions. Rather, they constitute our best
effort to adapt the methods presented in Bradley, Wikle and Holan (2015) and
Bradley, Wikle and Holan (2016) to model our data. While we had to modify both
models to accommodate estimates of proportions, we took care to do so in a way
that would not needlessly favor our model.

(a) Our proposed model (b) Standard binomial model (c) BWH Poisson space-time
(d) BWH Gaussian Delta Method

Figure 7: Comparison to other models. Predicted proportion of families in poverty vs. ACS 3-year estimates of the proportion of families in poverty in Michigan counties for the period 2010-2012 and 2011-2013 as yielded by: (a) our proposed model, (b) the Standard Binomial model, (c) the BWH Poisson space-time model and (d) the BWH Gaussian Delta method model.  Table 6: Comparison to other models. Bias, Mean Squared Predictive Error (MSPE), Mean Absolute Predictive Error (MAPE) of the out-of-sample predictions, as well as empirical coverage of the 50% and the 95% prediction intervals (PI) for our proposed model and the three competing models. | Bias | MSPE | MAPE | Coverage | Coverage
---|---|---|---|---|---
Model | $\times 10^{2}$ | $\times 10^{5}$ | $\times 10^{3}$ | 50% PI | 95% PI
Our proposed model | 0.03 | 4.16 | 4.83 | 52.3% | 93.0%
Standard Binomial |  $-$0.03 | 3.69 | 4.39 | 46.9% | 91.4%
BWH Poisson space-time |  $-$0.70 | 13.91 | 9.52 | 9.4% | 23.4%
BWH Gaussian delta method |  $-$1.70 | 52.49 | 18.11 | 26.6% | 52.3%

##  6 Discussion

This paper proposes a spatio-temporal Bayesian hierarchical model to
disaggregate estimates of proportions over areal units derived from sampling
surveys while accounting for the survey design. Previous to our work, Bradley,
Wikle and Holan (2016) formulated a stochastic model for ACS estimates
distributed according to a Poisson distribution. The model explicitly
accounted for the survey design, as it specified a lognormal distribution for
the ACS design-based variance; however it focused only on addressing the
change of support problem in a spatial setting for count variables. Other work
by Bradley, Wikle and Holan (2015) considered the space-time setting, but it
did not incorporate design effects as it postulated a Gaussian likelihood with
the ACS design-based variance taken as known and set equal to the variance of
the normal distribution.

The main motivation for the development of our modeling framework is the
ability to generate data on socio-economic indicators at fine spatial and
temporal resolution, thus responding to the needs of health researchers
investigating the effect of social determinants of health on health outcomes.
We have demonstrated the utility of our Bayesian hierarchical modeling
framework by applying it to the ACS estimates of families in poverty. This
application highlighted several advantages of our model, among which the fact
that it generates annual estimates at census tract spatial resolution. In
addition, due to the borrowing of information from neighboring units and from
ACS estimates at different spatial and temporal resolutions, these estimates
are characterized by smaller uncertainty. We use our disaggregated estimates
to examine trends over time of poverty in Michigan focusing on Detroit, for
which we could highlight yearly changes at small spatial scale. These changes
could not be detected easily using the 5-year ACS census tract estimates.

We recognize that a standard Binomial disaggregation model that did not
explicitly account for the design effect, applied to the same data, yields a
very comparable, or slightly better, predictive performance than our model
when evaluated in terms of Mean Squared Predictive Error, Mean Absolute
Predictive Error, and Empirical Coverage of the 50% and 95% pointwise
prediction intervals. However, our simulation study 2 also indicated that
while such a performance by the standard Binomial model is expected for small
design effect, the aforementioned model is subject to a worsening in
predictive performance as the design effect size increases. In these
situations, our model is preferable. While it is true that the ACS design
effect for the estimates of proportion of families in poverty in Michigan
during the period considered – 2006-2016 – is estimated to be around 2.5 on
average, our model has not been developed only for handling estimates
resulting from the ACS. Rather, our model has a wider applicability and it has
been formulated to disaggregate spatially and temporally any set of survey-
based multi-year estimates of proportions. Other surveys typically used in
epidemiological studies are characterized by larger design effects than the
ACS. For example, the national Behavior Risk Factor Surveillance System
(BRFSS) in 2013 had an average design effect of 4.45, with state BRFSS surveys
having design effects ranging from 1.47 to 5.16 (Iachan et al., 2016). For
estimates provided by these surveys, our model is expected to yield better
results than the standard Binomial model.

Our model is not the only one using the concept of design effect to yield
small-scale spatio-temporal estimates: Li et al. (2019) applied the model of
Mercer et al. (2014) to smooth the spatial distribution of 1-year estimates of
the under-5 mortality rate over 35 countries in Africa, producing subnational
estimates at the 1-year resolution. Although Li et al. (2019) are also
concerned with generating small scale spatio-temporal estimates, their work
did not address the spatio-temporal change of support problem in the same way
we do here. Specifically, Li et al. (2019) introduce a spatio-temporal process
that is discrete both in space and time, whereas our model employs a point-
referenced spatio-temporal process, discrete in time but continuous in space.
Additionally, our model accounts for the clustering units of the survey (e.g.
counties in our application) by including random effects specified at the
county-level spatial resolution.

To deal with the large dimensionality of the data, we approximate the latent
spatio-temporal process driving the true population proportions via a basis
function expansion. Multiple choices are available to alleviate the
computational burden associated with fitting a spatial statistical model to
large spatial data, as reviewed by Heaton et al. (2019). Here, acknowledging
the nested and multi-resolution geography of the ACS data, also noted by
Savitsky (2016), we elect to choose the Multi-Resolution Approximation (MRA)
of Katzfuss (2017), extending it to the space-time setting, an additional
contribution of our paper. However, other basis functions could be employed,
namely, wavelets, radial basis functions, and Moran’s basis functions as in
Bradley, Wikle and Holan (2016).

In our model, partly for computational considerations, we use a probit link to
relate the true areal-level proportions to the underlying Gaussian spatio-
temporal process, and we employ the data augmentation algorithm of Albert and
Chib (1993) for posterior computation. We believe that it is possible to
devise an MCMC algorithm based on the skew-normal posterior results for probit
regressions derived by Durante (2019). Additionally, we remark that one could
replace the probit link with a logit link. In this case, we encourage readers
to employ a Pólya-Gamma augmentation scheme (Polson, Scott and Windle, 2013)
for greater computational efficiency.

Much of our predictive performance evaluation is based on the empirical
probability that credible and/or prediction intervals cover the true value,
which inherently conflates a frequentist property (empirical coverage
probability) with Bayesian modeling frameworks. This type of assessment is in
line with the notion of calibrated Bayes (Little, 2006) and recommended in a
predictive context (Dawid, 1982); moreover, it is the authors’ experience that
coverage probabilities are frequently used in assessing Bayesian models,
particularly in a spatial context (see Entezari, Brown and Rosenthal (2019);
Gilani, Berrocal and Batterman (2019); Berrocal, Gelfand and Holland (2010) as
example), where prediction is the main goal.

We note a potential abuse of terminology in calling “out-of-sample validation”
the comparison of the 3-year county-level proportions yielded by our model
with the corresponding ACS estimates. Even though the 3-year county level ACS
estimates were not used in fitting the model, the microdata that is leveraged
to derive such ACS estimates is also employed to calculate the 1-year and
5-year estimates of proportions to which our model was fit. In adopting this
terminology, we follow previous examples in the literature on this topic, see
Bradley, Wikle and Holan (2015), where this type of assessment was performed
and this nomenclature was used.

Finally, a characteristic of our model is the assumption of conditional
independence between the 1-year ACS PUMA-level estimates and the 5-year census
tract estimates, conditional on the true areal proportions. Since both sets of
estimates are derived using the same microdata, it is possible that the
assumption of conditional independence is not realistic. Not having access to
the actual microdata, we have no means to determine whether this assumption is
violated. Future work could be devoted to relax the assumption of conditional
independence.

Supplementary Information In the Supplementary Material, we derive statistical
properties for the ST-MRA method, provide details on how predictions were
derived, and present results of the exploratory data analysis described in
Section 2.2. Specifically: Section 1 shows that the ST-MRA expression
presented in Section 3.3 provides an approximation to a Gaussian spatio-
temporal process with a separable covariance function, with an AR(1) structure
in time and a dependence structure in space encoded by a Matérn covariance
function. Section 2 discusses how to derive out-of-sample predictions under
the alternative models discussed in Section 3.6, while Section 3 shows results
of the exploratory data analysis that supports our modeling choices. Finally,
Section 4 concludes the Supplementary Material presenting results for the city
of Flint.

## References

  * Abrams and Szefler (2020) [author] Abrams, E. M.E. M. and Szefler, S. J.S. J. (2020). COVID-19 and the impact of social determinants of health. The Lancet - Respiratory Medicine 8 659–661.
  * Aguilar (2015) [author] Aguilar, L.L. (2015). Detroit’s Cass Corridor makes way for new era. The Detroit News, published April 2015.
  * Albert and Chib (1993) [author] Albert, J. H.J. H. and Chib, S.S. (1993). Bayesian analysis of binary and polychotomous data. Journal of the American Statistical Association 88 669–679.
  * Banerjee, Carlin and Gelfand (2004) [author] Banerjee, S.S., Carlin, B. P.B. P. and Gelfand, A. E.A. E. (2004). Hierarchical modeling and analysis for spatial data. Chapman & Hall/CRC Boca Raton, FL.
  * Benedetti, Berrocal and Little (2021) [author] Benedetti, M. H.M. H., Berrocal, V. J.V. J. and Little, R.R. (2021). Supplement to “Accounting for survey design in Bayesian disaggregation of survey-based areal estimates of proportions: an application to the American Community Survey”.
  * Berrocal, Gelfand and Holland (2010) [author] Berrocal, V. J.V. J., Gelfand, A. E.A. E. and Holland (2010). A bivariate space-time downscaler under space and time misalignment. The Annals of Applied Statistics 4 1942–1975.
  * Bradley, Holan and Wikle (2016) [author] Bradley, J. R.J. R., Holan, S. H.S. H. and Wikle, C. K.C. K. (2016). Multivariate spatio-temporal survey fusion with application to the American Community Survey and local area unemployment statistics. Stat 5 224–233.
  * Bradley, Wikle and Holan (2015) [author] Bradley, J. R.J. R., Wikle, C. K.C. K. and Holan, S. H.S. H. (2015). Spatio-temporal change of support with application to American Community Survey multi-year period estimates. Stat 4 255–270.
  * Bradley, Wikle and Holan (2016) [author] Bradley, J. R.J. R., Wikle, C. K.C. K. and Holan, S. H.S. H. (2016). Bayesian spatial change of support for count-valued survey data with application to the American Community Survey. Journal of the American Statistical Association 111 472–487.
  * Braverman, Egerter and Williams (2011) [author] Braverman, P.P., Egerter, S.S. and Williams, D. R.D. R. (2011). The social determinants of health: coming of age. Annual Reviews of Public Health 32 381–398.
  * U.S. Census Bureau (2008) [author] U. S. Census Bureau (2008). A Compass for Understanding and Using American Community Survey Data: What General Data Users Need to Know. U.S. Government Printing Office, Washington, DC.
  * U.S. Census Bureau (2014) [author] U. S. Census Bureau (2014). American Community Survey Design and Methodology. U.S. Government Printing Office, Washington, DC.
  * Chen, Wakefield and Lumley (2014) [author] Chen, C.C., Wakefield, J.J. and Lumley, T.T. (2014). The use of sampling weights in Bayesian hierarchical models for small area estimation. Spatial and Spatio-temporal Epidemiology 11 33 – 43.
  * Dawid (1982) [author] Dawid, A. P.A. P. (1982). The well-calibrated Bayesian. Journal of the American Statistical Association 77 605–610.
  * Diggle, Tawn and Moyeed (1998) [author] Diggle, P. J.P. J., Tawn, J. A.J. A. and Moyeed, R. A.R. A. (1998). Model-based geostatistics. Journal of the Royal Statistical Society: Series C (Applied Statistics) 47 29-9-350.
  * Durante (2019) [author] Durante, D.D. (2019). Conjugate Bayes for probit regression via unified skew-normal distributions. Biometrika 106 765–779.
  * Entezari, Brown and Rosenthal (2019) [author] Entezari, R.R., Brown, P. E.P. E. and Rosenthal, J. S.J. S. (2019). Bayesian spatial analysis of hardwood tree counts via MCMC. Environmetrics 31 e2608.
  * Fay and Herriot (1979) [author] Fay, R.R. and Herriot, R.R. (1979). Estimates of income for small places: an application of James-Stein procedure to census data. Journal of the American Statistical Association 74 269–277.
  * Finley, Banerjee and Carlin (2007) [author] Finley, A.A., Banerjee, S.S. and Carlin, B. P.B. P. (2007). spBayes: an R package for univariate and multivariate hierarchical point-referenced spatial models. Journal of Statistical Software 19 1–24.
  * Gelfand, Banerjee and Gamerman (2005) [author] Gelfand, A. E.A. E., Banerjee, S.S. and Gamerman, D.D. (2005). Spatial process modelling for univariate and multivariate dynamic spatial data. Environmetrics 16 465–479.
  * Gelfand, Zhu and Carlin (2001) [author] Gelfand, A. E.A. E., Zhu, L.L. and Carlin, B. P.B. P. (2001). On the change of support problem for spatio-temporal data. Biostatistics 2 31–45.
  * Geweke (1992) [author] Geweke, J.J. (1992). Evaluating the accuracy of sampling-based approaches to calculate posterior moments. In Bayesian Statistics 4 (J. M.J. M. Bernardo, J. O.J. O. Berger, A. P.A. P. Dawid and A. F. M.A. F. M. Smith, eds.) 169–193. Clarendon Press.
  * Ghitza and Gelman (2013) [author] Ghitza, Y.Y. and Gelman, A.A. (2013). Deep interaction with MRP: election turnout and voting patterns among small electoral subgroups. American Journal of Political Science 57 762-776.
  * Gilani, Berrocal and Batterman (2019) [author] Gilani, O.O., Berrocal, V. J.V. J. and Batterman, S.S. (2019). Nonstationary spatiotemporal Bayesian data fusion for pollutants in the near-road environment. Environmetrics 30 1-19.
  * Gotway and Young (2002) [author] Gotway, C. AC. A. and Young, L. JL. J. (2002). Combining incompatible spatial data. Journal of the American Statistical Association 97 632–648.
  * Heaton et al. (2019) [author] Heaton, M. J.M. J., Datta, A.A., Finley, A. O.A. O., Furrer, R.R., Guinness, J.J., Guhaniyogi, R.R., Gerber, F.F., Gramacy, R. B.R. B., Hammerling, D.D., Katzfuss, M.M., Lindgren, F.F., Nychka, D. W.D. W., Sun, F.F. and Zammit-Mangion, A.A. (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics 24 398-425.
  * Iachan et al. (2016) [author] Iachan, RonaldoR., Pierannunzi, CarolC., Healey, KristieK., Greenlund, KurtK. and Town, MachellM. (2016). National weighting of data from the Behavioral Risk Factor Surveillance System (BRFSS). BMC Medical Research Methodology 16 1–12.
  * Katzfuss (2017) [author] Katzfuss, M.M. (2017). A Multi-Resolution Approximation for massive spatial datasets. Journal of the American Statistical Association 112 201–214.
  * Kish (1965) [author] Kish, L.L. (1965). Survey Sampling. John Wiley & Sons, New York.
  * Kish (1995) [author] Kish, L.L. (1995). Methods for Design Effects. Journal of Official Statistics 11 55–77.
  * Korn and Graubard (1998) [author] Korn, E. L.E. L. and Graubard, B. I.B. I. (1998). Confidence intervals for proportions with small expected number of positive counts estimated from survey data. Survey Methodology 24 193–201.
  * Li et al. (2019) [author] Li, ZehangZ., Hsiao, YuanY., Godwin, JessicaJ., Martin, Bryan D.B. D., Wakefield, JonJ. and Clark, Samuel J.S. J. (2019). Changes in the spatial distribution of the under-five mortality rate: Small-area analysis of 122 DHS surveys in 262 subregions of 35 countries in Africa. PLOS ONE 14 1-17.
  * Little (2006) [author] Little, R. J.R. J. (2006). Calibrated Bayes: A Bayes/frequentist roadmap. The American Statistician 60 213 – 223.
  * Marmot et al. (2012) [author] Marmot, M.M., Allen, J.J., Bell, R.R., Bloomer, E.E., Goldblatt, P. on behalf of the Consortium for the Eurpopean Review of Social Determinants of HealthP. o. b. o. t. C. f. t. E. R. o. S. D. o. H. and the Health Divide (2012). WHO European review of social determinants of health and the health divide. Lancet 380 1011–1029.
  * Mercer et al. (2014) [author] Mercer, L.L., Wakefield, J.J., Chen, C.C. and Lumley, T.T. (2014). A comparison of spatial smoothing methods for small area estimation with sampling weights. Spatial Statistics 8 69–85.
  * Mitchell, Genton and Gumpertz (2005) [author] Mitchell, Matthew W.M. W., Genton, Marc G.M. G. and Gumpertz, Marcia L.M. L. (2005). Testing for separability of space?time covariances. Environmetrics 16 819-831.
  * Moehlman and Robins-Somerville (2016) [author] Moehlman, L.L. and Robins-Somerville, M.M. (2016). The new Detroit: How gentrification has changed Detroit’s economic landscape. Michigan Daily, published September 2016.
  * Pereira and Coelho (2010) [author] Pereira, L. N.L. N. and Coelho, P.P. (2010). Small area estimation of habitation transaction using time-series and cross sectional areal-level models. Journal of Applied Statistics 37 651–666.
  * Pfeffermann (2013) [author] Pfeffermann, D.D. (2013). New important developments in small area estimation. Statistical Science 28 40–68.
  * Polson, Scott and Windle (2013) [author] Polson, Nicholas G.N. G., Scott, James G.J. G. and Windle, JesseJ. (2013). Bayesian inference for logistic models Using Pólya?Gamma latent variables. Journal of the American Statistical Association 108 1339–1349.
  * Porter et al. (2014) [author] Porter, A. T.A. T., Holan, S. H.S. H., Wikle, C. K.C. K. and Cressie, N.N. (2014). Spatial Fay-Herriot models for small area estimation with functional covariates. Statistical Methods and Applications 10 27-42.
  * Pratesi and Salvati (2008) [author] Pratesi, M.M. and Salvati, N.N. (2008). Small area estimation: the EBLUP estimator based on spatially correlated random area effects. Statistical Methods and Applications 17 113-141.
  * Roberts, Gelman and Gilks (1997) [author] Roberts, G. O.G. O., Gelman, A.A. and Gilks, W. R.W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. The Annals of Applied Probability 7 110–120.
  * Rollston and Galea (2020) [author] Rollston, R.R. and Galea, S.S. (2020). COVID-19 and the social determinants of health. American Journal of Health Promotion 34 687-689.
  * Savitsky (2016) [author] Savitsky, T. D.T. D. (2016). Bayesian nonparametric multiresolution estimation for the American Community Survey. Annals of Applied Statistics 10 2157–2181.
  * Simpson et al. (2017) [author] Simpson, D.D., Rue, H.H., Riebler, A.A., Martins, T. G.T. G. and Sørbye, S. H.S. H. (2017). Penalising model component complexity: a principled, practical approach to constructing priors. Statistical Science 32 1–28.
  * Simpson et al. (2019) [author] Simpson, M.M., Holan, S. H.S. H., Wikle, C. K.C. K. and Bradley, J. R.J. R. (2019). Interpolating distributions for populations in nested geographies using public-use data with application to the American Community Survey. Preprint available at: arXiv:1802.02626.
  * Singh, Shukla and Kundu (2005) [author] Singh, B.B., Shukla, G.G. and Kundu, D.D. (2005). Spatio-temporal models in small-area estimation. Survey Methodology 31 183–195.
  * Singu et al. (2020) [author] Singu, S.S., Acharya, A.A., Challagundla, K.K. and Byareddy, S. B.S. B. (2020). Impact of social determinants of health on the emerging COVID-19 pandemic in the United States. Frontiers in Public Health 8 406\.
  * Sørbye and Rue (2011) [author] Sørbye, Sigrunn H.S. H. and Rue, HåvardH. (2011). Simultaneous Credible Bands for Latent Gaussian Models. Scandinavian Journal of Statistics 38 712-725.