The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
84,352 characters
=0pt =0pt plus .5=0pt plus .5=.3Identification and Semiparametric Estimation of Conditional Means from Aggregate Data
\allsectionsfont{\sffamily}
\maketitle
\begin{abstract}
We introduce a new method for estimating the mean of an outcome variable
within groups when researchers only observe the average of the outcome
and group indicators across a set of aggregation units, such as
geographical areas. Existing methods for this problem, also known as
ecological inference, implicitly make strong assumptions about the
aggregation process. We first formalize weaker conditions for
identification which hold conditionally on covariates. To efficiently
control for many covariates, we propose a debiased machine learning
estimator that is based on nuisance functions restricted to a partially
linear form. Our estimator admits a semiparametric sensitivity analysis
which allows researchers to evaluate the impact of violations of the key
identifying assumption. We also propose a nonparametric test for the
identifying assumption itself. Finally, we derive asymptotically valid
confidence intervals for local, unit-level estimates under additional
assumptions. Simulations and validation on real-world data where ground
truth is available demonstrate the advantages of our approach over
existing methods. Open-source software is available which implements the
proposed methods.
\end{abstract}
\textbf{\textit{Keywords}}\quad aggregate data~\textbullet~ecological
inference~\textbullet~double/debiased machine learning
\section{Introduction}\label{sec-intro}
One of the most common statistical tasks is estimating the mean of an
outcome \(Y\) within subgroups defined by a discrete variable \(X\). In
many settings, however, researchers do not jointly observe \(Y\) and
\(X\) for each unit, but only observe the average of each variable
within a grouping variable such as geography. Despite the simplicity of
the aggregation operation, the problem of estimating the conditional
mean \(\operatorname{\mathbb{E}}[Y\mid X]\) from marginal means \(\overline Y\) and \(\overline X\),
known as \emph{ecological inference}, is far from straightforward.
For example, consider estimating how different racial and income groups
(\(X\)) are exposed to pollution (\(Y\)), while only observing the
fraction of individuals in each racial or income group in a ZIP code
(\(\overline X\)) and the average pollution exposure in the ZIP code
(\(\overline Y\)) \citep{jbaily2022air}. Because joint information about
\(X\) and \(Y\) within ZIP codes is not observed, the conditional mean
\(\operatorname{\mathbb{E}}[Y\mid X]\) cannot be point-identified from the data without further
assumptions \citep{cho2008cross}. This is a common challenge in
epidemiology, where exposure and disease are measured separately
\citep{greenland1994invited}. Another well-studied instance arises in
voting rights litigation, which requires estimating the voting behavior
within racial subgroups, while only observing precinct-level election
returns and Census statistic on race \citep{greiner2006ecological}.
Finally, statistical agencies generally report demographic, economic,
and public health data only in aggregate form, due to data collection or
privacy constraints. A long line of work has therefore tackled the
ecological inference problem, beginning with
\citet{robinson1950ecological} and Goodman
\citetext{\citeyear{goodman1953ecological}; \citeyear{goodman1959some}}.
This paper develops a new understanding of and methodology for the
ecological inference problem. We propose a semiparametrically efficient
estimator that allows researchers to minimize bias by controlling for
many covariates. Our approach formalizes previously implicit
assumptions, and improves on the issues of confounding and computational
efficiency in past work.
Existing methods largely fall into two categories. One, applicable only
when \(Y\) is bounded, focuses on partial identification of
\(\operatorname{\mathbb{E}}[Y\mid X]\) by deriving bounds on the estimand under various
assumptions
\citep{duncan1953alternative, cross2002regressions, fan2016estimation, manski2018credible, jiang2020ecological}.
In practice these bounds are often too wide; intervals for different
levels of \(X\) almost always overlap in practice.\footnote{ Some
authors
\citep[e.g.,][]{judge20047, muzellec2017tsallis, bontemps2025functional}
have proposed selecting a point from the partial identification region
according to an ad-hoc criterion, such as entropy minimization or
divergences based on optimal transport. While providing a single
estimate, these proposals lack statistical or substantive
justification, and as such, it is not possible to quantify their bias
or uncertainty.}
The other approach, following Goodman, aims for point identification
under often-unstated assumptions, generally relying on parametric
regression models. The seminal work of \citet{king1997solution} combined
the regression framework with the partial identification bounds in a
Bayesian varying coefficient model, spurring numerous extensions
\citep{rosen2001bayesian, wakefield2004ecological, james2009r}, though
the approach was somewhat controversial \citep{freedman1998solution}.
Our proposed method addresses two general challenges in this existing
literature. First, the necessary identifying assumptions are rarely
stated explicitly and are often extremely strong
\citep{freedman1998solution}; in practice, users of ecological inference
methods rarely evaluate the plausibility of their
assumptions.\footnote{ An exception to the pattern is
\citet{imai2008bayesian}, who presented two possible parametric
identifying assumptions. Nevertheless, awareness among practitioners
of the necessary identifying assumptions remains low.} Second, most
methods rely on strong parametric assumptions, are computationally
intensive, and lack any inferential guarantees.\footnote{ The
computational limitations may have discouraged practitioners from
controlling for essential covariates. Existing methods have
complicated likelihood functions that require computationally
intensive inference methods such as Markov chain Monte Carlo (MCMC)
algorithms, and exhibit significant slowdowns as the number of
covariates increases. \citet{km2025review} discuss these computational
challenges in more detail.}
To address the identification ambiguity, we formalize the ecological
inference problem in Section~\ref{sec-ident} with an explicit model for
the aggregation process, which allows us to relate individual-level and
aggregate-level data. Under this framework, we state two main
identifying assumptions, one at the individual level and one at the
aggregate level, and prove they are sufficient for identification of
\(\operatorname{\mathbb{E}}[Y\mid X]\) in aggregate data. These identification results
highlight the importance of aggregate-level covariates \(Z\) in making
the identifying assumptions more plausible. These results also make
clear the role that the number of individuals in each aggregation unit
plays in identification and estimation, an aspect of the problem that
prior literature has largely ignored.
To address the restrictive parametric specifications, we propose in
Section~\ref{sec-est} a new double/debiased machine learning estimator
for \(\operatorname{\mathbb{E}}[Y\mid X]\) that allows researchers to minimize bias by
controlling for many covariates. Compared to existing estimation
approaches, our proposed semiparametric estimator is statistically and
computationally efficient without making strong parametric assumptions,
and achieves good accuracy and coverage in practice, as we demonstrate
in simulations and validation on real-world data where ground truth is
available (Section~\ref{sec-valid}).
Key to the estimator's development is a result on the conditional
expectation function (CEF) \(\operatorname{\mathbb{E}}[\overline Y\mid \overline X, Z]\) of the aggregate
outcome: under the identification assumptions, we show that the CEF
takes a partially linear form. This connection also enables application
of the Riesz representation theorem, which highlights the importance of
an additional \emph{positivity} assumption for ecological inference that
is rarely recognized in the literature. Analogously to causal inference,
positivity essentially requires sufficient variation in \(\overline X\) after
controlling for covariates.
We also introduce three new tools that make ecological inference more
useful and more reliable (Section~\ref{sec-extend}). First, in addition
to our estimation theory for the so-called \emph{global} estimand
\(\operatorname{\mathbb{E}}[Y\mid X]\), we develop asymptotically valid confidence intervals
for the \emph{local} estimands \(\operatorname{\mathbb{E}}[Y\mid X, G=g]\), the conditional
means within each aggregation unit \(g\). While not point-identified,
these local estimands are often of interest to practitioners. Second, we
introduce both a sensitivity analysis and a hypothesis test for the key
identifying assumption. These are particularly valuable for applied
researchers, who often are concerened about the plausibility of the
identifying assumption in practical settings. We adapt a sensitivity
analysis framework in causal inference \citep{chernozhukov2022sens} to
provide the first sensitivity analysis for aggregate data. Third, we
provide in the appendix a test for the identifying assumption that draws
on the testable implication of a partially linear CEF. This test, which
has no analogue in the missing data context, is nevertheless approximate
and has limited power, and so we recommend use of the sensitivy analysis
primarily.
Our proposed semiparametric estimator and these three tools are all
implemented in open-source software \citep{seine}, which we apply in a
demonstration in Section~\ref{sec-appl} to the air pollution data of
\citet{jbaily2022air}. Together with the formalization of identifying
assumptions, these methods place ecological inference on a more
practical and robust foundation.
\section{Identification}\label{sec-ident}
Consider a population of exchangeable individuals \(i=1,\dots,n\), each
belonging to an aggregation unit \(G_i\in\mathcal{G}\), which we will refer to as
\emph{geographies} herein for simplicity, since in most applications the
aggregation units correspond to geographic areas. We write the
population of each geography as \(N_g\coloneq |\{i:G_i=g\}|\). Each
individual has a continuous outcome variable \(Y_i\in\ensuremath{\mathbb{R}}\) and a
categorical predictor variable \(X_i\in\{0, 1\}^d\) with
\(d\coloneq |\mathcal{X}|\) levels, represented as a vector of \(d\) mutually
exclusive indicator variables for each possible level in \(\mathcal{X}\). For
cases where \(Y\) is discrete, one can apply the methods here to the
indicator variable for each level of \(Y\) separately. We assume
throughout that \(\operatorname{\mathbb{E}}[Y_i^2]<\infty\).
\subsection{Aggregation model and
estimand}\label{aggregation-model-and-estimand}
Rather than observing \(X_i\) and \(Y_i\) for each individual, the
researcher observes the aggregated variables \[
\overline{Y}_g \coloneq \frac{1}{N_g}\sum_{i\,:\,G_i=g} Y_i \qand
\overline{X}_g \coloneq \frac{1}{N_g}\sum_{i\,:\,G_i=g} X_i.
\]
One of the challenges in studying inference is the need to work with
both individual-level and aggregate-level data simultaneously;
observations that are i.i.d. at one level are not i.i.d. at the other
level, in general. To aid in working across levels, we introduce random
indices over individuals and geographies, which will allow us to
compactly write aggregations and regressions as expectations over these
random indices. For \(\omega\in\ensuremath{\mathbb{R}}^n\) an arbitrary vector of individual
weights with \(\operatorname{\mathbb{E}}[\omega_i]=1\), define a random index \(I^\omega\) by
\(\mathbb{P}(I^\omega=i)=\omega_i / n\). Then we can define
\(G^\omega\coloneq G_{I^\omega}\) to be a random index over geographies.
We will primarily work with two special cases. First, when
\(\omega_i\propto N_{G_i}^{-1}\), so that \(G^\omega\) is uniform over
the geographies, we will drop the superscript and simply write \(G\).
Second, when all \(\omega_i=1\), we will use \(I^n\) and \(G^n\), so
\(I^n\) is uniform over the \(n\) individuals, and \(\mathbb{P}(G^n=g)\) is
proportional to \(N_g\).
In this notation, we may write \(\overline Y_g=\operatorname{\mathbb{E}}_n[Y_I\mid G=g]\)
\footnote{We could have equivalently used \(I_n\) here, since both \(I\)
and \(I^n\) are uniform conditional on geography.} and the global mean
\(\overline Y=\operatorname{\mathbb{E}}_n[Y_{I^n}]=\operatorname{\mathbb{E}}_n[{\overline Y}_{G^n}]\), where \(\operatorname{\mathbb{E}}_n\) denotes an
expectation over the empirical measure. In addition to being more
compact, the random index notation will permit us to state the main
identification results without making an i.i.d. or superpopulation
assumption about the aggregate-level data.\footnote{ The estimation
results require an asymptotic framework, for which we adopt an i.i.d.
model of geographies.}
Associated with each observed \((\overline X_g, \overline Y_g)\) is a vector of
unobserved regression coefficients \[
B_g \coloneq \operatorname{\mathbb{E}}_n[X_I X_I^\top\mid G=g]^{-1}\operatorname{\mathbb{E}}_n[X_I Y_I\mid G=g],
\] which represent the (sample) mean value of \(Y\) for each group in
\(\mathcal{X}\). By definition, \(\overline Y_g\), \(\overline X_g\), and \(B_g\) are
connected by the law of total expectation, traditionally referred to in
ecological inference as the \emph{accounting identity}:
\begin{equation}\protect\phantomsection\label{eq-acct-id}{
\overline Y_g = B_g^\top\overline X_g.
}\end{equation}
Finally, there may be covariates \(Z_g\) available at the geography
level. Together, \((Z_g, \overline X_g, B_g)\) are the full data at the
aggregate level; the researcher observes only the coarsened
\((Z_g, \overline X_g, \overline Y_g)\).
The \emph{global estimand} is the vector of individual-level conditional
means \(\beta\),\footnote{Often, the sample equivalent of this estimand,
where \(\operatorname{\mathbb{E}}\) is replaced by \(\operatorname{\mathbb{E}}_n\), is of interest. The
identification arguments go through identically, but proving
estimation results requires a superpopulation and asymptotic
framework.} defined by \[
\begin{aligned}
\beta_j \coloneq&
\operatorname{\mathbb{E}}[Y_{I^n} \mid X_{I^nj}=1]
= \frac{\operatorname{\mathbb{E}}[N_{G_I} Y_I \mid X_{I^nj}=1]}{\operatorname{\mathbb{E}}[N_{G_I}\mid X_{I^nj}=1]} \\
=& \frac{\operatorname{\mathbb{E}}[\overline X_{G^nj} B_{G^nj}]}{\operatorname{\mathbb{E}}[\overline X_{G^nj}]}
= \frac{\operatorname{\mathbb{E}}[N_G \overline X_{Gj} B_{Gj}]}{\operatorname{\mathbb{E}}[N_G \overline X_{Gj}]}.
\end{aligned}
\] Here, we have written both representations---as a conditional average
of the individual \(Y_i\), and as a weighted average over the
\(B_g\)---in terms of both the individual-weighted random indices
\(G^n\) and the geography-weighted random indices \(G\). We next
investigate under what conditions \(\beta\) is identified from the
coarsened data \((Z_g, \overline X_g, \overline Y_g)\).
\subsection{Identification at the aggregate and individual
level}\label{identification-at-the-aggregate-and-individual-level}
Eq.~\ref{eq-acct-id} makes clear the fundamental identification
challenge: each observation \((\overline X_g, \overline Y_g)\) brings with it
\(d\) unknown parameters: the entries of \(B_g\). This makes
Eq.~\ref{eq-acct-id} a type of random-coefficient model, albeit one with
no error term. These models are well-studied
\citep[e.g.,][]{beran1992estimating}, and to identify \(\beta\), some
kind of regularity across the \(B_g\) must be assumed. For example, if
there is no variation in \(B_g\), so that each \(B_g=\beta\), then there
is a single \(d\)-dimensional unknown parameter, which can be estimated
via linear regression. In fact, assuming constancy across \(B_g\) is
stronger than necessary. What is required is that variation in \(B_g\)
be unrelated to variation in \(\overline X_g\); constancy is a special case
of this condition. The following assumption formalizes the condition;
while \citet{beran1992estimating} state the assumption without
covariates, it can be easily weakened to hold conditional on covariates.
\begin{assump}[Coarsening at random, uniform over individuals]{CAR-U}
For all \(\overline x\) and \(z\),
\(\operatorname{\mathbb{E}}[B_{G^n}\mid Z_{G^n}=z, \overline X_{G^n}=\overline x]=\operatorname{\mathbb{E}}[B_{G^n}\mid Z_{G^n}=z]\).\end{assump}
This assumption is a familiar analogue of the ignorability assumption in
causal inference or the missing-at-random assumption in missing data
analysis. As in \citet{heitjan1991ignorability}, coarsening at random
(CAR) means that the variable which determines the amount of coarsening
or information loss, \(\overline X_{G^n}\), is (mean) independent of the
unobserved data \(B_{G^n}\), given covariates. Since \(B_{G^n}\) is not
observed, in general it is not possible to directly check whether
Assumption~\ref{asm-car-u} holds in the data at hand: researchers should
rely on their substantive knowledge. However, as we will see,
Assumption~\ref{asm-car-u} implies a certain modeling restriction which
may be testable from data. We propose a test for the assumption and
discuss its limitations in Appendix~\ref{sec-id-test}.
Because of the weighting by \(N_g\), Assumption~\ref{asm-car-u} is best
interpreted at the individual level: that for an individual \(i\)
selected uniformly at random, knowing the average \(\overline X_{G_i}\) in
their geography \(G_i\) does not change the expectation of the
individual's corresponding \(B_{G_i}\), given the covariates
\(Z_{G_i}\). Since the assumption is an individual-level one stated in
terms of aggregate variables, it may be difficult to interpret. The
following weighted version of the assumption yields a more helpful
interpretation.
\begin{assump}[Coarsening at random]{CAR}
For all \(\overline x\), \(k\), and \(z\),
\(\operatorname{\mathbb{E}}[B_G\mid Z_G=z, \overline X_G=\overline x, N_G=k]=\operatorname{\mathbb{E}}[B_G\mid Z_G=z]\).\end{assump}
This assumption can of course be stated in two stages: first, that
\(\overline X\) is mean-independent of \(B\) given \(Z\) and \(N\), and
second, that \(B\) is mean-independent of \(N\) given \(\overline X\) and
\(Z\). Although Assumption~\ref{asm-car} is slightly stronger than
Assumption~\ref{asm-car-u}, we will use it throughout the rest of the
paper due to its easier interpretability and the flexibility it provides
in estimation. If only Assumption~\ref{asm-car-u} but not
Assumption~\ref{asm-car} holds, then estimation can proceed identically,
but with observations weighted by \(N_g\) throughout.
Previous work which used an aggregate-level setup often took
\(\overline X_g\) and \(N_g\) as fixed, and so did not consider the ways in
which \(N_g\) could be correlated with other variables. For example,
\citet{ansolabehere1995bias} claimed that weighting by \(N_g\) was
necessary for unbiased estimation, but note that in practice weighting
did not seem to make a large difference. This is the case because
weighting is only required when \(N_g\) is related to \(B_g\) even after
controlling for \(\overline X_g\) and \(Z_g\). As the next result shows, in
general, either Assumption~\ref{asm-car-u} or Assumption~\ref{asm-car}
is sufficient for identification of \(\beta\). All proofs are deferred
to Appendix~\ref{sec-app-proofs}.
\begin{theorem}[Nonparametric
identification]\protect\hypertarget{thm-id}{}\label{thm-id}
For all \(j\in\mathcal{X}\), under Assumption~\ref{asm-car-u}, \[
\beta_j = \frac{\operatorname{\mathbb{E}}[\overline X_{G^nj}\operatorname{\mathbb{E}}[\overline Y_{G^n}\mid Z_{G^n}, \overline X_{G^nj}=1]]}{\operatorname{\mathbb{E}}[\overline X_{G^nj}]},
\] and under Assumption~\ref{asm-car}, \[
\beta_j = \frac{\operatorname{\mathbb{E}}[N_G\overline X_{Gj}\operatorname{\mathbb{E}}[\overline Y_G\mid Z_G, \overline X_{Gj}=1]]}{\operatorname{\mathbb{E}}[N_G\overline X_{Gj}]},
\]
\end{theorem}
Note that when there are no covariates, i.e., \(Z\) is null,
Assumption~\ref{asm-car} is strong and implausible. In the air pollution
example, Assumption~\ref{asm-car} without covariates would imply that a
low-income resident of Los Angeles and a low-income resident of rural
Montana would have the same average \(\text{PM}_{2.5}\) exposure. In a
setting where \(Y\) is vote choice and \(X\) is race,
Assumption~\ref{asm-car} without covariates would imply that white
voters' preferences are identical between Seattle, Wash. and a heavily
Republican rural county such as Fairmount, Ga. Thus, in most
applications, it will be critical to include relevant covariates that
explain variation in the \(B_g\), so that Assumption~\ref{asm-car} is
more plausible. Where covariates are available at the individual level,
they can be aggregated, either marginally or jointly, to form \(Z_g\).
For example, in the voting setting, individual-level age and sex may be
available from the voter file, and their contingency table at the
precinct level could be included as a covariate.
One additional difficulty in evaluating the plausibility of
Assumption~\ref{asm-car} is that it is stated in terms of the aggregated
data itself, while the estimand itself is defined at the individual
level. In some contexts, it may be more straightforward to make
identifying assumptions at the individual level. As
Theorem~\ref{thm-car-ind-agg} records, the following assumption is
sufficient for Assumption~\ref{asm-car}.
\begin{assump}[Coarsening at random at the individual level]{CAR-IND}
For every individual \(i\), and for each \(g\), \(x\), and \(z\),
\(\operatorname{\mathbb{E}}[Y_i\mid G_i=g, X_i=x, Z_{G_i}=z]=\operatorname{\mathbb{E}}[Y_i\mid X_i=x, Z_{G_i}=z]\).\end{assump}
\begin{theorem}[Identification at individual
level]\protect\hypertarget{thm-car-ind-agg}{}\label{thm-car-ind-agg}
Assumption~\ref{asm-car-ind} \(\implies\) Assumption~\ref{asm-car}.
\end{theorem}
Because \(G\) appears directly in Assumption~\ref{asm-car-ind}, it may
be particularly helpful when researchers have substantive knowledge of
the process that assigns individuals to aggregation units (geographies).
However, it is a stronger assumption than Assumption~\ref{asm-car}.
There may be situations where Assumption~\ref{asm-car-ind} does not hold
while Assumption~\ref{asm-car} does.
\section{Estimation}\label{sec-est}
In this section, we apply Assumption~\ref{asm-car} and
Theorem~\ref{thm-id} to develop a semiparametrically efficient estimator
for \(\beta\). We begin with an observation about the form of the
conditional expectation function (CEF) \(\gamma_0\) of \(\overline Y\) under
Assumption~\ref{asm-car}:
\begin{equation}\protect\phantomsection\label{eq-cef}{
\gamma_0(\overline X, Z) := \operatorname{\mathbb{E}}[\overline Y_G\mid Z_G, \overline X_G]
= \operatorname{\mathbb{E}}[B_G^\top\overline X_G\mid Z_G, \overline X_G]
= \eta_0(Z_G)^\top \overline X_G,
}\end{equation} where \(\eta_0(Z_G):=\operatorname{\mathbb{E}}[B_G\mid Z_G]\). Thus, without
any parametric assumptions, \(\gamma_0\) belongs to a restricted class
of partially linear functions \[
\Gamma := \{(\overline x, z) \mapsto \eta(z)^\top \overline x: \{\eta_j\}_{j\in\mathcal{X}}\in L^2(Z)\}.
\] Clearly, \(\Gamma\) is a linear subspace of \(L^2(Z_G, \overline X_G)\);
below, we will show that under an additional assumption, \(\Gamma\) is
in fact a \emph{closed} linear subspace. First, however, we discuss
estimation when Assumption~\ref{asm-car} holds without covariates.
\subsection{Estimation without
covariates}\label{estimation-without-covariates}
The first ecological inference methods were based on simple linear
regression of \(\overline Y_G\) on \(\overline X_G\)
\citep{goodman1953ecological, goodman1959some}. When \(Z\) is null, we
can express \(Y_G\) as \[
\overline Y_G =\eta^\top \overline X_G + \varepsilon_G^\top\overline X_G,
\] where \(\varepsilon_G=B_G-\operatorname{\mathbb{E}}[B_G\mid Z_G]\) is the projection residual from
Eq.~\ref{eq-cef}. Because
\(\operatorname{\mathbb{E}}[B_G\mid Z_G]=\operatorname{\mathbb{E}}[B_G\mid Z_G, \overline X_G, N_G]\), \(\varepsilon_G\) is
orthogonal to any function of \((Z_G, \overline X_G, N_G)\), and we have
immediately that \[
\operatorname{\mathbb{E}}[\varepsilon_G^\top\overline X_G \mid\overline X_G]
=\operatorname{\mathbb{E}}[\varepsilon_G^\top \mid\overline X_G]\overline X_G=0.
\] Thus \(\eta\) can be estimated efficiently by least squares, and
since it is constant, \(\beta=\eta\). When only
Assumption~\ref{asm-car-u} holds, the least-squares regression must be
weighted by \(N_G\), optionally multiplied by a function of
\(\overline X_G\), to guarantee unbiasedness; when Assumption~\ref{asm-car}
holds, any weights which are a function of \(N_G\) and \(\overline X_G\) can
be used. Slightly weaker conditions for finite-sample unbiasedness of
least squares in this setting are possible; see
\citet{ansolabehere1995bias} for an analysis when \(d=2\).
When \(Y\) is binary, so \(\overline Y\) is bounded, a least squares
estimator does not incorporate information contained in these bounds,
which \citet{duncan1953alternative} and \citet{king1997solution} argue
can be substantial. When \(d=2\), \citet{king1997solution} explicitly
models the random coefficients \(B_G\) in order to incorporate the
bounds on \(\overline Y_G\). Specifically, he takes
\(B_G=\eta+\varepsilon_G\sim \operatorname{\mathcal{N}}_{[0,1]^2}(\mu,\Sigma)\), where the subscript
indicates truncation to the unit square. This ensures that
\(0\le \overline Y_G\le 1\), and allows for Bayesian inference for each
\(B_G\). However, it does impose a strong parametric assumption on the
distribution of \(\varepsilon_G\).
Thus it is clear that both Goodman's regression and King's method are
fully consistent with the accounting identity Eq.~\ref{eq-acct-id}, and
they both implicitly assume Assumption~\ref{asm-car} holds
unconditionally. The key difference is in the treatment of the error
term \(\varepsilon_G\): Goodman's regression is semiparametric, in that it is
agnostic to the distribution of the error term; King, by contrast, makes
a distributional assumption. However, this assumption provides several
benefits: while Goodman regression only estimates
\(\varepsilon_G^\top\overline X_G\), King's estimates \(\varepsilon_G\) directly and
ensures it respects any bounds on \(Y\). This allows for estimates of
each geography's \(B_G\) which are consistent with the accounting
identity Eq.~\ref{eq-acct-id} and may be partially identified due to
bounds on \(Y\).
King's model runs into computational difficulties when \(d>2\), since
the normalizing constant and moments of a truncated Normal distribution
are not easily available in higher dimensions \citep[see][ for
details]{km2025review}. Moreover, while King allows for the linear
inclusion of a covariate \(Z\), he does not discuss the challenges of
modeling \(Z\) flexibly, as is required to avoid misspecification bias.
\subsection{Semiparametric estimation with
covariates}\label{sec-est-semi}
By Theorem~\ref{thm-id}, we can write the global estimand using the
notation in Eq.~\ref{eq-cef} as
\begin{equation}\protect\phantomsection\label{eq-id-functional}{
\beta_j = \operatorname{\mathbb{E}}[\gamma_0(e_j, Z_G)u(N_G, \overline X_{Gj})],
}\end{equation} where \(e_j\) is a standard basis vector and
\(u(N_G, \overline X_{Gj})=N_GX_{Gj}/\operatorname{\mathbb{E}}[N_GX_{Gj}]\) weights by the size of
group \(j\) in each geography. Our overall estimation strategy,
following \citet{chernozhukov2022riesz}, is to rewrite \(\beta_j\) in
Neyman-orthogonal form using a Riesz representation of \(\beta_j\),
estimate nuisance functions flexibly using a series estimator, and then
combine the nuisance estimates to form a semiparametrically efficient
estimator for \(\beta_j\).
Consider the Eq.~\ref{eq-id-functional} as a mapping \(\Gamma\to\ensuremath{\mathbb{R}}\). We
can easily see that \(\gamma\mapsto \gamma(e_j, Z)u(N, \overline X_j)\) is
linear in \(\gamma\), so \(\beta_j\) is a linear functional of
\(\gamma\). To apply the Riesz representation theorem to this
functional, we require two additional assumptions.
\begin{assump}[Positivity]{POS}
The random variables \((\overline X_G, Z_G)\) have joint density
\(f(\overline x, z)\) with respect to a dominating measure on
\(\Delta^d\times\mathrm{supp}(Z_G)\), where \(\Delta^d\) is the
\(d\)-dimensional simplex, and there exists a \(\delta>0\) such that
\(f(\overline x, z)>\delta f(z)\) for all
\((\overline x, z)\in\Delta^d\times\mathrm{supp}(Z_G)\), where \(f(z)\) is the
marginal density of \(Z_G\).\end{assump}
\begin{assump}[Bounded N]{BND}
There exists a \(C_N<\infty\) with \(N_GX_{Gj}\le C_N\operatorname{\mathbb{E}}[N_GX_{Gj}]\) for
each \(j\).\end{assump}
Assumption~\ref{asm-bnd} bounds \(u\), ensuring that no single geography
dominates the estimand. Assumption~\ref{asm-pos}, while more involved to
state, essentially requires that there be sufficient variation in
\(\overline X_G\) after controlling for \(Z_G\). It is analogous to the
positivity or overlap assumption in causal inference. Note that the
existence of a joint density \(f(x, z)\) and its positivity at the
vertices of the simplex is also sufficient for the uniqueness of the
conditional expectations \(\gamma(e_j, Z)\), which are evaluated at a
measure-zero set; Assumption~\ref{asm-pos} is in practical terms,
therefore, a mild strengthening of this requirement.
These two assumptions establish two key results: first, that the
partially-linear function class \(\Gamma\) is a \emph{closed} linear
subspace of \(L^2(\overline X_G, Z_G)\); and second, that \(\beta_j\) is
mean-square continuous in \(\gamma\). In what follows, we write
\(\norm{\cdot}\) for the \(L^2(\overline X_G, Z_G)\) norm. These results
yield an immediate corollary, due to the Riesz representation theorem.
\begin{proposition}[]\protect\hypertarget{prp-subspace}{}\label{prp-subspace}
Under Assumption~\ref{asm-pos}, \(\Gamma\) is a closed linear subspace
of \(L^2(\overline X_G, Z_G)\).
\end{proposition}
\begin{proposition}[]\protect\hypertarget{prp-msc}{}\label{prp-msc}
Under Assumption~\ref{asm-pos} and Assumption~\ref{asm-bnd}, for each
\(x\in\mathcal{X}\) the mapping
\(\gamma\mapsto\operatorname{\mathbb{E}}[\gamma(e_j, Z_G)u(N_G, \overline X_{Gj})]\) is mean-square
continuous in \(\gamma\), i.e.,
\(\operatorname{\mathbb{E}}[\gamma(e_j, Z_G)u(N_G, \overline X_{Gj})]\le C_{\mathbb{P}}\norm{\gamma}^2\)
for a \(C_{\mathbb{P}}<\infty\) depending on \(\mathbb{P}\).
\end{proposition}
\begin{corollary}[]\protect\hypertarget{cor-rr}{}\label{cor-rr}
For each \(j\in\mathcal{X}\), there exists a unique
\(\alpha_{0j}(\overline X_G, Z_G)=\zeta_{0j}(Z_G)^\top\overline X_G \in\Gamma\)
with \(\norm{\alpha_{0j}}^2<\infty\) and satisfying
\(\beta_j=\operatorname{\mathbb{E}}[\alpha_{0j}\gamma_0]=\operatorname{\mathbb{E}}[\alpha_{0j} \overline Y]\).
\end{corollary}
We refer to \(\alpha_{0j}\) as the \emph{Riesz representer} of
\(\beta_j\); critically, it also belongs to the restricted class
\(\Gamma\). While it is defined implicitly, in
Appendix~\ref{sec-app-riesz} we present a closed-form expression for
\(\alpha_{0j}\) as a weighted log-derivative of the conditional density
\(f(\overline x\mid z)\). We further discuss its interpretation in the
context of our sensitivity analysis in Section~\ref{sec-sens}.
The second moment of the Riesz representer is tied directly to the
modulus of continuity that establishes Proposition~\ref{prp-msc}, which
in turn hinges critically on Assumption~\ref{asm-pos}. The larger the
second moment of \(\alpha_{0j}\), the less variation in \(\overline X_{gj}\)
there is conditional on \(Z_g\), and the greater the risk of
Assumption~\ref{asm-pos} not holding. This is analogous to causal
inference, where the distribution of the propensity scores plays a
similar role in assessing overlap.
Corollary~\ref{cor-rr} implies that we could estimate \(\beta_j\) in two
ways: either by estimating \(\gamma_0\) and then plugging into
Eq.~\ref{eq-id-functional}, or by estimating \(\alpha_{0j}\) and
plugging into \(\operatorname{\mathbb{E}}[\alpha_{0j} \overline Y]\). However, since both
\(\gamma_0\) and \(\alpha_{0j}\) are functions which must be estimated,
either approach can lead to significant regularization biases. Instead,
a Neyman-orthogonal representation of \(\beta_j\) can be formed based on
the efficient influence function of \(\beta_j\)
\citep{newey1994asymptotic}:
\begin{equation}\protect\phantomsection\label{eq-neyman}{
\beta_j(\gamma_0,\alpha_{0j}) = \operatorname{\mathbb{E}}[\gamma_0(e_j, Z_g)u(N_G,\overline X_G) +
\alpha_{0j}(\overline X_G, Z_G)(\overline Y_G - \gamma_0(\overline X_G, Z_G))].
}\end{equation} This representation is robust to small errors in either
nuisance function, in the sense that its Gateaux derivative with respect
to the nuisance functions vanishes, i.e.,
\(\partial_\gamma \beta_j(\gamma_0,\alpha_{0j}) = \partial_{\alpha_j} \beta_j(\gamma_0,\alpha_0) = 0\)
\citep{chernozhukov2022riesz}. This Neyman orthogonality property is
closely related to double robustness: in fact, if \(\gamma\) is properly
specified in Eq.~\ref{eq-neyman}, then even with a misspecified
\(\alpha\), the score in Eq.~\ref{eq-neyman} is still unbiased for
\(\beta_j\), and vice versa. This can be easily seen by applying
iterated expectations and the representing property of \(\alpha_{0j}\).
Let \(\widehat{\gamma}_m\) and \(\widehat{\alpha}_{mj}\) be estimates of
\(\gamma_0\) and \(\alpha_{0j}\) based on \(m:=|\mathcal{G}|\) geographies. Then
the proposed estimator for \(\beta_j\) is
\begin{equation}\protect\phantomsection\label{eq-beta-hat}{
\begin{aligned}
\widehat\beta_{mj} = \beta_{mj}(\widehat\gamma_m, \widehat\alpha_{mj})
&:= \frac{1}{m}\sum_{g\in\mathcal{G}} \psi_{gj}(\overline Y_g, N_g, \overline X_g, Z_g, \widehat\gamma_m, \widehat\alpha_{mj});\\
\psi_{gj}(\overline Y_g, N_g, \overline X_g, Z_g, \widehat\gamma_m, \widehat\alpha_{mj})
&:= \widehat\gamma_m(e_j, Z_g)u(N_g,\overline X_g) +
\widehat\alpha_{mj}(\overline X_g, Z_g)(\overline Y_g - \widehat\gamma_m(\overline X_g, Z_g)).
\end{aligned}
}\end{equation}
In the next subsections, we discuss estimation of the nuisance functions
\(\gamma_0\) and \(\alpha_{0j}\) and the statistical properties of
\(\widehat\beta_{mj}\). As we noted in Section~\ref{sec-ident}, to state the
asymptotic results for the estimation of the nuisance function and for
\(\widehat\beta_{mj}\), we need an asymptotic framework that treats the
geographies as i.i.d. Consequently, in the remainder of the paper we
index the geographies by \(g\) rather than \(G\) and rely on the
following assumption.
\begin{assump}[]{IID}
\((\overline Y_g, \overline X_g, Z_g, N_g)\overset{\text{iid}}{\sim}\mathbb{P}\).\end{assump}
This differs from the random-index setup of Section~\ref{sec-ident},
where the individuals are i.i.d. but the geographies are in general not.
We stress that the assumption is introduced only to allow a meaningful
asymptotic analysis of the estimation procedure; it is not required for
identification. Rather than interpreting the i.i.d. assumption
literally, one should think of the remaining results as meaningfully
characterizing the behavior of our estimators insofar as the set of
observed geographies can be treated \emph{as if} i.i.d. Our validation
studies demonstrate that this is a reasonable assumption in practice.
\subsection{\texorpdfstring{Estimation of \(\gamma_0\) and
\(\alpha_0\)}{Estimation of \textbackslash gamma\_0 and \textbackslash alpha\_0}}\label{sec-nuisance}
As discussed in \citet{chernozhukov2022riesz}, \(\alpha_{0j}\) can be
estimated via the following representation:
\begin{equation}\protect\phantomsection\label{eq-min-alpha}{
\begin{aligned}
\alpha_{0j} &= \arg\min_{\alpha\in\Gamma} \operatorname{\mathbb{E}}[(\alpha - \alpha_{0j})^2]
= \arg\min_{\alpha\in\Gamma} \operatorname{\mathbb{E}}[\alpha^2 - 2\alpha_{0j}\alpha + \alpha_{0j}^2] \\
&= \arg\min_{\alpha\in\Gamma} \operatorname{\mathbb{E}}[\alpha^2(\overline X_g, Z_g) - 2\alpha(e_j, Z_g)u(N_g, \overline X_{Gj})].
\end{aligned}
}\end{equation} where the final step follows by the representation
property of \(\alpha_{0j}\) and since \(\alpha_{0j}\) is fixed.
In principle, \(\gamma_0\) and \(\alpha_{0j}\) could be therefore
estimated using any nonparametric regression or machine learning method,
with \(\gamma_0\) estimated by minimizing a squared-error loss, and
\(\alpha_{0j}\) estimated by minimizing the loss in
Eq.~\ref{eq-min-alpha}. However, these methods would not produce
estimates \(\widehat\gamma_m\) and \(\widehat\alpha_{mj}\) which belong to the
restricted class \(\Gamma\). To fully leverage the restriction to
\(\Gamma\), we propose series estimators that are in \(\Gamma\) by
construction.
A \emph{linear sieve basis} is a sequence \(\{\Phi_m\}_{m=1}^\infty\) of
vectors of uniformly bounded functions
\(\Phi_m=(\phi_{mk}\in L^\infty(Z_g))_{k=1}^{J_m}\) of dimension
\(J_m\). We discuss three possible sieve bases in
Appendix~\ref{sec-bases}, including interactions of polynomials in each
variable, tensor-product splines, and a more recent tensor-product
cosine basis \citep{zhang2023regression}, as well as the number of basis
functions needed to attain the rates discussed here.
By interacting the elements of a particular sieve basis with the
components of \(\overline X\), we form a basis for a subspace of \(\Gamma\):
\[
\Gamma_m := \{(\overline x, z)\mapsto (\overline x\otimes \Phi_m(z))^\top\theta : \theta\in\ensuremath{\mathbb{R}}^{dJ_m}\},
\] where \(\otimes\) is the Kronecker product, i.e., all pairwise
interactions between the elements of \(\Phi_m\) and \(\overline x\). Our
proposed estimators for \(\gamma_0\) and \(\alpha_{0j}\) are then the
series ridge regression estimators
\begin{equation}\protect\phantomsection\label{eq-sieve-est}{
\begin{aligned}
\widehat\gamma_m(\lambda) &:= \arg\min_{\theta}\{ \operatorname{\mathbb{E}}_m[(\overline Y_g -
(\overline X_g\otimes\Phi_m(Z_g))^\top\theta)^2] + \lambda \theta^\top\theta \} \qand \\
\widehat\alpha_{mj}(\lambda) &:= \arg\min_{\theta}\{ \operatorname{\mathbb{E}}_m[
((\overline X_g\otimes\Phi_m(Z_g))^\top\theta)^2
- 2u(N_g, \overline X_{Gj})(e_j\otimes \Phi_m(Z_g))^\top\theta
] + \lambda \theta^\top\theta \},
\end{aligned}
}\end{equation} where \(\operatorname{\mathbb{E}}_m\) denotes the empirical expectation. It is
clear that
\(\widehat\gamma_m(\lambda), \widehat\alpha_{mj}(\lambda)\in\Gamma_m\subseteq\Gamma\).
The closed-form solution for \(\widehat\gamma_m(\lambda)\) is well-known; we
present a closed-form solution for \(\widehat\alpha_{mj}(\lambda)\) in
Appendix~\ref{sec-app-riesz}. To estimate \(\lambda\), we employ
leave-one-out cross-validation (LOOCV) for the loss of \(\widehat\gamma_m\).
The LOOCV loss can be efficiently computed for a range of \(\lambda\)
using the hat matrix and the singular value decomposition of the design
matrix. The asymptotic considerations in \citet{singh2024kernel} suggest
that the same \(\lambda\) can be used for both \(\widehat\gamma_m\) and
\(\widehat\alpha_{mj}\); we do so here, since it is not computationally
convenient to compute LOOCV errors for \(\widehat\alpha_{mj}\). We discuss
in Appendix~\ref{sec-app-impl} several other practical considerations in
implementing the estimator: bounded \(Y\), the use of weights in
estimation, and the choice of sieve basis.
\subsubsection{Convergence rates}\label{convergence-rates}
For the estimate \(\widehat\beta_{mj}\) to be asymptotically normal, we will
need \(\widehat\gamma_m\) and \(\widehat\alpha_{mj}\) to converge quickly enough
to their targets, which requires conditions on the true \(\eta_0(z)\) as
well as on the sieve basis \(\Phi_m\). Sieve estimators for varying
coefficient models have been proposed and analyzed before
\citep{park2015varying}. However, most treatments focus on regression
estimators and do not directly apply to estimating \(\alpha_{0j}\).
Moreover, Assumption~\ref{asm-pos} and Assumption~\ref{asm-bnd} can be
sufficient for various regularity conditions required by other
estimators, and the fact that \(\overline X\) is supported on the simplex
\(\Delta^d\) is specific to this case as well. Thus we state and prove
the necessary conditions and results in full here.
\begin{assump}[Sieve estimation regularity conditions]{SR}
We assume the following about the data-generating process:
\begin{enumerate}
\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item
\(Z_g\) is supported on a bounded subset of \(\ensuremath{\mathbb{R}}^p\) for some \(p\),
and has a bounded density \(0<f(z)<\infty\) with respect to Lebesgue
measure on its support.
\item
\(\operatorname{\mathbb{V}}[Y_g\mid Z_g=z, \overline X_g=\overline x]\) is bounded on
\(\Delta^d\times\mathrm{supp}(Z_g)\).
\item
There exists a function class \(\mathcal{F}\subseteq L^2(Z_g)\) with
\(\eta_0,\zeta_{0j}\in\mathcal{F}^d\) and \(\norm{\eta_{0j}}_\infty<\infty\)
and \(\norm{\zeta_{0jk}}_\infty<\infty\) for each \(j,k\in\mathcal{X}\), where
\(\eta_0\) and \(\zeta_{0j}\) are the true component functions for
\(\gamma_0\) and \(\alpha_{0j}\), respectively.
\end{enumerate}
For the conditions on the sieve basis \(\Phi_m\), let \(\nu\) denote
Lebesgue measure on \(\mathrm{supp}(Z_g)\), rescaled to have unit mass, and let
\[
\rho_{m} := \sup_{f'\in\mathcal{F}}\inf_{f\in\operatorname{span}\Phi_m} \norm{f - f'}_{2,\nu} \qand
A_m := \sup_{f\in\operatorname{span}\Phi_m, \norm{f}_{2,\nu}\neq 0} \norm{f}_\infty / \norm{f}_{2,\nu},
\] where \(\norm{f}_{2,\nu}^2 := \int f^2 d\nu\). We assume
\begin{enumerate}
\setcounter{enumi}{3}
\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item
\(\operatorname{span}\Phi_m\) is identifiable, i.e.,
\(\norm{f}_{2,\nu}=0 \implies f=0\) for \(f\in\operatorname{span}\Phi_m\);
\item
\(A_m\rho_{m}\to 0\); and
\item
\(A_m^2J_m/m\to 0\).
\end{enumerate}\end{assump}
Condition (1) ensures that the \(L^2(\overline X_G,Z_G)\) norm is equivalent
to the Lebesgue norm \(\norm{\cdot}_{2,\nu}\) on \(\mathrm{supp}(Z_G)\), which
means that properties of the sieve basis can be checked on the latter,
independent of the data distribution, and carry over to the former.
Conditions (4)--(6) hold for many common sieve bases and function
classes, as we discuss below. These conditions would generally permit
estimation of functions in \(\mathcal{F}\) at the rate \(\rho_m + \sqrt{J_m/m}\);
we, however, need to estimate functions in the space
\(\Gamma^\mathcal{F} := \{(\overline x, z)\mapsto f(z)^\top\overline x : f\in\mathcal{F}^d\}\subseteq \Gamma\)
which contains \(\gamma_0\) and \(\alpha_{0j}\) under
Assumption~\ref{asm-sr} (3). The next theorem shows that this is
possible at the same rate.
\begin{theorem}[]\protect\hypertarget{thm-sieve}{}\label{thm-sieve}
Under Assumption~\ref{asm-bnd}, Assumption~\ref{asm-pos}, and
Assumption~\ref{asm-sr}, \(\widehat\gamma_m\) and \(\widehat\alpha_{mj}\) exist
uniquely with probability approaching one as \(m\to\infty\), and for all
\(j\in\mathcal{X}\), we have for the unpenalized estimators that \[
\norm{\widehat\gamma_m(0) - \gamma_0} = O_{\mathbb{P}}(\rho_m + \sqrt{J_m/m}) \qand
\norm{\widehat\alpha_{mj}(0) - \alpha_{0j}} = O_{\mathbb{P}}(\rho_m + \sqrt{J_m/m}).
\] If additionally, the eigenvalues of the matrix
\(\Phi_m(\vb Z)^\top\Phi_m(\vb Z)\) (where \(\vb Z\) is the matrix of
\(Z_g\)) are uniformly bounded away from zero, and
\(\lambda_m=O(\sqrt{J_m/m})\), then \[
\norm{\widehat\gamma_m(\lambda_m) - \gamma_0} = O_{\mathbb{P}}(\rho_m + \sqrt{J_m/m}) \qand
\norm{\widehat\alpha_{mj}(\lambda_m) - \alpha_{0j}} = O_{\mathbb{P}}(\rho_m + \sqrt{J_m/m}).
\]
\end{theorem}
The rate on \(\lambda_m\) ensures that the penalty does not bias the
estimator asymptotically. In practice, as mentioned above, we pick the
penalty using LOOCV. When the design matrix satisfies certain
conditions, such as being a linear transformation of i.i.d. variables,
\citet{patil2021uniform} show that LOOCV is uniformly consistent for the
optimal penalty, even when the number of predictors grows with the
sample size. While their setup differs from the one here, LOOCV is very
likely to reduce the estimation error in finite samples from the
unpenalized estimator, and so \(\widehat\gamma_m(\widehat\lambda_\mathrm{LOO})\) and
\(\widehat\alpha_{mj}(\widehat\lambda_\mathrm{LOO})\) should converge at the same
rate.
\subsection{Properties of the proposed estimator}\label{sec-conv}
The proposed estimator \(\widehat\beta_{mj}\) in Eq.~\ref{eq-beta-hat} is a
double/debiased machine learning (DML) estimator
\citep{chernozhukov2022riesz}, but unlike many DML estimators, the
nuisance functions \(\gamma_0\) and \(\alpha_{0j}\) are estimated on the
same data as \(\widehat\beta_{mj}\) rather than being cross-fitted. While
cross-fitting makes theoretical analysis easier, in finite samples cross
fitting over just \(k\) folds can substantially increase variance. Some
authors recommend generating multiple sets of \(k\) folds, which reduces
variance but increases computational cost. Here, \(\gamma_0\) and
\(\alpha_{0j}\) are estimated using a ridge penalty, and so we might
expect better behavior than an estimator using a black-box machine
learning method that could be arbitrarily sensitive to the data used to
fit it. We therefore establish the asymptotic normality and
semiparametric efficiency of \(\widehat\beta_{mj}\) without cross-fitting,
using the results of \citet{chen2022debiased}, which rely on an
algorithmic stability condition that is satisfied by the ridge penalty
used here.
\begin{theorem}[]\protect\hypertarget{thm-dml}{}\label{thm-dml}
Suppose that \(\widehat\gamma_m\) and \(\widehat\alpha_{mj}\) are the
ridge-penalized series estimators in Eq.~\ref{eq-sieve-est} that achieve
estimation error rates
\(\norm{\widehat\gamma_m - \gamma_0} = o_{\mathbb{P}}(m^{-1/4})\) and
\(\norm{\widehat\alpha_{mj} - \alpha_{0j}} = o_{\mathbb{P}}(m^{-1/4})\) for each
\(j\in\mathcal{X}\). Assume also that \(\lambda\asymp \sqrt{J_m/m}\), that the
eigenvalues of the matrix \(\Phi_m(\vb Z)^\top\Phi_m(\vb Z)\) are
uniformly bounded away from zero, and that \(\norm{\overline Y}_{2r}<\infty\)
for some \(r>1\). Then \(\widehat\beta_m\) is asymptotically normal with
limiting distribution \[
\sqrt{m}(\widehat\beta_m - \beta) \Rightarrow \operatorname{\mathcal{N}}(0, \operatorname{\mathbb{E}}[\psi_g\psi_g^\top]),
\] where \(\psi_g\) is the vector of scores \(\psi_{gj}\) for each
\(j\in\mathcal{X}\).
\end{theorem}
The semiparametric efficiency of \(\widehat\beta_m\) follows because
\(\psi_g\) is the efficient influence function for \(\beta\) and
\(\operatorname{\mathbb{E}}[\psi_g\psi_g^\top]\) is the semiparametric efficiency bound.
Theorem~\ref{thm-dml} means that in practice we can easily construct
asymptotically valid confidence regions for \(\beta\) using the sample
covariance of the scores \(\psi_g\). This also allows for asymptotically
valid confidence intervals for linear contrasts of \(\beta\), which are
often of interest in applications measuring differences in \(Y\) between
groups.
\section{Further Extensions}\label{sec-extend}
In this section, we discuss two extensions of the proposed method:
estimation of the \emph{local} quantities \(B_g\) for each geography
\(g\), and a sensitivity analysis for violations of
Assumption~\ref{asm-car}. A third extensions, a possible hypothesis test
for Assumption~\ref{asm-car}, is presented in
Appendix~\ref{sec-id-test}.
\subsection{Local estimates}\label{sec-local}
Often, researchers are interested not just in the global estimand
\(\beta=\operatorname{\mathbb{E}}[Y\mid X]\), but how this relationship varies by geography.
For example, in political science, \(\beta\) may describe the voting
preference (\(Y\)) of a racial group (\(X\)) nationally, but researchers
may also be interested in how this relationship varies across counties
or precincts. These local quantities are exactly the missing data
\(B_G\). Since there is a single (unobserved) \(B_G\) per geography, it
is of course not possible to consistently estimate the \(B_G\)
themselves. However, it is possible to construct valid confidence
regions \(B_G\), under additional assumptions.
Under Assumption~\ref{asm-car}, \(B_G := \eta_0(Z_G) + \varepsilon_G\), with
\(\varepsilon_G\) mean-zero and mean-independent of \((N_G, \overline X_G, Z_G)\).
Thus, a natural point estimate for \(B_G\) is \(\widehat\eta(Z_G)\). The
more variation in \(B_G\) (and thus \(\overline Y_G\)) explained by \(Z_G\),
the more accurate this point estimate will be. However, it will not be
consistent for \(\beta_G\) since \(\varepsilon_G\) has non-zero variance.
Additionally, while \(B_G\) must satisfy the accounting identity
(Eq.~\ref{eq-acct-id}), i.e., \(\overline Y_G=B_G^\top \overline X_G\), in general
the estimates \(\widehat\eta(Z_G)\) will not. We aim to develop a point
estimate and confidence region for \(B_G\) that addresses these two
issues. Doing so will require consistently estimating the covariance
matrix \(\operatorname{\mathbb{V}}[\varepsilon_G]\), which requires additional assumptions.
\begin{assump}[Coarsening at random, second moments]{CAR2}
For all \(\overline x\), \(k\), and \(z\),
\(\operatorname{\mathbb{E}}[B_GB_G^\top\mid X_G=\overline x, Z_G=\overline z, N_G=k]=\operatorname{\mathbb{E}}[B_GB_G^\top\mid Z_G=z]\).\end{assump}
Assumption~\ref{asm-car2} could be equivalently written in terms of the
residuals \(\varepsilon_G\). Of course, both Assumption~\ref{asm-car2} and
Assumption~\ref{asm-car} are implied by the stronger condition that
\(B_G\) is conditionally independent of \(N_G\) and \(\overline X_G\) given
\(Z_G\). A version of Assumption~\ref{asm-car2} that does not condition
on \(N_G\) could also be applied, analogously to
Assumption~\ref{asm-car-u}.
Denote the covariance matrix as \(\Sigma\) and let
\(\Sigma(z):=\operatorname{\mathbb{V}}[\varepsilon_G\mid Z_G=z]\). Under Assumption~\ref{asm-car2},
\(\Sigma(z)\) fully describes the conditional variance structure of
\(\varepsilon_G\). Additionally, let
\(\kappa_0(x,z) :=\operatorname{\mathbb{E}}[(\overline Y_G - \operatorname{\mathbb{E}}[\overline Y_G \mid \overline X_G,Z_G])^2\mid \overline X_G=x,Z_G=z]\).\\
We then have the following identification result.
\begin{proposition}[]\protect\hypertarget{prp-var-id}{}\label{prp-var-id}
Under Assumption~\ref{asm-car2}, for any \(j,k\in\mathcal{X}\) we have
\(\Sigma_{jk}(z) = 2(\kappa_0(\tfrac{1}{2} e_j + \tfrac{1}{2} e_k, z) - \tfrac{1}{4}\kappa_0(e_j, z) - \tfrac{1}{4}\kappa_0(e_k, z))\).
\end{proposition}
The form of this result is due to the polarization identity for
recovering a bilinear form from a quadratic form. The result in
Proposition~\ref{prp-var-id} means that a consistent estimate of
\(\kappa_0\) can be used along with a consistent estimate of
\(\gamma_0\) (which includes \(\eta_0\)) to form an asymptotically valid
confidence region for \(B_G\) using the multivariate Chebyshev
inequality. As discussed above, however, the point estimate
\(\widehat\eta(Z_G)\) will not in general satisfy the accounting identity.
To further improve the point estimate and confidence region, we can
project the estimate onto the \(d-1\)-dimensional region implied by the
accounting identity. Specifically, let
\(H(\overline x, \overline y) := \{b\in \ensuremath{\mathbb{R}}^d : b^\top \overline x = \overline y\}\) be the
set of possible values of \(B_G\) given \(\overline X_G=\overline x\) and
\(\overline Y_G=\overline y\) for an unbounded \(\overline Y_G\); we let \([H_g]\)
denote the matrix with columns forming an orthonormal basis
\(H(\overline X_g, \overline Y_g)\). Such a basis can be efficiently computed by
taking the QR decomposition of the block matrix \((\overline X_g\quad I_d)\)
and discarding the first column of \(Q\). Also let \(\widehat\Sigma(z)\) be
the estimate of \(\Sigma(z)\) obtained from the expression in
Proposition~\ref{prp-var-id} using an estimated
\(\widehat\kappa\).\footnote{ In finite samples, the \(\widehat\Sigma(z)\)
yielded by Proposition~\ref{prp-var-id} may not be positive
semidefinite. In these cases, we recommend projecting
\(\widehat\Sigma(z)\) onto the space of positive semidefinite matrices,
e.g., by setting all negative eigenvalues to zero. This will not
affect the consistency of \(\widehat\Sigma(z)\).} Then we can obliquely
project \(\widehat\eta(Z_g)\) onto \(H(\overline X_g, \overline Y_g)\) along
\(\widehat\Sigma(Z_g)\) to obtain a point estimate \(\tilde B_g\) that
satisfies the accounting identity: \[
\widehat B_g := \widehat\Pi_g \widehat\eta(Z_g); \quad
\widehat\Pi_g := [H_g]([H_g]^\top \widehat\Sigma(Z_g)^{-1}[H_g])^{-1}[H_g]^\top \widehat\Sigma(Z_g)^{-1}.
\] where \(\widehat\Pi_g\) is the oblique projection matrix. We can then
further obliquely project \(\widehat B_g\) onto
\(\mathrm{supp}(B_G)=\mathrm{supp}(\overline Y_G)^d\), which is convex, yielding a point
estimate \(\widehat B'_g\) that lies in \[
H'(\overline x, \overline y) := H(\overline x, \overline y)\cap \mathrm{supp}(\overline Y_G)^d.
\] Of course, when \(Y\) is unbounded, \(H'=H\). Now define a confidence
region \[
\begin{aligned}
{R'_g}^\alpha &:= \{b\in H'(\overline X_g,\overline Y_g) :
(b - \widehat B'_g)^\top (\widehat\Pi_g\widehat\Sigma(Z_g))^+ (b - \widehat B'_g) \le \frac{d-1}{\alpha} \},
\end{aligned}
\] where \(A^+\) denotes the Moore-Penrose pseudoinverse of \(A\). Then
we have the following result.
\begin{theorem}[]\protect\hypertarget{thm-local-ci}{}\label{thm-local-ci}
Suppose that \(\widehat\eta\to^p\eta_0\) and \(\widehat\kappa\to^p\kappa_0\)
pointwise. Then for \(0<\alpha<1\), as \(m\to\infty\),
\(\mathbb{P}(B_g\in {R'_g}^\alpha) \ge 1-\alpha + o(1)\).
\end{theorem}
In practice, we might expect these confidence regions to be
conservative, especially for bounded \(Y\) when a second projection is
used. Additional distributional assumptions on \(\varepsilon_G\), such as
unimodality or Normality, can be used to further tighten the confidence
regions in practice. For example, for confidence intervals for a single
component \(B_{Gj}\), if the distribution of
\(\varepsilon_{G}^\top\overline X_G\mid Z_G\) is assumed unimodal, then the width of
the confidence interval can be reduced by a factor of \(2/3\)
\citep{vysochanskij1980justification}.
\citet{breunig2021varying} discusses estimation of other aspects of the
distribution of \(\varepsilon\) in varying coefficient models, such as higher
moments or quantiles, in more detail, and develops sieve estimators for
efficiently estimating these quantities. These estimators may be applied
to build intervals for \(B_G\), which in some cases may be narrower, but
more work is needed to apply these in a way that respects the accounting
identity, as the intervals developed in this section do.
\subsection{Sensitivity analysis}\label{sec-sens}
Every result so far has relied critically on the
Assumption~\ref{asm-car} assumption. In practice, it is unlikely that
Assumption~\ref{asm-car} holds \emph{exactly}, and so it is important to
understand how sensitive the proposed estimates of \(\beta\) are to
violations of this assumption. To do so, we can apply results of
\citet{chernozhukov2022sens}, who develop a nonparametric sensitivity
analysis for estimands for which a Riesz representer exists.
Rather than assume that Assumption~\ref{asm-car} holds, the sensitivity
analysis assumes that it holds conditional on an unobserved variable
\(A_G\). This is not really an additional assumption, since we can
always take \(A_G=\operatorname{\mathbb{E}}[B_G\mid N_G, \overline X_G, Z_G]\), which then makes
Assumption~\ref{asm-car} conditional on \(A_G\) hold trivially.
Let \(\gamma^A_0\) and \(\alpha^A_{0j}\) be the regression function and
Riesz representer, respectively, defined conditional on \(A_G\). Unlike
\(\gamma_0\) and \(\alpha_{0j}\), which can be estimated from the data,
\(\gamma^A_0\) and \(\alpha^A_{0j}\) cannot be estimated because \(A_G\)
is not observed. However, if Assumption~\ref{asm-car} holds conditional
on \(A_G\) only, then \(\beta_j\) can only be consistently estimated
using \(\gamma^A_0\) and \(\alpha^A_{0j}\). The estimate from the data,
\(\widehat\beta_j\), will converge to some \(\beta^*_j\).
\citet{chernozhukov2022sens} (Theorem 2 and Corollary 2) then establish
the following result.
\begin{theorem}[]\protect\hypertarget{thm-sens}{}\label{thm-sens}
When Assumption~\ref{asm-car} holds conditional on \(A_G\), then
\(|\beta^*_j - \beta_j| \le \rho S C_\gamma C_\alpha\), where \[
\begin{aligned}
\rho &:= |\mathrm{Cor}(\gamma^A_0 - \gamma_0, \alpha^A_{0j} - \alpha_{0j})|, \qquad
S^2 := \operatorname{\mathbb{E}}[(\overline Y - \gamma_0(\overline X, Z))^2]\operatorname{\mathbb{E}}[\alpha_{0j}(\overline X, Z)^2], \\
C_\gamma^2 &:= \frac{\operatorname{\mathbb{E}}[(\gamma^A_0 - \gamma_0)^2]}{\operatorname{\mathbb{E}}[(\alpha^A_{0j} - \alpha_{0j})^2]}
= R^2_{\overline Y \sim A \mid \overline X, Z}, \qand
C_\alpha^2 := \frac{\operatorname{\mathbb{E}}[{\alpha^A_{0j}}^2] - \operatorname{\mathbb{E}}[\alpha_{0j}^2]}{\operatorname{\mathbb{E}}[\alpha_{0j}^2]}
= \frac{1 - R^2_{\alpha^A_{0j}\sim \alpha_{0j}}}{R^2_{\alpha^A_{0j}\sim \alpha_{0j}}}.
\end{aligned}
\]
\end{theorem}
In other words, the bias due to violations of Assumption~\ref{asm-car}
is bounded by the product of four terms. The first, \(\rho\), can be
upper bounded by 1, which represents adversarial confounding. It can
also be benchmarked to observed covariates, as discussed below. The
second, \(S\), is a scaling factor which can be estimated from the data.
The third, \(C_\gamma\), measures the proportion of the residual
variation in \(\overline Y_g\) explained by the unobserved confounder
\(A_G\). The fourth, \(C_\alpha\), decreases with the proportion of the
residual variation in the Riesz representer \(\alpha_{0j}\) explained by
the unobserved confounder \(A_G\).
Theorem~\ref{thm-sens} can also be directly applied to differences of
the form \(\beta_j-\beta_k\) (or, more generally, any linear contrast),
since the Riesz representer for such differences is simply
\(\alpha_{0j}-\alpha_{0k}\). When researchers are primarily interested
in differences between groups, this approach can yield tighter bounds
than applying Theorem~\ref{thm-sens} to each group separately and then
using the triangle inequality.
Researchers can vary the \emph{sensitivity parameters} \(C_\gamma\) and
\(C_\alpha\) to understand how sensitive their estimates are to
violations of Assumption~\ref{asm-car}. In fact, the entire sensitivity
analysis can be visualized on a single plot, by plotting contours of the
bound against \(C_\gamma\) and \(C_\alpha\) as contour lines. This type
of plot is familiar to causal inference researchers, who use it to
visualize sensitivity to confounding in observational studies.
As a minimal alternative to a sensitivity plot, researchers can
calculate the \emph{robustness value}, which measures the minimum
assumption violation (in terms of \(C_\gamma\) and \(C_\alpha\)) needed
to cause a bias of a specified amount. Formally, \(RV(\delta)\) is the
maximum value \(RV\) such that
\(R^2_{\overline Y \sim A \mid \overline X, Z}\le RV\) and
\(1-R^2_{\alpha^A_{0j}\sim \alpha_{0j}}\le RV\) imply
\(|\beta^*_j - \beta_j|< \delta\). In other words, if either
\(R^2_{\overline Y \sim A \mid \overline X, Z}\) and
\(1-R^2_{\alpha^A_{0j}\sim \alpha_{0j}}\) are both smaller than
\(RV(\delta)\), then the bias is less than \(\delta\). Possible values
for \(\delta\) include a certain multiple of the standard error of
\(\widehat\beta_j\), or a substantively meaningful threshold. For example,
in comparing groups \(X=1\) and \(X=2\), \(RV(\widehat\beta_2-\widehat\beta_1)\)
would measure the minimum confounding needed to explain away the entire
estimated difference between the two groups.
For inference, \citet{chernozhukov2022sens} propose a DML estimateof the
bounds, \(\widehat\beta_j \pm \widehat\sigma\widehat\nu|\rho|C_\gamma C_\alpha\),
where \[
\widehat\sigma^2 := \operatorname{\mathbb{E}}_m[(\overline Y_g - \widehat\gamma(\overline X_g, Z_g))^2] \qand
\widehat\nu^2 := \operatorname{\mathbb{E}}_m[2\widehat\alpha_j(e_j, Z_g) - \widehat\alpha_j(\overline X_g, Z_g)^2],
\] Because these use the same nuisance functions \(\widehat\gamma\) and
\(\widehat\alpha_j\), and both estimators are based on Neyman-orthogonal
representations, these estimates will be semiparametrically efficient by
the same argument as for Theorem~\ref{thm-dml} under slightly modified
regularity conditions. The full conditions are stated in
\citet{chernozhukov2022sens}, who also propose DML confidence bounds for
these bounds which involve further computation.
\subsubsection{Interpretation}\label{interpretation}
Interpreting \(C_\gamma\) is relatively straightforward as a
(nonparametric) partial \(R^2\) of \(\overline Y_G\) on \(A_G\), conditional
on \(\overline X_G\) and \(Z_G\). Interpreting \(C_\alpha\) is more
difficult, since \(\alpha_{0j}\) is defined implicitly by
Corollary~\ref{cor-rr}. Appendix~\ref{sec-app-riesz} derives an explicit
representation of \(\alpha_{0j}\) as a weighted log derivative of the
conditional density of \(\overline X_G\) given \(Z_G\), \[
\alpha_{0j} = -u(N_g, \overline X_{gj})\partial_{\overline x_j}\log f_{\overline x\mid z}(\overline X_G, Z_G),
\] with \(\alpha^A_{0j}\) defined analogously but conditional on \(A_G\)
as well. When \(\overline X_{Gj}\mid Z_G\) is homoskedastic Gaussian, then we
have \[
\alpha_{0j}\propto N_g\overline X_{gj}(\overline X_{gj} - \operatorname{\mathbb{E}}[\overline X_{gj}\mid Z_g]),
\] and if \(\operatorname{\mathbb{E}}[\overline X_{gj}\mid Z_g]\) is not particularly variable
(i.e., \(R^2_{\overline X_{j}\sim Z}\) is small), then \(C_\alpha^2\) is
approximately upper bounded by
\(R^2_{\overline X_j\sim A\mid Z}/(1 - R^2_{\overline X_j\sim A\mid Z})\), which
is increasing in \(R^2_{\overline X_j\sim A\mid Z}\) So, in \emph{very} rough
terms, \(C_\alpha\) measures how much of the variation in
\(\overline X_{Gj}\) is explained by \(A_G\), conditional on \(Z_G\). In
practice, we recommend that researchers benchmark \(C_\alpha\) to
observed covariates to help in judging the plausibility of different
values of \(C_\alpha\). This benchmarking, described in
Appendix~\ref{sec-bench} and demonstrated in the application, is used in
causal inference as well.
\section{Validation}\label{sec-valid}
This section validates the proposed method in a simulation study and on
real-world data where the ground truth is known. In simulations, the
estimator outperforms alternatives, and both the global and local
confidence intervals achieve nominal coverage. In real-world data where
standard linear regression badly misses the ground truth, our method
that models a basis expansion of dozens of covariates reduces the
estimation error to within a percentage point for most groups.
\subsection{Simulation studies}\label{sec-sim}
We examine the performance of our method on data simulated from the data
generating process assumed by the now-standard method of
\citet{king1997solution}. This data-generating process draws \(\overline X\)
and \(Z\) in a correlated manner, with \(\overline X\in\Delta^d\), and then
draws \(B\) conditional on \(Z\) from a Normal distribution truncated to
the unit hypercube, so that each \(B_j\in[0, 1]\). The aggregate outcome
\(\overline Y\) is then directly calculated as \(B^\top \overline X\). For
simplicity, the size of each geography is assumed uniform, i.e.,
\(N=1\). We simulate different levels of confounding by changing both
the correlation between \(\overline X\) and \(Z\), and with the correlation
between \(B\) and \(Z\) fixed at 0.2. The entries in \(B\) are also
correlated, with a pairwise \(R^2=0.25\). Full details of the data
generating process are in Appendix~\ref{sec-valid-detail}.
In the first simulation study, we generate 1,000 datasets with \(m=500\)
geographies, \(d=2\) predictors, \(p=3\) covariates, and moderate
confounding: \(R^2_{B\sim Z}=R^2_{\overline X\sim Z}=0.5\). This 2-by-2
caseallows us to compare our method to existing methods which only
support \(d=2\). On each of the 1,000 data replicates, we applied (1)
our proposed method, including covariates entered linearly, (2) linear
regression without covariates \citep{goodman1953ecological}, (3) the
truncated-normal model of \citet{king1997solution}, from the R package
\texttt{ei}, both with and without covariates, and (4) the
Multinomial-Dirichlet count model of \citet{rosen2001bayesian},
implemented as \texttt{ei.MD.Bayes} in the R package \texttt{eiPack},
both with and without covariates.
\begin{table}
\centering{
\centering
\begin{tabular}[t]{lccccc}
\toprule
& Covariates? & RMSE & Coverage ($50\%$) & Coverage ($95\%$) & Time (s)\\
\midrule
Proposed method & Yes & 0.013 & 0.62 & 0.99 & 0.05\\
Goodman (linear regression) & & 0.054 & 0.11 & 0.31 & 0.01\\
Goodman (linear regression) & Yes & 0.014 & 0.41 & 0.86 & 0.00\\
King (1997; ei) & & 0.051 & 0.08 & 0.22 & 10.00\\
King (1997; ei) & Yes & 0.016 & 0.30 & 0.74 & 83.00\\
RJKT (2002; eiPack) & & 0.064 & 0.48 & 0.86 & 2.40\\
RJKT (2002; eiPack) & Yes & 0.083 & 0.39 & 0.75 & 9.80\\
\bottomrule
\end{tabular}
}
\caption{\label{tbl-sims}\textbf{Comparison of existing methods on
simulated data, 2×2 case}. Root mean squared error (RMSE), coverage of
nominal 50\% and 95\% confidence intervals, and average computation time
(in seconds) for different methods on simulated data. Each data
replicate contained \(m = 500\) precincts. A `X' in the covariates
column indicates that the method controlled for confounding covariates.
RJKT refers to \citet{rosen2001bayesian}.}
\end{table}
The proposed method achieved the lowest root mean square error (RMSE) in
estimating the global parameters \(\beta\), with King's
\citeyearpar{king1997solution} model with covariates a close second.
Table~\ref{tbl-sims} presents the results. The three methods that did
not control for confounding all had similar error, around 3--4 times
higher than the proposed method. The confidence intervals for the
proposed achieved nominal coverage and in fact moderately over-covered.
None of the other methods achieved close to nominal coverage, despite
the data being drawn from a model that is exactly consistent with the
model fit by King's method. Even more concerningly, the model of
\citet{rosen2001bayesian}, which is the only method implemented in
public software that can handle \(d>2\), suffers \emph{higher} error and
\emph{lower} coverage rates when covariates are included. Finally,
estimation in competing methods is two orders of magnitude slower than
our method when covariates are not used, and even more when covariates
are included.
\begin{figure}[t]
\begin{minipage}{0.50\linewidth}
\centering{
\sbox\pandoc@box{\includegraphics[keepaspectratio]{figures/sims-seine.pdf}}
\Gscale@div\@tempa{\textheight}{\dimexpr\ht\pandoc@box+\dp\pandoc@box\relax}
\Gscale@div\@tempb{\linewidth}{\wd\pandoc@box}
\ifdim\@tempb\p@<\@tempa\p@\let\@tempa\@tempb\fi
\ifdim\@tempa\p@<\p@\scalebox{\@tempa}{\usebox\pandoc@box}
\else\usebox{\pandoc@box}
\fi
}
\subcaption{\label{fig-sim-seine}\textbf{Global estimates}}
\end{minipage}
\begin{minipage}{0.50\linewidth}
\centering{
\sbox\pandoc@box{\includegraphics[keepaspectratio]{figures/sims-seine-local.pdf}}
\Gscale@div\@tempa{\textheight}{\dimexpr\ht\pandoc@box+\dp\pandoc@box\relax}
\Gscale@div\@tempb{\linewidth}{\wd\pandoc@box}
\ifdim\@tempb\p@<\@tempa\p@\let\@tempa\@tempb\fi
\ifdim\@tempa\p@<\p@\scalebox{\@tempa}{\usebox\pandoc@box}
\else\usebox{\pandoc@box}
\fi
}
\subcaption{\label{fig-sim-local}\textbf{Local estimates}}
\end{minipage}
\caption{\label{fig-sims}\textbf{Error and coverage on simulated data}.
RMSE and coverage of 95\% nominal confidence intervals for (a)
\textbf{global parameters} \(B\) and (b) \textbf{local parameters}
\(B_G\). Proposed method is run on 1000 replicates of simulations with
different sample sizes, numbers of predictors (columns; \(d\)), number
of covariates (rows; \(p\)), and strength of confounding (colors;
measured as the \(R^2_{\overline X\sim Z}\) with \(R^2_{B\sim Z}\) is fixed
at 0.2).}
\end{figure}
In the second simulation study, we vary
\(m\in\{50, 100, 500, 1\,000, 10\,000\}\) (with \(m=50\) mimicking a
50-state regression), \(d\in\{2,5,10\}\), \(p\in\{1,3,10\}\), and
\(R^2_{\overline X\sim Z}\in\{0, 0.2, 0.5\}\), while fixing
\(R^2_{B\sim Z}=0.2\), with 1,000 simulated datasets for each
combination. We applied the proposed method, with covariates entering
linearly, to each simulated dataset, and also calculated local
confidence intervals using the method in Section~\ref{sec-local}.
Figure~\ref{fig-sims} shows the RMSE and coverage results for both
global and local estimates.
As our theoretical results predict, error in the global estimates
converges to \(0\) as the number of geographies increased; error in the
local estimates decreases but is lower-bounded by the intrinsic variance
of the local parameters. Error was little affected by the number of
covariates \(p\) or the strength of confounding (correlation with
\(\overline X\)), but did increase substantially with the number of
predictors \(d\). This indicates that in many ecological inference
applications, the main statistical challenge is that of many predictors,
not many covariates. This further supports the routine use of many
covariates. Across combinations of \(m\) and \(p\), coverage rates were
close to their nominal levels for \(d>2\), but above nominal levels for
\(d=2\). This is somewhat surprising given that the error grows with
\(d\). Coverage of the global and local confidence intervals is close to
the nominal level, though the coverage of the local intervals falls
somewhat below 95\% coverage for large \(m\) and \(d>2\).\footnote{ We
observe undercoverage as well for regression-based estimates of the
global parameter without the Riesz representer adjustment, suggesting
finite-sample estimation error in the regression may be to blame.}
\subsection{Voter file validation}\label{sec-valid-real}
The simulation study results, while encouraging, have the virtue of a
data-generating process that exactly satisfies the required assumptions
here. We therefore turn next to a much more challenging real-world
setting, where we cannot verify that the assumptions hold exactly. This
also provides an opportunity to test the sieve estimation methods for
\(\gamma\) and \(\alpha\); in the simulation studies, the true models
were linear in the covariates.
Our data consist of 1,759 precincts in the Miami metropolitan area. The
quantity of interest is the proportion of a racial group's party
registrants who register for the Republican party. The Miami area has a
mix of different racial groups, including Cuban Americans, who are
well-known to political observers as having systematically more
Republican political preferences than other Hispanic groups, so it
serves as a good test case for our method.\footnote{ For example, the
registration file reveals that 40 percent of Hispanic registrants
living in Census tracts where the majority of Hispanic voters are of
Cuban origin are Republicans, but only 24 percent of Hispanic
registrants living in other Census tracts are Republican. This
correlation between a covariate and the outcome of interest would lead
to bias unless one can properly adjust for confounding covariates.}
Data on party registration come from Florida voter registration records,
and we augment this data with Census data on the racial composition of
each precinct, along with other covariates such as Hispanic origin,
population density, income, age, and past election results.
Section~\ref{sec-valid-detail} describes the voter file data and
covariates in more detail. Crucially, in Florida, voter registration
records record both a voter's party registration and their racial
affiliation, so we observe the true value of the estimand.
\begin{figure}[t]
\centering{
\sbox\pandoc@box{\includegraphics[keepaspectratio]{figures/miami.pdf}}
\Gscale@div\@tempa{\textheight}{\dimexpr\ht\pandoc@box+\dp\pandoc@box\relax}
\Gscale@div\@tempb{\linewidth}{\wd\pandoc@box}
\ifdim\@tempb\p@<\@tempa\p@\let\@tempa\@tempb\fi
\ifdim\@tempa\p@<\p@\scalebox{\@tempa}{\usebox\pandoc@box}
\else\usebox{\pandoc@box}
\fi
}
\caption{\label{fig-miami-topline}\textbf{Accuracy of predicting
Republican registration by racial group.} \textbf{(a) Global accuracy.}
The labeled vertical lines indicate the true value of Republican
registration from the voter file in the Miami metropolitan area, and the
horizontal grey lines indicate the partial identification bounds for
each group. \textbf{(b) Local accuracy.} A scatterplot of the local
estimates versus true values of Republican registration for each racial
group for 100 randomly sampled precincts. 95\% local confidence
intervals are shown for each point.}
\end{figure}
We apply our proposed method using three different sets of covariates
for fitting \(\gamma\) and \(\alpha\). The main specification controls
for all 16 continuous covariates and dummy variables for the county and
subdivision (around 30 levels), and uses a \(J_m=1000\) tensor-product
cosine basis as described in Section~\ref{sec-bases}. A second
specification uses only one covariate, the percentage of Hispanic adults
in the Census tract that are of Cuban origin, modeled using the same
basis expansion with \(J_m=100\). Finally, we also fit our method with
no covariates and thus without penalization, which is equivalent to a
simple linear regression. The estimates from all three specifications
are displayed in the left panel of Figure~\ref{fig-miami-topline} for
the four major racial groups, along with the true values from the voter
file.
A linear regression with no covariates overestimated White GOP
registration by \(9\) percentage points (pp), overestimated Hispanic GOP
registration by \(9\)pp, and produced impossible, negative estimates for
Black and Asian voters. Controlling for covariates with the proposed
method moves all of these estimates in the correct direction. In the
more complex model with all covariates, the estimate for White voters is
only \(2\)pp off, and the estimate of Hispanic voters is only \(0.4\)pp
off. Estimates for Black voters are also no longer negative and only
\(2.5\)pp off. The estimate for Asian voters is quite variable, given
the small fraction of Asian voters in Miami, but the error is still a
double-digit improvement over the simple regression. Importantly, all
four confidence intervals for the full specification cover the true
value (just barely, for White voters).
Finally, we also evaluate the accuracy of the precinct-level local
estimates obtained using the methods in Section~\ref{sec-local}. The
right panel of Figure~\ref{fig-miami-topline} shows a scatterplot of the
local estimates versus the true values for 100 randomly sampled
precincts. For larger racial groups like White and Hispanic voters, the
local estimates are quite accurate, with an overall RMSE of \(7.3\)pp
and \(9.7\)pp, respectively. For smaller racial groups, the estimates
are shrunk towards a global mean, and the RMSE is higher: \(14\)pp for
Asian voters, for instance. Critically, however, the local confidence
intervals for all groups cover at least the nominal rate. Even using the
narrower intervals implied by a unimodality assumption, coverage of
\(95\%\) intervals (averaged across precincts) is \(96\%\), \(97.6\%\),
\(99.8\%\), and \(97.9\%\) for White, Hispanic, Black, and Asian voters,
respectively.
\section{Application}\label{sec-appl}
We now apply our method to the problem studied by \citet{jbaily2022air}:
estimating exposure to fine particulate matter (\(\text{PM}_{2.5}\)) by
racial and income groups. This application also illustrates our
sensitivity analysis.
Our data consist of 31,853 ZIP Code Tabulation Areas (ZCTA) in the U.S.
The outcome is the average \(\text{PM}_{2.5}\) exposure in 2016 for each
ZCTA, and our main predictor variable is race by income combination. The
predictors are coded as seven household income bins and two racial
groups, White and Other. Our covariates consists of ZCTA-level
population density, fraction of the over-65 population in poverty,
fraction of the over-65 population without a high-school degree, an
indicator for urbanity, as well as the latitude and longitude of the
centroid of each ZCTA. We apply the proposed estimator using a
tensor-product cosine basis on the non-geographic covariates (18 terms)
and the geographic coordinates (400 terms). All in all, the ridge
regression procedure fits 5,852 coefficients; this takes around 30
seconds on a modern laptop.
As discussed in Section~\ref{sec-est-semi}, the second moments of the
fitted Riesz representer serve as a way to assess the positivity
assumption which underlies estimation. Here, they range from \(26.1\) to
\(74.9\), which is larger than the minimum possible value of 1. This
reflects the fact that few ZCTAs comprise entirely one race-income
group: most of the \(\overline X_j\) are closer to 0 than to 1. As a result,
more extrapolation is needed to estimate the conditional mean for each
group. Exploratory analysis confirms reasonable variation in each
\(\overline X_j\), however, and so we believe that Assumption~\ref{asm-pos}
is plausible. The larger second moments will lead to more variable
estimates, however.
\begin{figure}[t]
\centering{
\sbox\pandoc@box{\includegraphics[keepaspectratio]{figures/pollution.pdf}}
\Gscale@div\@tempa{\textheight}{\dimexpr\ht\pandoc@box+\dp\pandoc@box\relax}
\Gscale@div\@tempb{\linewidth}{\wd\pandoc@box}
\ifdim\@tempb\p@<\@tempa\p@\let\@tempa\@tempb\fi
\ifdim\@tempa\p@<\p@\scalebox{\@tempa}{\usebox\pandoc@box}
\else\usebox{\pandoc@box}
\fi
}
\caption{\label{fig-air}\textbf{(a) Estimates of pollution exposure by
race and income group.} Circles are centered at the estimate and have
area proportional to the size of each group in the U.S. population.
Vertical lines display 95\% confidence intervals. \textbf{(b)
Sensitivity analysis for the racial disparity in exposure among people
earning less than \$20,000.} The two sensitivity parameters plotted
along each axis, and the contours indicate the bias in the estimated
difference (Other - White) that would arise from the specified degree of
confounding. The blue contour corresponds to bias that would be
sufficient for the estimated disparity to be zero. Benchmarked
sensitivity parameters for observed covariates are also plotted.}
\end{figure}
Figure~\ref{fig-air} (a) displays the estimates for each race and income
group along with 95\% confidence intervals. There is no clear income
disparity within either racial group, but there are clear disparities
across racial groups within the lower income categories. These
disparities are statistically significant but not particularly large:
monthly variation in \(\text{PM}_{2.5}\) exposure can be on the order of
\(10\ \mu g/m^3\) \citep{rao2011understanding}. The direction of this
disparity is consistent with the findings of \citet{jbaily2022air}.
The estimates in Figure~\ref{fig-air} (a) rely on
Assumption~\ref{asm-car} holding: that conditional on the population
density, education, poverty, urbanity, and approximate geographic
location of a ZCTA, air pollution exposure is unrelated to the racial
composition of the ZCTA. While this assumption seems plausible,
especially due to the control for geographic location, it is important
nonetheless to assess the sensitivity of the estimates to violations of
this assumption.
For simplicity, we show only the sensitivity analysis for the difference
in exposure between Other and White residents earning less than \$20,000
per year (leftmost points in Figure~\ref{fig-air} (a)). The point
estimate is that the non-White population is exposed to \(2\ \mu g/m^3\)
more pollution than the White population. We first calculate the
robustness value for bias equal to this point estimate, which is
\(0.0314\). This means that if either
\(R^2_{\overline Y\sim A\mid \overline X, Z}\) or \(1-R^2_{\alpha^A\sim \alpha}\)
is larger than this value, then the bias could be large enough to
explain away the entire estimated disparity.
To better interpret these robustness values, Figure~\ref{fig-air} (b)
shows a sensitivity contour plot. For each combination of sensitivity
parameters, the contour lines indicate the size of the bias that would
arise from confounding of that magnitude. The contour labeled
``Estimated difference'' marks the dividing line at which the disparity
estimate would change sign. This contour is rather close to the origin,
indicating substantial sensitivity, which agrees with the high
sensitivity implied by the small robustness value.
Figure~\ref{fig-air} (b) also displays benchmarked values of the
sensitivity parameters for observed covariates (detailed in
Appendices~\ref{sec-bench} and \ref{sec-pollution}). These benchmarks
show that if an omitted confounder is of similar strength to ZCTA
education, population density, urbanity, or poverty, it would likely not
change the sign of the disparity estimate.
In contrast, the location variable has a much larger benchmarking value.
The value is closer to 40, which is far larger than the estimated
difference. In other words, if the omitted confounder is of similar
strength to geographic location, then it would easily change the sign of
the estimate and create substantial bias. In a more in-depth analysis,
these findings would prompt us to consider collecting other covariates,
and more carefully evaluate the model specification as regards the
critical geography covariate.
\section{Conclusion}\label{sec-concl}
We have formalized the identification assumptions for ecological
inference and proposed a new set of tools for estimation and sensitivity
analysis. We stress for practitioners the importance of thinking
carefully about the identifying assumptions, rather than blindly
applying existing ecological inference methods without any covariates.
With the tools presented here, the plausibility of and sensitivity to
these identifying assumptions can be directly assessed, and when they
are judged reasonable, estimation may be carried out efficiently. We
strongly recommend the routine use of sensitivity analysis in performing
ecological inferences.
One drawback of the proposed estimator is that when \(Y\) is bounded,
the regression \(\widehat\gamma\) can be fit to respect these bounds, but
the overall estimate \(\widehat\beta\) may not, due to the form of the
efficient influence function (Eq.~\ref{eq-neyman}). Future work could
explore ways to modify the estimator to respect these bounds, which may
further reduce error in finite samples.
There are other possible extensions of the methods proposed here. One
interesting case is when only \(Y\) but not \(X\) is aggregated, such as
when a voter's ballot \(Y\) is secret and can only be observed in
aggregate at the precinct level, but many individual-level covariates
\(X\) are available from voter files or surveys.
\citet{flaxman2015supported} and \citet{fishman2024estimating} have
proposed methods for this setting, but results on identification, and
estimation guarantees, remain limited. Related to this case is the
challenge of estimating conditional means for a variety of \(X\) at once
(e.g., race and income and education), or for a high-dimensional \(X\).
Another possible extension is to leverage spatial correlation in the
data, which is likely to be present in many applications (such as in
Section~\ref{sec-appl}), and may allow for estimating a latent spatial
confounder. Finally, we have focused on estimating conditional means
here (and conditional variances in Section~\ref{sec-local}), but other
estimands may be of interest, such as quantiles.
\section*{References}\label{references}
\addcontentsline{toc}{section}{References}
\bibliography{references.bib}
\newpage{}