EconBase
← Back to paper

Partial Identification of Individual-Level Parameters Using Aggregate Data in a Nonparametric Model

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.

86,983 characters · 11 sections · 39 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Partial Identification of Individual-Level Parameters Using Aggregate Data in a Nonparametric Model10pt

titlepage\onehalfspacing \abstract{ \selectfont I develop a methodology to partially identify linear combinations of conditional mean outcomes when the researcher only has access to aggregate data. Unlike the existing literature, I only allow for marginal, not joint, distributions of covariates in my model of aggregate data. Bounds are obtained by solving an optimization program and can easily accommodate additional polyhedral shape restrictions. I provide a procedure to construct confidence intervals on the identified set and demonstrate performance of my method in a simulation study. In an empirical illustration of the method using Rhode Island standardized exam data, I find that conditional pass rates vary across student subgroups and across counties. \\[5pt] Keywords: Aggregate data, partial identification, ecological inference, nonparametric \\ JEL Codes: C14 \\ }

\onehalfspacing

\setcounter{page}{2}

Introduction

It has been long known that the relationship between variables at the individual level can be different from the relationship between those same variables aggregated over individuals. The ecological inference literature, which presents potential solutions to this problem, dates back to the seminal work of robinson1950ecological, duncan1953ecological, and theil1954linear. Much of the literature is concerned with point identification of individual-level parameters, which requires strong assumptions that may be implausible in applied settings (see, e.g., cho1998iff, cho2004limits, freedman1998solution and kousser2001ecological). Throughout this paper I use the phrase “individual-level parameters” to refer to parameters that rely on the joint distribution of the individual-level variables.

Absent such strong assumptions, I can still partially identify individual-level parameters based on aggregate data. The resulting identified set precisely quantifies the extent to which individual-level results are sensitive to assumptions. In this paper I build upon the methods of cross2002regressions and cho2008cross to derive partial identification of individual-level parameters when only aggregate data is available. I consider a model of aggregate data where marginal distributions of covariates are observed together with a marginal average outcome across different groups. The individual-level parameters of interest are linear combinations of conditional mean outcomes $\mathbb{E}[Y_i|X_{1i},\dots, X_{Li}]$, which encompass parameters like average predictive effects. I construct bounds by solving an optimization problem that considers all joint distributions of individual-level variables that are consistent with the observed marginal information. The optimization problem formulation allows for easy accommodation of additional model restrictions.

I provide a consistent plug-in estimator for the bounds and a valid, albeit conservative, inference procedure that works with aggregate data. I conduct a simulation study to demonstrate the performance of the estimation and inference procedures and what features of the aggregate data drive the width of the bounds. I find that when the aggregate data provides stronger restrictions on the joint distribution of covariates, bounds are more informative.

To examine the informativeness of these bounds in practice, I apply this methodology to a Rhode Island standardized exam dataset. I find that bounds are very wide on conditional pass rate gaps of interest. Imposing monotonicity shape restrictions narrows the bounds. With restrictions implied by additional subgroup pass rate data, bounds are more informative, especially under monotonicity and subgroup pass rate restrictions. A homogeneity restriction on the underlying individual-level model results in empty estimated identified sets.

Relevant to this paper is the literature on partial identification when combining multiple data sets, also known as ecological inference, which deals with similar issues of needing to infer joint information that is unobserved cross2002regressions, molinari2006generalization, ridder2007econometrics, fan2014identifying, fan2016estimation, buchinsky2022estimation, haultfoeuille2024partial, elzayn2025monotone, mccartan2025identification. Most relevant is cho2008cross, who use the methods of cross2002regressions to derive bounds on conditional mean outcomes when the observed data is the marginal distribution of $Y_i$ over groups together with the joint distribution of covariates $X_i = (X_{1i},\dots, X_{Li})$ for many groups. In contrast, in my model of aggregate data the researcher observes the marginal, not the joint, distribution of the covariates $X_{1i},\dots, X_{Li}$ for groups. My model better fits the aggregate data seen in practice. For example, aggregate data can provide ethnicity distributions, gender distributions, education distributions, and income distributions but often does not provide the joint distribution of all of these variables due to privacy concerns. I argue that the existing cross2002regressions and cho2008cross methods do not immediately apply for the kind of aggregate data I consider. Furthermore, I demonstrate that approaches that are valid with aggregate data and directly use cross2002regressions bounds produce bounds that are not sharp.

The rest of the paper proceeds as follows. Section (ref) presents identified sets on the parameters of interest. Section (ref) presents a discussion of consistent estimation and inference procedures. In Section (ref) I perform three simulation studies to evaluate the performance of the inference procedure and explore what features of the data drive bound width. Section (ref) presents an empirical illustration of the methodology. Section (ref) concludes.

Identification

In this section I construct identified sets for linear combinations of conditional mean outcomes using aggregate data. I model aggregate data as marginal distributions of variables for many different groups. Suppose there exists a sequence of latent random variables $(Y_i, X_{1i}, \dots, X_{Li}, G_i), i = 1, \dots, n$, where $Y_i$ with support $[y_\ell, y_u]$ denotes individual $i$'s outcome, $X_i \equiv (X_{1i}, \dots, X_{Li})'$ with support $\{x_k\}_{k=1}^K$ denotes individual $i$'s covariates, and $G_i$ with support $\{1,\dots, G\}$ denotes individual $i$'s group. Aggregate data identifies the average outcome within each group $\mathbb{E}[Y_i|G_i = g]$, the marginal distributions of each of the $L$ covariates within each group $\mathbb{P}[X_{\ell i} = x_{k,\ell}|G_i = g]$, and the relative size of each group $\mathbb{P}[G_i = g]$. In practice aggregate data consists of the sample equivalents of these values. I deal with sampling uncertainty in Section (ref).

I thus impose the following assumptions on the random variables:

assumption\begin{enumerate}[label=(\roman*)] • $Y_i$ is a random variable with bounded support $[y_\ell, y_u]$.\footnote{The assumption of bounded support of $Y_i$ is for tractability of computing the identified set. Similar to cross2002regressions, the main results in this section will still hold with unbounded support.} • $G_i$ is a discrete random variable with finite support $\{1,\dots,G\}$. • $X_i$ is an $L$-dimensional discrete random vector with finite support $\{x_k\}_{k=1}^K$ for $K \geq 2$.\footnote{The case of continuous covariates is beyond the scope of this paper. The assumption $K \geq 2$ ensures the problem is not trivial.} • The joint distribution of the random variables $(Y_i, G_i, X_i)$ is latent. Instead, the researcher observes $\mathbb{E}[Y_i|G_i = g], \mathbb{P}[X_{\ell i} = x_{k,\ell}|G_i = g],$ and $\mathbb{P}[G_i = g]$ for every $\ell = 1, \dots, L, k = 1, \dots, K, $ and $g = 1, \dots, G$. Furthermore, $n$ is observed. \end{enumerate}

As an example, consider a dataset of standardized exam results and demographics. $G_i$ denotes student $i$'s school district, $Y_i$ is an indicator for whether student $i$ passed the exam or not, and $X_i$ are student $i$'s demographics, like race and socioeconomic status. The researcher observes (sample estimates of) the pass rate for every school district $\mathbb{E}[Y_i|G_i=g]$, marginal demographics of every school district $\mathbb{P}[X_{\ell i} = x_{k,\ell}|G_i = g]$, and the number of students enrolled in each school district, with which one can obtain $\mathbb{P}[G_i = g]$. Alternatively, $Y_i$ could be student $i$'s exam score, and the researcher observes (sample estimates of) the average exam score for every school district $\mathbb{E}[Y_i|G_i=g]$.

The goal is to construct bounds on linear combinations of $\mathbb{E}[Y_i|X_i = x_k]$ for given weights $\{\lambda_k\}_{k=1}^K$, $\sum_{k=1}^K \lambda_k \mathbb{E}[Y_i|X_i = x_k]$. For example, if the researcher is interested in the average predictive effect of changing $X_i$ from $x_{k_1}$ to $x_{k_2}$ on $Y_i$, the researcher can choose $\lambda_{k_2} = 1, \lambda_{k_1} = -1,$ and $\lambda_k = 0$ for all other $k$. Note that this expectation is taken over the population of individuals $i$; the number of groups $G$ is fixed. I will construct identified sets using only what is observed in aggregate data.

Construction of the identified set relies on bounding joint distributions of variables using marginal distributions. The object of interest, linear combinations of $\mathbb{E}[Y_i|X_i = x_k]$, depends both on the joint distribution of covariates $X_i$ and on the joint distribution of outcome $Y_i$ with covariates $X_i$. Neither of these distributions are observed in aggregate data. Instead, I construct bounds on these joint distributions similar to cross2002regressions, and discuss how they relate to the cross2002regressions bounds in Remark (ref). I first use the marginal distributions of the $L$ covariates $\mathbb{P}[X_{\ell i} = x_{k,\ell}|G_i = g]$ to obtain bounds on the joint distribution of the covariates $\mathbb{P}[X_i = x_k|G_i = g]$. I then use the average of the outcome $\mathbb{E}[Y_i|G_i=g]$ together with bounds on the joint distribution of covariates to obtain bounds on the conditional mean group outcomes $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$. Finally, I aggregate over groups to obtain bounds on (linear combinations of) the conditional mean outcomes $\mathbb{E}[Y_i|X_i = x_k]$.

I can bound joint distributions using marginal distributions by applying the law of iterated expectations and the law of total probability. The following three equations must hold:

align[align omitted — 465 chars of source]

From Bayes' rule and the law of total probability I also know

align[align omitted — 263 chars of source]

I can use equations (ref) and (ref) to rewrite (ref) as

align*[align* omitted — 214 chars of source]

which can be rearranged as

align[align omitted — 236 chars of source]

This rearrangement is equivalent when $\mathbb{P}[X_i = x_k] = \sum_{g=1}^G \mathbb{P}[G_i = g]\mathbb{P}[X_i = x_k|G_i = g] > 0$. If $\mathbb{P}[X_i = x_k] = 0$ then $\mathbb{E}[Y_i|X_i = x_k]$ is undefined. If $\mathbb{E}[Y_i|X_i = x_k]$ were to exist then the only information available is that it lies in $[y_\ell, y_u]$ by Assumption (ref).i. Thus I will work with equation (ref), which is still well-defined and valid when $\mathbb{P}[X_i = x_k] = 0$ because both sides are zero when $\mathbb{E}[Y_i|X_i = x_k]$ is finite.

To make clear what is known and what is unknown under Assumption (ref), I will replace all unknown objects with variables in equations (ref), (ref), and (ref). Let $p_{kg}$ denote $\mathbb{P}[X_i = x_k|G_i = g]$, $c_{kg}$ denote $\mathbb{E}[Y_i|X_i = x_k,G_i = g]$, and $d_{k}$ denote $\mathbb{E}[Y_i|X_i = x_k]$. Then I have the following three equations:

align[align omitted — 345 chars of source]

In addition to these equations, I also know that all probabilities $p_{kg}$ are non-negative and that for each $g$, $p_{1g}, \dots, p_{Kg}$ sum to 1. By Assumption (ref).i, any conditional expectation of $Y_i$ must be between $y_\ell$ and $y_u$. As in cross2002regressions, these are the only relationships I can use to relate joint information of interest to the observed marginal information in the data without any further assumptions.

I first characterize the identified set $P$ for the probability distribution of $X_i|G_i$ under Assumption (ref) from equations (ref) and (ref) and the additional restrictions discussed above:

multline[multline omitted — 344 chars of source]

Then from set $P$ and equations (ref) and (ref), I can characterize the identified set $D$ for $\sum_{k=1}^K \lambda_k \mathbb{E}[Y_i|X_i = x_k]$, the parameter of interest, under Assumption (ref):

multline[multline omitted — 433 chars of source]

From this characterization it is not immediately clear how to compute the identified set. For this it is useful to re-characterize the identified set as a constrained programming problem, similar to honore2006bounds. In the following proposition I show that $D$ can be written as an interval $[L,U]$, where $L$ and $U$ are the solutions to bilevel optimization problems over $\{d_k\}$, $\{c_{kg}\}$, and $\{p_{kg}\}$. Furthermore I formally state and prove that $D$ is sharp for parameter $\sum_{k=1}^K \lambda_k \mathbb{E}[Y_i|X_i=x_k]$.

propositionUnder Assumption (ref), $D$ is the sharp identified set for parameter $\sum_{k=1}^K \lambda_k \mathbb{E}[Y_i|X_i=x_k]$. In addition, $D$ is given by $ D = [L,U]$, where \begin{multline} L = \inf_{\{p_{kg}\},\{c_{kg}\} , \{d_k\}} \sum_{k=1}^K \lambda_k d_k s.t. \{p_{kg}\}_{k=1,g=1}^{K,G} \in P, \{c_{kg}\}_{k=1,g=1}^{K,G} \in [y_\ell,y_u]^{KG}, \{d_k\}_{k=1}^{K} \in [y_\ell, y_u]^K, \\ d_k \sum_{g=1}^G \mathbb{P}[G_i = g]p_{kg} = \sum_{g=1}^G \mathbb{P}[G_i = g]p_{kg}c_{kg} \forall k, and \mathbb{E}[Y_i|G_i = g] = \sum_{k=1}^K c_{kg} p_{kg} \forall g, \end{multline} \begin{multline} U = \sup_{\{p_{kg}\},\{c_{kg}\}, \{d_k\}} \sum_{k=1}^K \lambda_k d_k s.t. \{p_{kg}\}_{k=1,g=1}^{K,G} \in P, \{c_{kg}\}_{k=1,g=1}^{K,G} \in [y_\ell,y_u]^{KG}, \{d_k\}_{k=1}^{K} \in [y_\ell, y_u]^K, \\ d_k \sum_{g=1}^G \mathbb{P}[G_i = g]p_{kg} = \sum_{g=1}^G \mathbb{P}[G_i = g]p_{kg}c_{kg} \forall k, and \mathbb{E}[Y_i|G_i = g] = \sum_{k=1}^K c_{kg} p_{kg} \forall g, \end{multline} and \begin{multline} P = \underset{\{p_{kg}\}}{arg\min} \sum_{g=1}^G \sum_{r=1}^{LK} v_{rg}^+ + v_{rg}^- \;\; s.t. v_{rg}^+, v_{rg}^- \geq 0 \forall r,g, p_{kg} \geq 0 \forall k,g, \sum_{k=1}^K p_{kg} = 1 \forall g, \\ \textnormal{and } \mathbb{P}[X_{\ell i}=x_{k,\ell}| G_i = g] - \sum_{j=1}^K \mathbbm{1}\{x_{j,\ell} = x_{k,\ell}\} p_{jg} = v_{K(\ell-1)+k, g}^+ - v_{K(\ell-1)+k, g}^- \forall \ell, k, g. \end{multline}

The proof of this proposition consists in showing that the characterization (ref) is equivalent to that of (ref) because $P$ is nonempty, and then showing that $D$ as defined in (ref) is an interval. I defer all proof details to Appendix (ref).

While the optimization problems of (ref) and (ref) have nonconvex constraints, this bilevel formulation is helpful because it suggests how computation of the lower and upper bound might be performed. In particular, the problems of (ref) and (ref) given any particular $p = \{p_{kg}\} \in P$ is a linear program in $\{c_{kg}\}$ and $\{d_k\}$. Letting $L(p)$ and $U(p)$ denote the solutions to (ref) and (ref) given $p \in P$ respectively, solving for $L(p)$ and $U(p)$ for a given $p \in P$ is a linear program and thus fast. Then I can solve for $D$ by searching for the infimum and supremum of $L(p)$ and $U(p)$, respectively, over all $p \in P$. I recommend solving for $\inf_{p \in P} L(p)$ and $\sup_{p \in P} U(p)$ using a nonconvex solver with a coarse grid of different starting points.\footnote{This may be computationally intensive when $X_i$ or $G_i$ has large support.} Checking if a candidate starting point $p$ is in $P$ is fast because the optimization problem of (ref) is a linear program.

remarkSharp bounds on $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$ jointly over the support of $X_i$ can be derived under Assumption (ref) using the method of Section 2.2 of cross2002regressions, as discussed in Appendix (ref). However, I cannot use these bounds on $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$ to directly obtain sharp bounds on (linear combinations of) $\mathbb{E}[Y_i|X_i = x_k]$ in my setting. An explanation for why is that, as can be seen in (ref), $\mathbb{E}[Y_i|X_i = x_k]$ depends on both $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$ and $\mathbb{P}[G_i = g|X_i = x_k]$. Both of these objects depend on the distribution of $X_i|G_i$, which is partially identified as set $P$ in my setting. Following (ref), one might think to obtain bounds on $\mathbb{E}[Y_i|X_i = x_k]$ from the set of all possible products of $\mathbb{P}[G_i = g|X_i = x_k]$ and $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$, summed over $g$. However, the distribution of $X_i|G_i = g$ in $P$ that minimizes (maximizes) $\mathbb{P}[G_i = g|X_i = x_k]$ may not be the same distribution that minimizes (maximizes) $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$ subject to (ref). Thus the product of the lower (upper) bounds of the sharp identified sets for $\mathbb{P}[G_i = g|X_i = x_k]$ and $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$ summed over $g$ may be strictly smaller (bigger) than the smallest (largest) possible value that $\mathbb{E}[Y_i|X_i = x_k]$ can take on. So this product approach does not produce sharp bounds. The programs of (ref) and (ref) do not have this issue because they directly minimize and maximize the parameter of interest. If I observe the joint distribution of $X_i|G_i$, as in cross2002regressions, instead of each of the marginal distributions, then $\mathbb{P}[G_i = g|X_i = x_k]$ is point identified, as can be seen in (ref). Then I can use the sharp bounds on $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$ derived in Appendix (ref) from the cross2002regressions method to simply construct bounds on linear combinations of $\mathbb{E}[Y_i|X_i = x_k]$ following (ref). Note that if there is only a single covariate in the data, $L=1$, then I trivially observe the joint distribution of $X_i|G_i$ and thus $\mathbb{P}[G_i = g|X_i = x_k]$ is point-identified, as can be seen from equations (ref) and (ref). This means that in the special case that $L=1$, my proposed bounds are equal to bounds on $\mathbb{E}[Y_i|X_i = x_k]$ obtained using the cross2002regressions method. Instead of constructing an identified set for a linear combination of $\mathbb{E}[Y_i|X_i = x_k]$, I could also construct a (sharp) identified set for a parameter like the one considered in cross2002regressions, namely any linear combination of $\mathbb{E}[Y_i|X_i = x_k, G_i = g]$, using my constrained optimization approach. In particular, using the same set $P$ in Proposition (ref) I can define a new linear program that minimizes/maximizes a linear combination of $\{c_{kg}\}_{k=1,g=1}^{K,G}$ subject to the constraints that $\{c_{kg}\}_{k=1,g=1}^{K,G} \in [y_\ell, y_u]^{KG}, \{p_{kg}\}_{k=1,g=1}^{K,G} \in P$, and $\mathbb{E}[Y_i|G_i = g] = \sum_{k=1}^K c_{kg}p_{kg}$ for all $g$. This procedure is a bilevel constrained optimization problem where both levels are linear, for which solution methodologies have been widely studied in the optimization theory literature. The nonconvexity of the first level of the bilevel optimization program of Proposition (ref) is because I target a different, albeit related, parameter.
remarkAs noted by obradovic2024identification, the formulation of the identified set as an optimization problem is convenient because additional assumptions on the data generating process can easily be accommodated by adding additional restrictions to the optimization problems, especially if the restrictions are linear or inequality restrictions. Depending on the restriction, the sharp identified set may not be an interval, but the interval defined by the lower and upper bound of the optimization problem will be the smallest closed interval containing the sharp identified set. One such restriction of interest is a polyhedral shape restriction on the conditional expectation function over groups. This can be expressed as $S_g c_{g} \leq a_g$ for each $g = 1, \dots, G$, where $c_{g} \equiv \left(c_{1g}, \dots, c_{Kg} \right)'$, $S_g \in \mathbb{R}^{s_g \times K}$ are known fixed matrices, and $a_g$ are known fixed vectors. Another restriction of interest is a homogeneity restriction on the conditional average outcome over groups, which can be expressed as $c_{kg} = c_{kg'}$ for all $k$ and $g,g'$. It may also be the case that additional data at a finer level of aggregation is available to the researcher. For example, the researcher may also observe average group outcomes conditional on one of the $L$ covariates, $\mathbb{E}[Y_i|X_{\ell i} = x_{k,\ell},G_i=g]$. This additional data implies the restriction \begin{align*} \mathbb{P}[X_{\ell i} = x_{k,\ell}|G_i=g]\mathbb{E}[Y_i|X_{\ell i} = x_{k,\ell},G_i=g] &= \sum_{j=1}^K \mathbbm{1}\{x_{j,\ell} = x_{k,\ell}\}c_{jg}p_{jg}. \end{align*} I will consider the extent to which imposing these restrictions reduces the width of the identified set in the simulations and empirical illustration in Sections (ref) and (ref). Note that if imposing a restriction makes the identified set empty, that restriction is inconsistent with the observed data.
remarkIn this paper I consider parameters that are linear combinations of conditional mean outcomes $\mathbb{E}[Y_i|X_i]$. In principle I could target other different but related parameters that can be written in a format similar to equations (ref), (ref), and (ref) using my proposed bilevel constrained optimization approach. For example, if I observed objects like $\mathbb{P}[Y_i \leq y|G_i = g]$ instead of $\mathbb{E}[Y_i|G_i = g]$ for some fixed $y$, then I could partially identify (linear combinations of) $\mathbb{P}[Y_i \leq y|X_i = x_k]$ by taking $c_{kg}$ to be $\mathbb{P}[Y_i \leq y|X_i = x_k, G_i = g]$ and using the law of total probability statements corresponding to equations (ref), and (ref). Such a parameter would speak to conditional quantiles of outcomes, instead of conditional mean outcomes. Another extension of my proposed method is to allow the linear combination weights $\lambda_k$ to possibly be data-dependent. If the data-dependent weights are point-identified, identification follows exactly as before and estimation can proceed using a plug-in approach similar to Section (ref) below. However, one would need to take into account the statistical uncertainty of these data-dependent weights when performing inference using the inference procedure of Section (ref). I do not pursue detailed exploration of these extensions in this paper and leave them to future work.

Estimation and Inference

In this section I propose a consistent estimation method for the identified set and a method for valid inference. Throughout this section I will condition on group variable $G_i$ and take it to be fixed. Thus $\mathbb{P}[G_i=g]$ is known without any sampling uncertainty by the researcher. Under this conditioning, the uncertainty captured by my inference method captures sampling uncertainty over individual outcomes and covariates, but not individual group assignments.

Estimation

In practice the researcher observes sample analogs of the population values $\mathbb{E}[Y_i|G_i = g]$ and $\mathbb{P}[X_{\ell i} = x_{k,\ell} | G_i = g]$ in the aggregate data. For all $\ell = 1, \dots, L, j = 1, \dots, K, g = 1, \dots, G$, denote the observed sample analogs as

align[align omitted — 310 chars of source]

I maintain the following assumption in order to apply a law of large numbers.\footnote{The assumption that the variables are i.i.d. can be relaxed as long as a law of large numbers still holds.}

assumption\begin{enumerate}[label=(\roman*)] • $(Y_i, X_i): i \in \mathcal{N}_g$ are i.i.d. for each $g$, where $\mathcal{N}_g \equiv \{i \in \{1, \dots, n\}: G_i = g\}$. • $\mathbb{P}[G_i = g] > 0$ for all $g = 1,\dots,G$. \end{enumerate}

Under Assumption (ref), $\bar{Y}_g$ and $P_n[X_{\ell i}=x_{j,\ell}|G_i = g]$ both converge in probability to their respective population values as $n \to \infty$ by the law of large numbers and continuous mapping theorem, since group probabilities $\mathbb{P}[G_i = g]$ are positive. Thus I can construct a plug-in estimator, denoted $\hat D_n$, for the identified set by replacing all population values in the optimization problems of Proposition (ref) with their sample estimates. The following proposition shows that the lower and upper bounds of the plug-in estimated set $\hat D_n$ are consistent.

propositionSuppose Assumptions (ref) and (ref) hold. Define $\hat D_n$ with respect to Proposition (ref). Then the lower and upper bounds of $\hat D_n$ converge in probability to the lower and upper bounds of $D$ as $n \to \infty$.

The proof relies on Berge's maximum theorem, which requires that the set of parameters satisfying the restrictions of the optimization problem of Proposition (ref) is compact-valued and continuous in the sample data. Note that both of the additional restrictions discussed in Remark (ref) are either linear restrictions or weak inequality restrictions on $\mathbb{E}[Y_i|X_i,G_i]$. Thus a simple extension of the proof shows that under either of the two additional restrictions, as long as the identified set is nonempty the lower and upper bounds of the estimated set are consistent.

Inference

Existing inference methods for estimated partially identified sets usually require knowledge of the joint distribution of the individual-level data to estimate a covariance matrix used in constructing critical values or test statistics for valid coverage. Examples of such methods include horowitz2000nonparametric, imbens2004confidence, and hsieh2022inference, in addition to standard delta method or bootstrap approaches. However, in my setting I only observe marginal distributions of each variable, so I must consider inference methods that require only marginal information of each variable.

I choose to construct marginal confidence intervals on each sample observation and use the Bonferroni correction to make the intervals jointly valid. Then a valid confidence interval for the identified set is the union of the bounds from the optimization programs of Proposition (ref) solved by plugging in each combination of sample observations lying within the jointly valid marginal confidence intervals. While this method produces confidence intervals that are conservative due to the Bonferroni correction, it is not immediately clear how to obtain confidence intervals that are meaningfully tighter. This approach has the advantage that constructing Bonferroni-corrected marginal confidence intervals is computationally simple, as described below. Thus to the extent that solving the programs of Proposition (ref) is computationally feasible, constructing confidence intervals with this method is as well.

Sample statistics of the form $P_n[X_{\ell i} = x_{j,\ell} | G_i = g]$ are sample averages of a binary random variable, as can be seen from equation (ref). Thus I can construct Clopper-Pearson marginal confidence intervals, which are finite-sample valid, for each of these statistics given that I observe $n$. Note that any other binomial proportion confidence interval that obtains asymptotically nominal coverage will provide asymptotically valid coverage; Clopper-Pearson has the advantage of being finite-sample valid, even if conservative because of the finite-sample validity.

For statistics of the form $\bar{Y}_g$, to construct construct confidence intervals on the population quantities I require standard errors for these group sample means. In the case that $Y_i$ is binary, I can again construct Clopper-Pearson intervals (or any other asymptotically valid interval). If $Y_i$ is not binary, I can use the inequality from bhatia2000better to bound the variance, and hence the standard errors, of the mean of each $Y_i|G_i = g$. The Bhatia-Davis inequality states that for a random variable $X$ with mean $\mu$, variance $\sigma^2$, and support with $x_\ell = \min(\textnormal{supp}(X)), x_u = \max(\textnormal{supp}(X))$,

equation*[equation* omitted — 85 chars of source]

Thus I can estimate an upper bound for the standard error of $\bar{Y}_g$ given that I know the sample means, the support of $Y_i$, and the number of observations in each group $g$. An asymptotically valid marginal confidence interval can be constructed with the typical Gaussian limiting distribution approach, but using this upper bound instead of the standard error. Alternatively, if in the aggregate data I directly observe the standard errors for each $\bar{Y}_g$, I can construct shorter marginal confidence intervals by using the actual standard errors.

Let $M$ be the total number of statistics in the aggregate data, that is, the total number of $\bar{Y}_g$ and $P_n[X_{\ell i} = x_{k,\ell}|G_i = g]$ statistics in the data across all groups $g$, support points $k$ and covariates $\ell$.

The inference procedure is as follows:

enumerate• For every sample statistic $\hat{p}$ construct two-sided level $1-\frac{\alpha}{M}$ asymptotically valid CIs that contain $\hat p$, denoted $\left[\hat p_L, \hat p_U \right]$, as discussed above. The resulting confidence intervals are \begin{itemize} • $\left[P_n[X_{\ell i}=x_{k,\ell}| G_i = g]_L, P_n[X_{\ell i}=x_{k,\ell}| G_i = g]_U \right]$ for each $P_n[X_{\ell i}=x_{k,\ell}| G_i = g]$, • $\left[\bar{Y}_{g,L}, \bar{Y}_{g,U} \right]$ for each $\bar{Y}_g$. \end{itemize} • Solve the optimization programs of $\hat D_n$ for all values of sample statistics within the marginal confidence intervals constructed in step 1: \begin{multline*} a) \quad \hat{P}_{CI,n} \equiv \underset{\{p_{kg}\}}{arg\min} \sum_{g=1}^G \sum_{r=1}^{2LK} v_{rg}^+ + v_{rg}^- \;\; s.t. p_{kg}, v_{rg}^+, v_{rg}^- \geq 0 \forall k, r, g, \; \sum_{k=1}^K p_{kg} = 1 \forall g, \\ P_n[X_{\ell i}=x_{k,\ell}| G_i = g]_L - \sum_{k=1}^K \mathbbm{1}\{x_{j,\ell} = x_{k,\ell}\} p_{jg} \leq \; v_{K(\ell-1)+k}^+ - v_{K(\ell-1)+k}^- \forall \ell, k, g, and \\ P_n[X_{\ell i}=x_{k,\ell}| G_i = g]_U - \sum_{k=1}^K \mathbbm{1}\{x_{j,\ell} = x_{k,\ell}\} p_{jg} \geq v_{K(L+\ell-1)+k}^+ - v_{K(L+\ell-1)+k}^- \forall \ell, k, g. \end{multline*} \begin{multline*} b) \quad \hat L_{CI,n} \equiv \inf_{\{p_{kg}\},\{c_{kg}\}, \{d_k\}} \sum_{k=1}^K \lambda_k d_k s.t. \{p_{kg}\}_{k=1,g=1}^{K,G} \in \hat P_{CI,n}, \{c_{kg}\}_{k=1,g=1}^{K,G} \in [y_\ell,y_u]^{KG}, \\ \{d_k\}_{k=1}^{K} \in [y_\ell, y_u]^K, \; \bar{Y}_{g,L} \leq \sum_{k=1}^K c_{kg}p_{kg} \forall g, \; \bar{Y}_{g,U} \geq \sum_{k=1}^K c_{kg}p_{kg} \forall g, \\ and d_k \sum_{g=1}^G \mathbb{P}[G_i = g] p_{kg} = \sum_{g=1}^G \mathbb{P}[G_i = g] p_{kg}c_{kg} \forall k. \end{multline*} \begin{multline*} c) \quad \hat U_{CI,n} \equiv \sup_{\{p_{kg}\},\{c_{kg}\}, \{d_k\}} \sum_{k=1}^K \lambda_k d_k s.t. \{p_{kg}\}_{k=1,g=1}^{K,G} \in \hat P_{CI,n}, \{c_{kg}\}_{k=1,g=1}^{K,G} \in [y_\ell,y_u]^{KG}, \\ \{d_k\}_{k=1}^{K} \in [y_\ell, y_u]^K, \; \bar{Y}_{g,L} \leq \sum_{k=1}^K c_{kg}p_{kg} \forall g, \; \bar{Y}_{g,U} \geq \sum_{k=1}^K c_{kg}p_{kg} \forall g, \\ \textnormal{and } d_k \sum_{g=1}^G \mathbb{P}[G_i = g] p_{kg} = \sum_{g=1}^G \mathbb{P}[G_i = g] p_{kg}c_{kg} \forall k. \end{multline*} • The confidence interval is given by $\hat D_{CI,n} \equiv [\hat L_{CI,n}, \hat U_{CI,n}]$.
propositionSuppose Assumption (ref) holds. Define $\hat D_{CI,n}$ as in the inference procedure described above. Then $\displaystyle \lim_{n \to \infty} \mathbb{P}[D \subseteq \hat D_{CI,n}] \geq 1-\alpha$.

Proposition (ref) says that the confidence interval $\hat D_{CI,n}$ has correct asymptotic coverage for the identified set. Since $\sum_{k=1}^K \lambda_k \mathbb{E}[Y_i|X_i = x_k] \in D$, $ \lim_{n \to \infty} \mathbb{P}\left[\sum_{k=1}^K \lambda_k \mathbb{E} [Y_i|X_i = x_k] \in \hat D_{CI,n} \right] \geq 1-\alpha$, and thus the confidence interval also has correct asymptotic coverage for the target parameter, as discussed in imbens2004confidence. Note that if outcome $Y_i$ is binary and one uses Clopper-Pearson confidence intervals in step 1 for all statistics, coverage in the above proposition is finite-sample valid.

Simulations

In this section I present results from three simulation studies in order to illustrate the coverage properties of the inference procedure of Section (ref) and what features of the aggregate data drive the width of the identified set. I also compare the sharp bounds obtained using my method to bounds obtained using the product approach with the cross2002regressions method in Remark (ref).

I consider a binary outcome $Y_i \in \{0,1\}$, three binary covariates $X_i = (X_{1i}, X_{2i}, X_{3i}) \in \{0,1\}^3$, and five groups $G_i \in \{1, \dots, 5\}$. Across all three simulation studies I assume the following model for how latent individual outcome $Y_i$ relates to latent covariates $X_i$:

align*[align* omitted — 121 chars of source]

I present the true population aggregate features of the data for the first simulation study in Table (ref), for the second simulation study in Table (ref), and for the third simulation study in Table (ref). In each study $\mathbb{E}[Y_i|G_i]$ is implied by the choice of marginal distributions for covariates $X_i$ and the above model. The data-generating process of the first study is calibrated to approximate the data used in the empirical illustration. In the second study, I maintain the same marginal distribution for $X_{3i}$ but push the distributions for $X_{1i}$ and $X_{2i}$ much closer to 1 and 0 respectively. In the third study I maintain the same marginal distribution for $X_{2i}$ and $X_{3i}$ as in the second study, but instead enforce that some groups have almost all individuals with $X_{1i} = 0$ and other groups have almost all individuals with $X_{1i} = 1$.

table[table omitted — 980 chars of source]
table[table omitted — 980 chars of source]
table[table omitted — 979 chars of source]

In accordance with the sampling thought experiment used for inference in Section (ref), in which I condition on group assignment $G_i$ and thus take it to be fixed, I create 500 aggregate data sets from each of the three true data-generating processes by, for each individual in group $g$, randomly drawing $X_{1i}, X_{2i},$ and $X_{3i}$ each according to the true marginal distributions for $X_{1i}|G_i=g, X_{2i}|G_i=g$, and $X_{3i}|G_i=g$ respectively. I then compute $Y_i$ according to the model above. Once $(X_i, Y_i)$ is drawn for all individuals in group $g$, I aggregate variables across individuals in each group to create a sample aggregate data set. Note that while I draw each sample $X_i$ as if $X_{1i}, X_{2i}$, and $X_{3i}$ are independent within each group, this is without loss of generality for the simulation because the covariance structure of $X_i$ within each group does not matter for the aggregate sample data set.

I consider two kinds of parameters across the simulation studies: conditional mean outcomes $\mathbb{E}[Y_i|X_i]$ and conditional outcome differences across $X_{1i}$, $\mathbb{E}[Y_i|X_{1i}=1,X_{2i},X_{3i}] - \mathbb{E}[Y_i|X_{1i}=0,X_{2i},X_{3i}]$. These are the same as the parameters that I consider in the empirical illustration. For each simulation study and parameter of interest I compute the population identified set using the population aggregate features without any additional restrictions and then with the monotonicity restrictions used later in the empirical illustration:

equation[equation omitted — 278 chars of source]

I then compute estimated bounds and 90% confidence intervals for each of the 500 sample data sets in each simulation study.

In Table (ref) I report results from the first simulation study. In particular, I report the population identified set, average coverage of the 90% confidence intervals, the ratio of the average width of the estimated bounds to the width of the identified set, and the ratio of the average width of the 90% confidence intervals to the width of the identified set, both without restrictions and under monotonicity. Bounds are not very informative across all parameters, although monotonicity does have some identifying power for certain parameters. As expected, coverage is conservative for all parameters. However, the average width of the confidence intervals is no more than 3% more than the width of the identified set. Similarly, the average width of the estimated bounds is essentially the same as the width of the identified set.

table[table omitted — 4,006 chars of source]

In Table (ref) I report results from the second simulation study analogous to those of Table (ref). Coverage is again conservative, although confidence intervals are not too much wider than the identified set on average, and estimated bounds are the same width as the identified set on average. Note from Table (ref) that the data is most representative of groups with high proportions of $X_{1i}$ and low proportions of $X_{2i}$ and $X_{3i}$. Bounds on the parameters corresponding to this subgroup---$\mathbb{E}[Y_i|X_{1i}=1,X_{2i}=0,X_{3i}=0]$ and $\mathbb{E}[Y_i|X_{1i}=1,X_{2i}=0,X_{3i}=0] - \mathbb{E}[Y_i|X_{1i}=0,X_{2i}=0,X_{3i}=0]$---are more informative in the second study compared to the first simulation study. Bounds on all other parameters are weakly less informative than the first simulation study.

This suggests that a feature of the aggregate data that shortens the width of bounds on conditional outcome parameters involving $X_i = x_k$ is how well-represented individuals with $X_{\ell i} = x_{k,\ell}$ for each covariate $\ell$ are in the data. For example, bounds on parameters involving $X_i = (1,0,0)$ are more informative in the second study than the first study because more individuals in the second study have each of $X_{1i} = 1, X_{2i} = 0,$ and $X_{3i} = 0$. Theoretically this feature helps because if marginal probabilities $\mathbb{P}[X_{\ell i} = x_{k,\ell}|G_i = g]$ are close to 1 across all $\ell$, there are fewer possible values the joint probability $\mathbb{P}[X_i=x_k|G_i=g]$ can take on. This reduces the width of the identified set for $\mathbb{P}[X_i=x_k|G_i=g]$ which then helps to narrow the bounds, as can be seen from the bilevel program representation of the identified set.

table[table omitted — 3,984 chars of source]

In Table (ref) I report results analogous to those above for the third simulation study. The third simulation study demonstrates the potential of my partial identification method to recover informative bounds. For example, bounds on the conditional outcome difference parameter $\mathbb{E}[Y_i|X_{1i}=1,X_{2i}=0,X_{3i}=0]-\mathbb{E}[Y_i|X_{1i}=0,X_{2i}=0,X_{3i}=0]$ do not contain zero. As before coverage is conservative although confidence intervals are not too much wider than the identified set on average, while estimated bounds are the same width as the identified set on average.

In this simulation study the data is representative of groups with both high and low proportions of $X_{1i}$, but only low proportions of $X_{2i}$ and $X_{3i}$, where high proportions of $X_{1i}$ occur with high expected outcome $Y_i$ and low proportions of $X_{1i}$ occur with low expected $Y_i$. As suggested by the results of the second simulation study, the parameters corresponding to these subgroups---$\mathbb{E}[Y_i|X_{1i}=0,X_{2i}=0,X_{3i}=0], \mathbb{E}[Y_i|X_{1i}=1,X_{2i}=0,X_{3i}=0]$, and $\mathbb{E}[Y_i|X_{1i}=1,X_{2i}=0,X_{3i}=0] - \mathbb{E}[Y_i|X_{1i}=0,X_{2i}=0,X_{3i}=0]$---have bounds that are much more informative relative to both the first and second simulation studies.

table[table omitted — 4,042 chars of source]

\paragraph{Comparison with cross2002regressions bounds} In Appendix (ref) I show how to obtain sharp bounds on $\mathbb{E}[Y_i|X_i=x_k,G_i=g]$ with aggregate data using the method proposed in cross2002regressions. I note in Remark (ref) that although one can use the cross2002regressions approach on $\mathbb{E}[Y_i|X_i=x_k,G_i=g]$ together with the identified set $P$ for the distribution of $X_i|G_i$ to obtain valid bounds on linear combinations of $\mathbb{E}[Y_i|X_i=x_k]$, these bounds are not sharp. Note that this alternative approach is also a nonconvex optimization problem for similar reasons that the approach I propose is nonconvex.

I demonstrate this numerically in Table (ref) for the three aggregate data sets used in the simulation studies. For all parameters considered in the simulation study, I display the sharp identified set obtained from my method (columns 1, 3, and 5) and the bounds obtained from following the approach described in Remark (ref) using cross2002regressions bounds. The sharp identified set is contained in the bound from the cross2002regressions approach for all parameters. In fact, for all but one of the parameters for which the identified set is strictly informative, the sharp identified set is a strict subset of the bound from the cross2002regressions approach.

table[table omitted — 3,393 chars of source]

Empirical Illustration

One setting in which publicly available data are in aggregate form is standardized exam data. In this section I apply the methodology developed in the previous sections to construct bounds on conditional exam pass rates and conditional white/non-white exam pass rate gaps. I also impose the three additional restrictions discussed in Remark (ref). I find, as suggested by cho2008cross, that without any additional assumptions, aggregate data does not have much identifying power for the individual-level parameter. Imposing monotonicity shape restrictions provides some identifying power, although bounds on pass rate gaps are still wide. Using additionally available pass rates by subgroup helps to narrow bounds more than imposing monotonicity shape restrictions. Adding monotonicity shape restrictions to the pass rates by subgroup further tightens bounds. A homogeneity restriction on conditional exam pass rates in different counties results in empty identified sets, suggesting that homogeneity of pass rates across counties is inconsistent with the observed aggregate data.

In this application I focus on exam pass rates for English and math Rhode Island Comprehensive Assessment System (RICAS) exams and student demographic information for the state of Rhode Island in the spring of 2019 over all students in grades 3-8. Data are obtained from the state of Rhode Island Department of Education's Public Assessment Data Portal, available at \url{https://www3.ride.ri.gov/ADP}. Although data is available at the school district level, for computational simplicity I aggregate the data up to the county level. There are 5 counties in Rhode Island, and correspondingly 5 groups in the data.

I estimate identified sets for pass rates and white/non-white pass rate gaps conditional on three covariates: race (indicator $white_i$ for being white, where non-white students are Hispanic, Black, Asian, Pacific Islander, Native American, or two or more races), economically disadvantaged status (indicator $econ_i$ for being economically disadvantaged, as defined by the Rhode Island Department of Education), and English-language learner (ELL) status (indicator $ELL_i$ for students who are currently English-language learners or were English-language learners in the past 3 years).

In addition to the aggregate data of Assumption (ref).iv, the Rhode Island Department of Education also makes available pass rates conditional on each of the three covariates for school districts with a sufficiently high number of students in each subgroup. I use the following three subgroup pass rates: $\mathbb{E}[pass_i|white_i=1,G_i=g], \mathbb{E}[pass_i|econ_i=0,G_i=g]$, and $\mathbb{E}[pass_i|ELL_i=0,G_i=g]$. I drop all districts that are missing any of these subgroup pass rates before aggregating to the county level. The subgroup pass rates imply additional restrictions as discussed in Remark (ref).

I first estimate bounds without any further assumptions, then impose additional monotonicity shape restrictions on the conditional pass rate. Motivated by test score gaps that have been documented between rich and poor students NYT2012, NYT2015, I consider the monotonicity restriction that for each value of $(white_i, ELL_i)$ and each county $g$ that the average pass rate is lower for economically disadvantaged students:

equation[equation omitted — 152 chars of source]

For English exams, I impose an additional monotonicity restriction that for each value of $(white_i, econ_i)$ and each county $g$ the average pass rate is lower for English-language learner students than for non-English-language learner students:

equation[equation omitted — 152 chars of source]

I then estimate bounds using the additional subgroup pass rates, without imposing monotonicity. Finally, I estimate bounds using the additional subgroup pass rates under monotonicity.

I first present estimated bounds on math exam white/non-white pass rate gaps in Table (ref) and on English exam white/non-white pass rate gaps in Table (ref). 90% confidence intervals constructed as in Section (ref) are displayed below the bound estimates in curly brackets. For both types of exam, bounds on the white/non-white pass rate gaps reported in Column 1 are wide and either uninformative or close to uninformative. Imposing the additional monotonicity restrictions helps narrow the bounds for some parameters, reported in Column 2, but bounds are still wide. Using additional subgroup pass rates without monotonicity in Column 3 further narrows the bounds for those same parameters. Adding monotonicity to additional subgroup pass rates in Column 4 produces the tightest bounds, especially for white/non-white pass rate gaps among students who are economically disadvantaged but not English-language learners. For example, among economically disadvantaged students who are not English-language learners, the white/non-white English exam pass rate gap is estimated to be no lower than -30% and no higher than 64%. Bounds on pass rate gaps among English-language learners are never informative, as expected due to the small percentage of English-language learners in Rhode Island.

table[table omitted — 2,647 chars of source]
table[table omitted — 2,680 chars of source]

In Table (ref) I present estimated bounds on math exam conditional pass rates and in Table (ref) I present estimated bounds on English exam conditional pass rates. 90% confidence intervals constructed as in Section (ref) are displayed below the bound estimates in curly brackets. Just as for the white/non-white pass rate gaps, the additional restrictions help to narrow the bounds on pass rates for students who are not English-language learners. For example, the math exam pass rate among non-white students who are economically disadvantaged and not English-language learners is estimated to be bounded above by 29%, while the math exam pass rate among white students who are not economically disadvantaged and not English-language learners is estimated to be between 37% and 72%. Similarly, the English exam pass rate among non-white students who are economically disadvantaged and not English-language learners is estimated to be bounded above by 36%, while the math exam pass rate among white students who are not economically disadvantaged and not English-language learners is estimated to be between 50% and 80%.

table[table omitted — 3,133 chars of source]
table[table omitted — 3,171 chars of source]

While each of these bounds are sharp for each of the corresponding parameters, as proved in Proposition (ref), note that the bounds are not jointly sharp across all parameters in the tables. In fact, the Minkowski difference of bounds on conditional pass rates in Tables (ref) and (ref) are sometimes wider than the bounds on the corresponding pass rate gaps in Tables (ref) and (ref), especially with additional subgroup pass rates under monotonicity. If the bounds were jointly sharp across all parameters in Tables (ref) and (ref) then the Minkowski difference of the bounds on conditional pass rates should be equal to the bounds on the corresponding pass rate gaps. That the Minkowski difference is wider highlights an advantage of my method, which is that my method directly produces sharp bounds for the linear combination of conditional pass rates.

Another restriction that I could impose, as discussed in Remark (ref), is a homogeneity restriction on the conditional pass rate across counties. This restriction imposes that for every value of $(white_i, econ_i, ELL_i)$, $\mathbb{E}[pass_i|white_i, econ_i, ELL_i, G_i = g] = \mathbb{E}[pass_i|white_i, econ_i, ELL_i, G_i = g']$ for all pairs of counties $g$ and $g'$. This assumption, which is very strong, can be interpreted as imposing that county exam pass rates differ only because of differences in student characteristics and not because of, for example, differences in county value-added to student education. When I estimate bounds for both math and English exam white/non-white pass rate gaps, I obtain empty sets. This means that homogeneity is inconsistent with the observed data, as discussed in Remark (ref).

Conclusion

In this paper I consider the problem of identifying linear combinations of conditional mean outcomes when the data that is observed is marginal information on each individual-level variable over many groups, which I call aggregate data. This model of aggregate data necessitates a nontrivial extension of existing methods in the ecological inference literature. I develop a partial identification methodology to construct bounds by solving an optimization problem that considers all joint distributions of individual-level variables that are consistent with the observed marginal information. This approach easily accommodates additional model restrictions like polyhedral shape restrictions. I propose a plug-in estimator of the identified set that is consistent and a valid method to construct confidence intervals that uses aggregate data only. I demonstrate the performance of the method in simulation studies and apply the method to Rhode Island standardized exam data.

Acknowledgements

I am especially grateful to Edward Vytlacil, Isaiah Andrews, and Anna Mikusheva for their guidance and advice. I thank Editor Yuya Sasaki, the Associate Editor, two anonymous referees, Alberto Abadie, Hyungsik Roger Moon, Whitney Newey, Jesse Shapiro, Elie Tamer, Vod Vilfort, and participants in the MIT econometrics lunch seminar for their helpful comments and suggestions. All errors are my own.

Funding

This paper is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. 2141064.

Disclosure Statement

There are no competing interests to declare.

\singlespacing