Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
84,352 characters · 20 sections · 56 citation commands
=0pt =0pt plus .5=0pt plus .5=.3Identification and Semiparametric Estimation of Conditional Means from Aggregate Data
\allsectionsfont
Keywords\quad aggregate data \textbullet ecological inference \textbullet double/debiased machine learning
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 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\)) 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 cho2008cross. This is a common challenge in epidemiology, where exposure and disease are measured separately 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 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 robinson1950ecological and Goodman \citetext{goodman1953ecological; 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 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 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 king1997solution combined the regression framework with the partial identification bounds in a Bayesian varying coefficient model, spurring numerous extensions rosen2001bayesian, wakefield2004ecological, james2009r, though the approach was somewhat controversial 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 freedman1998solution; in practice, users of ecological inference methods rarely evaluate the plausibility of their assumptions.\footnote{ An exception to the pattern is 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. km2025review discuss these computational challenges in more detail.}
To address the identification ambiguity, we formalize the ecological inference problem in Section (ref) 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) 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)).
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 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)). First, in addition to our estimation theory for the so-called global estimand \(\operatorname{\mathbb{E}}[Y\mid X]\), we develop asymptotically valid confidence intervals for the 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 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 seine, which we apply in a demonstration in Section (ref) to the air pollution data of jbaily2022air. Together with the formalization of identifying assumptions, these methods place ecological inference on a more practical and robust foundation.
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 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\).
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 accounting identity:
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 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 \[
\] 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)\).
Eq. (ref) 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) a type of random-coefficient model, albeit one with no error term. These models are well-studied 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 beran1992estimating state the assumption without covariates, it can be easily weakened to hold conditional on covariates.
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 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) holds in the data at hand: researchers should rely on their substantive knowledge. However, as we will see, Assumption (ref) 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).
Because of the weighting by \(N_g\), Assumption (ref) 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.
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) is slightly stronger than Assumption (ref), 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) but not Assumption (ref) 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, 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) or Assumption (ref) is sufficient for identification of \(\beta\). All proofs are deferred to Appendix (ref).
Note that when there are no covariates, i.e., \(Z\) is null, Assumption (ref) is strong and implausible. In the air pollution example, Assumption (ref) 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) 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) 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) 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) records, the following assumption is sufficient for Assumption (ref).
Because \(G\) appears directly in Assumption (ref), 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). There may be situations where Assumption (ref) does not hold while Assumption (ref) does.
In this section, we apply Assumption (ref) and Theorem (ref) 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):
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 closed linear subspace. First, however, we discuss estimation when Assumption (ref) holds without covariates.
The first ecological inference methods were based on simple linear regression of \(\overline Y_G\) on \(\overline X_G\) 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). 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) 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) 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 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 duncan1953alternative and king1997solution argue can be substantial. When \(d=2\), 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), and they both implicitly assume Assumption (ref) 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) 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 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.
By Theorem (ref), we can write the global estimand using the notation in Eq. (ref) as
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 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) 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.
Assumption (ref) bounds \(u\), ensuring that no single geography dominates the estimand. Assumption (ref), 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) 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 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.
We refer to \(\alpha_{0j}\) as the Riesz representer of \(\beta_j\); critically, it also belongs to the restricted class \(\Gamma\). While it is defined implicitly, in Appendix (ref) 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).
The second moment of the Riesz representer is tied directly to the modulus of continuity that establishes Proposition (ref), which in turn hinges critically on Assumption (ref). 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) not holding. This is analogous to causal inference, where the distribution of the propensity scores plays a similar role in assessing overlap.
Corollary (ref) implies that we could estimate \(\beta_j\) in two ways: either by estimating \(\gamma_0\) and then plugging into Eq. (ref), 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\) newey1994asymptotic:
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\) chernozhukov2022riesz. This Neyman orthogonality property is closely related to double robustness: in fact, if \(\gamma\) is properly specified in Eq. (ref), then even with a misspecified \(\alpha\), the score in Eq. (ref) 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
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), 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.
This differs from the random-index setup of Section (ref), 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 as if i.i.d. Our validation studies demonstrate that this is a reasonable assumption in practice.
As discussed in chernozhukov2022riesz, \(\alpha_{0j}\) can be estimated via the following representation:
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). 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 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), including interactions of polynomials in each variable, tensor-product splines, and a more recent tensor-product cosine basis 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
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). 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 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) several other practical considerations in implementing the estimator: bounded \(Y\), the use of weights in estimation, and the choice of sieve basis.
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 park2015varying. However, most treatments focus on regression estimators and do not directly apply to estimating \(\alpha_{0j}\). Moreover, Assumption (ref) and Assumption (ref) 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.
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) (3). The next theorem shows that this is possible at the same rate.
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, 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.
The proposed estimator \(\widehat\beta_{mj}\) in Eq. (ref) is a double/debiased machine learning (DML) estimator 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 chen2022debiased, which rely on an algorithmic stability condition that is satisfied by the ridge penalty used here.
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) 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.
In this section, we discuss two extensions of the proposed method: estimation of the local quantities \(B_g\) for each geography \(g\), and a sensitivity analysis for violations of Assumption (ref). A third extensions, a possible hypothesis test for Assumption (ref), is presented in Appendix (ref).
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), \(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)), 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.
Assumption (ref) could be equivalently written in terms of the residuals \(\varepsilon_G\). Of course, both Assumption (ref) and Assumption (ref) 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) that does not condition on \(N_G\) could also be applied, analogously to Assumption (ref).
Denote the covariance matrix as \(\Sigma\) and let \(\Sigma(z):=\operatorname{\mathbb{V}}[\varepsilon_G\mid Z_G=z]\). Under Assumption (ref), \(\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.
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) 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) using an estimated \(\widehat\kappa\).\footnote{ In finite samples, the \(\widehat\Sigma(z)\) yielded by Proposition (ref) 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 \[
\] where \(A^+\) denotes the Moore-Penrose pseudoinverse of \(A\). Then we have the following result.
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\) vysochanskij1980justification.
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.
Every result so far has relied critically on the Assumption (ref) assumption. In practice, it is unlikely that Assumption (ref) holds 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 chernozhukov2022sens, who develop a nonparametric sensitivity analysis for estimands for which a Riesz representer exists.
Rather than assume that Assumption (ref) 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) 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) 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\). chernozhukov2022sens (Theorem 2 and Corollary 2) then establish the following result.
In other words, the bias due to violations of Assumption (ref) 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) 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) to each group separately and then using the triangle inequality.
Researchers can vary the sensitivity parameters \(C_\gamma\) and \(C_\alpha\) to understand how sensitive their estimates are to violations of Assumption (ref). 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 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, 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) under slightly modified regularity conditions. The full conditions are stated in chernozhukov2022sens, who also propose DML confidence bounds for these bounds which involve further computation.
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). Appendix (ref) 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 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) and demonstrated in the application, is used in causal inference as well.
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.
We examine the performance of our method on data simulated from the data generating process assumed by the now-standard method of 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).
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 goodman1953ecological, (3) the truncated-normal model of king1997solution, from the R package ei, both with and without covariates, and (4) the Multinomial-Dirichlet count model of rosen2001bayesian, implemented as ei.MD.Bayes in the R package eiPack, both with and without covariates.
The proposed method achieved the lowest root mean square error (RMSE) in estimating the global parameters \(\beta\), with King's king1997solution model with covariates a close second. Table (ref) 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 rosen2001bayesian, which is the only method implemented in public software that can handle \(d>2\), suffers higher error and 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.
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). Figure (ref) 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.}
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) 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.
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). 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) 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). The right panel of Figure (ref) 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.
We now apply our method to the problem studied by 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), 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) is plausible. The larger second moments will lead to more variable estimates, however.
Figure (ref) (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\) rao2011understanding. The direction of this disparity is consistent with the findings of jbaily2022air.
The estimates in Figure (ref) (a) rely on Assumption (ref) 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) (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) (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) (b) also displays benchmarked values of the sensitivity parameters for observed covariates (detailed in Appendices (ref) and (ref)). 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.
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)). 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. flaxman2015supported and 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)), and may allow for estimating a latent spatial confounder. Finally, we have focused on estimating conditional means here (and conditional variances in Section (ref)), but other estimands may be of interest, such as quantiles.
\addcontentsline{toc}{section}{References}