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.
108,306 characters · 21 sections · 42 citation commands
Partial Identification under Missing Data Using Weak Shadow Variables from Pretrained Models
\RUNAUTHOR{Chen, Simchi-Levi, and Xiong} \RUNTITLE{Partial Identification under MNAR Data}
\TITLE{Partial Identification under Missing Data Using Weak Shadow Variables from Pretrained Models}
\ARTICLEAUTHORS{ \AUTHOR{Hongyu Chen\textsuperscript{1} \quad\quad David Simchi-Levi\textsuperscript{1} \quad\quad Ruoxuan Xiong\textsuperscript{2}} \AFF{\textsuperscript{1}Massachusetts Institute of Technology, Cambridge, MA 02139\texorpdfstring{\\} \textsuperscript{2}Emory University, Atlanta, GA 30322\texorpdfstring{\\} \EMAIL{[email removed], [email removed], [email removed]}} }
\ABSTRACT{ Estimating population quantities such as mean outcomes from user feedback is fundamental to platform evaluation and social science, yet feedback is often missing not at random (MNAR): users with stronger opinions are more likely to respond, so standard estimators are biased and the estimand is not identified without additional assumptions. Existing approaches typically rely on strong parametric assumptions or bespoke auxiliary variables that may be unavailable in practice. In this paper, we develop a partial identification framework in which sharp bounds on the estimand are obtained by solving a pair of linear programs whose constraints encode the observed data structure. This formulation naturally incorporates outcome predictions from pretrained models, including large language models (LLMs), as additional linear constraints that tighten the feasible set. We call these predictions weak shadow variables: they satisfy a conditional independence assumption with respect to missingness but need not meet the completeness conditions required by classical shadow-variable methods. When predictions are sufficiently informative, the bounds collapse to a point, recovering standard identification as a special case. In finite samples, to provide valid coverage of the identified set, we propose a set-expansion estimator that achieves slower-than-$\sqrt{n}$ convergence rate in the set-identified regime and the standard $\sqrt{n}$ rate under point identification. In simulations and semi-synthetic experiments on customer-service dialogues, we find that LLM predictions are often ill-conditioned for classical shadow-variable methods yet remain highly effective in our framework. They shrink identification intervals by 75--83% while maintaining valid coverage under realistic MNAR mechanisms. }
\KEYWORDS{partial identification; missing not at random; shadow variables; large language models; linear programming}
Missing data is pervasive in economic and social research as well as on digital platforms. In household and health surveys, respondents often skip questions perceived as sensitive or irrelevant. On digital platforms, users often choose whether or not to leave feedback based on their experiences. As noted by abrevaya2017gmm, nearly 40% of top economics papers report data missingness, with about 70% dropping observations as a result.
In many of these settings, data is missing not at random (MNAR): the probability of observing an outcome depends on its possibly unobserved value. For example, bollinger2019trouble shows that nonresponse across the earnings distribution is U-shaped, where left-tail “strugglers” and right-tail “stars” are least likely to report earnings. Conversely, an inverse U-shaped missingness pattern can be found in online product reviews, where users with more extreme opinions are more likely to leave reviews hu2017self. This dependence between missingness and actual outcome creates a fundamental challenge for accurate estimation and decision-making. Estimators based solely on observed data without modeling the missing mechanism can be severely biased, motivating approaches that explicitly account for the missingness mechanism.
In this paper, we study the problem of identifying population quantities, such as the mean outcomes, when data are MNAR. Such questions are prevalent in service platforms or social surveys, e.g., when a platform seeks to evaluate average customer satisfaction, or when a researcher aims to estimate average income in a particular region. For those questions, one class of classical methods addresses this MNAR problem by imposing strong parametric structural assumptions, such as those in the Heckman selection model heckman1979sample or the Pattern-Mixture model rubin1987calculation, little1994class. Another common approach introduces auxiliary variables, including instrumental variables d2010new or shadow variables miao2016varieties, which need to satisfy restrictive independence or completeness conditions for identification. Both strategies can face practical limitations: structural parametric models may be misspecified, and identifying valid auxiliary variables may require substantial domain expertise or serendipity.
We therefore take a different approach and study: under realistically minimal assumptions, what can we still learn about population quantities like the mean outcomes? Instead of seeking point identification, we adopt a partial identification perspective manski2003partial and aim to characterize sharp upper and lower bounds on the estimand (e.g., mean outcome). Our key insight is that this problem can be reformulated as a pair of linear programs (LPs). In this formulation, the objective corresponds to the estimand, while constraints encode the probabilistic structure implied by the observed data. This yields a transparent and tractable framework for estimation under MNAR.
While bounds obtained under minimal assumptions are valid, they can be wide, especially when a large portion of the data is missing. To tighten the bounds, we propose incorporating auxiliary predictions from modern machine learning systems, such as large language models (LLMs). Recent work suggests that LLMs exhibit human-like reasoning abilities and can approximate decision-making in complex settings horton2023large,goli2024frontiers,brand2024using, making them promising candidates for predicting unobserved outcomes, e.g. predicting user satisfaction from chat transcripts. Importantly, because such predictions are generated by an external model rather than by the individuals themselves, they do not directly influence whether an outcome is observed; that is, they naturally satisfy an exclusion-type condition with respect to the missingness mechanism. At the same time, discrepancies between LLM outputs and actual human behavior have been documented gui2023challenge,li2024frontiers,gao2025take, and researchers have cautioned against assuming that model predictions can perfectly substitute for human judgments. These observations suggest that LLM predictions can serve as useful auxiliary signals for tightening identification bounds, but that an approach robust to prediction imperfections is needed.
In light of both the promise and limitations of these predictions, we treat LLM-generated outputs as weak shadow variables. Specifically, we assume an exclusion-type condition: conditional on the true outcome and observed covariates, the prediction is independent of the missingness indicator. However, we do not require strong relevance or completeness conditions as in classical shadow variables miao2016varieties; the predictions may only weakly correlate with the outcome. Even so, incorporating them introduces additional linear constraints into our identification framework, tightening the feasible region. When the predictions are sufficiently informative, the bounds may collapse to a single point, yielding point identification as a special case.
We summarize our three main contributions below.
First, we propose a novel linear programming framework for partial identification under MNAR that applies both with and without auxiliary predictions. In the baseline setting without auxiliary inputs, the formulation yields closed-form solutions for the identification region of the mean outcomes. When incorporating auxiliary predictions from LLMs, we derive analytical results that quantify how these predictions tighten the feasible set and shrink the identification region. This formulation offers a unified and tractable approach to understanding how predictive signals impact identification under minimal assumptions.
Second, we develop a set-expansion estimator that asymptotically converges to the identification region while accounting for estimation error in the probability constraints that define the bounds. We establish convergence rates for this estimator. In the partially identified setting, the convergence rate is slower than the usual $\sqrt{n}$ rate (for example, on the order of $\sqrt{n}/\log n$), as a result of the additional uncertainty inherent in set identification. When the auxiliary information is sufficiently informative and point identification is achieved, the estimator recovers the standard $\sqrt{n}$ rate, matching classical results in the shadow variable literature miao2016varieties. We further extend our framework to randomized experiments and derive bounds on treatment effects, and provide sufficient conditions under which reliable treatment decisions can still be made despite missing outcomes.
Third, we evaluate the proposed methods through simulation studies and semi-synthetic experiments based on real customer-service dialogue data. To construct auxiliary signals, we generate outcome predictions using LLMs under several prompting and training regimes, including zero-shot, few-shot, chain-of-thought prompting, and fine-tuning. Our results reveal two key insights. First, we show that LLM predictions can fail to meet the strong completeness conditions required for classical shadow variable methods, rendering point identification unstable or infeasible. This underscores the need for our partial identification framework. Second, despite these limitations, the predictions remain informative: incorporating LLM-based weak shadow variables reduces the width of identification intervals by 75--83% across prompting strategies, while preserving valid coverage under realistic MNAR mechanisms.
Our work contributes to the rich literature on identification and estimation under MNAR mechanisms. Classical approaches include parametric selection models, such as the Heckman correction, which jointly models the outcome and missingness process heckman1979sample, and Pattern-Mixture models that parameterize outcome distributions within each missingness stratum little1994class,rubin1987calculation. Other strands of work leverage graphical models to represent missing data processes fay1986causal, or use auxiliary variables such as instrumental variables that affect missingness but not outcomes das2003nonparametric,tchetgen2017general,sun2018semiparametric. Our approach is most closely aligned with recent developments in the shadow variable literature d2010new,miao2016varieties,miao2024identification, which typically uses the odds ratio for the identification of the distribution of missing outcomes. We contribute to this line of research in three key ways. First, we introduce a novel linear programming framework that characterizes the identification region for mean outcomes under MNAR. Second, we generalize the shadow variable approach by allowing weak shadow variables; this enables the use of auxiliary signals, e.g., from LLMs, that may violate classical completeness assumptions and thus do not yield point identification but can still significantly tighten bounds. Third, we establish convergence rates for our estimated identification region under both partial and point identification regimes.
Our linear programming formulation connects to the broader literature on inference for partially identified models manski2003partial, imbens2004confidence. ChernozhukovHongTamer2007 proposed a criterion-function approach with set expansion to construct confidence regions for identified sets, which directly inspires our estimator. beresteanu2008asymptotic connect identified sets to LP optimal values through a support function characterization, and mogstad2018using and kaido2019confidence develop LP-based inference for treatment effect bounds and subvector projections, respectively. Our work derives a specific LP structure from the shadow variable assumption under MNAR and shows that auxiliary predictions generate additional constraints that tighten the identified set, with convergence rates that adapt to whether the shadow variable yields partial or point identification.
Our work also relates to the growing literature on leveraging pretrained models as auxiliary signals to improve identification or statistical efficiency. Prediction-powered inference (PPI) methods angelopoulos2023prediction,angelopoulos2023ppi++ assume true labels are observed for only a random subset of the data, when predictions from an external model are available for the remainder, and aim to combine the two sources to enable valid inference. ji2025predictions propose PPI with “recalibrated” prediction, learning a map from the model prediction and covariates to the true outcome to correct bias. wang2025efficient further propose optimal sample allocation strategies that first fine-tune LLMs and then apply PPI to correct for prediction bias. From a different perspective, wang2024large explore how LLM-generated simulations, when grounded in real data, can support accurate conjoint analysis. chen2025utilizing further examine how to design data collection and efficient inference strategies in the presence of such LLM-based predictors. Our work differs in two key ways. First, we explicitly account for MNAR missingness. Second, we interpret auxiliary predictions as weak shadow variables, leading to a framework that provides valid bounds on population quantities, rather than relying on point estimates that require stronger missingness assumptions.
The remainder of the paper is organized as follows. Section (ref) introduces the problem setup. Section (ref) develops the linear programming framework for partial identification, both with and without weak shadow variables. Section (ref) presents the set-expansion estimator and its convergence properties. Section (ref) extends the framework to randomized experiments. Section (ref) reports simulation and semi-synthetic experiments, and Section (ref) concludes.
Suppose we are evaluating an economic, social, or digital system that solicits discrete feedback, such as survey responses, program satisfaction ratings, or customer-service reviews on online platforms. Outcomes are observed only when individuals choose to respond. Because not all individuals provide feedback, outcomes are partially observed. Let $R \in \{0,1\}$ indicate whether a user's rating is observed ($R = 1$) or missing ($R = 0$). The rating is denoted by $Y \in [M]$, where $[M] = \{1, \ldots, M\}$ represents a discrete set of possible scores.
We focus on the setting where the missingness is not at random, meaning that the probability of observing a rating may depend on its value, that is, $R \not\indep Y$. This is motivated by empirical evidence that users with extreme, either positive or negative, ratings are more likely to write reviews than users with moderate product ratings hu2017self. Thus, our observation is a set of i.i.d. data $\{R_i, R_iY_i\}_{i=1}^n$ with $Y_i$ observed if and only if $R_i=1$, and $n$ is the number of observations.
Our primary objective is to estimate the mean outcome \[\theta = \E[Y] \] such as the average rating across all customers. While we focus on the population mean for concreteness, our framework extends directly to other population quantities, including functionals of the form $\E[g(Y)]$ and other distributional summaries, as discussed later.
However, when outcomes are missing not at random, the observed ratings may not represent the full population, and a naive empirical average over observed data is typically biased. As a result, point identification of $\theta$ is generally infeasible without additional assumptions on the missingness mechanism, e.g., via introducing additional variables or a parametric model on missing pattern heckman1979sample, little1994class. These underlying structural assumptions are generally untestable from observed data and may also lead to biased estimates. In this paper, we study this problem in a fully non-parametric way with minimal structural assumptions. In particular, we adopt a partial identification perspective: rather than aiming for a single estimate under tenuous assumptions, we characterize the range of values $\theta$ can plausibly take given the observed data.
In this section, we develop a linear-programming characterization of partial identification under MNAR. We first show that the sharp identification set for the estimand $\theta$ reduces to an explicit interval whose endpoints solve a simple pair of linear programs. We then show how predictions from a pretrained model can be incorporated as weak shadow variables and yield additional linear constraints to tighten the identification set.
Our partial identification strategy is based on the following decomposition of the population mean: \[\theta = \sum_{y=1}^M y \cdot \P(Y = y) = \sum_{y=1}^M y \cdot \underbrace{\P(Y=y, R=1)}_{:=\alpha(y)} \big/ \underbrace{\P(R=1\mid Y=y)}_{:=\pi(y)} \,. \] If we are interested in other population quantities (e.g., $\E[g(Y)]$), then we replace $y \cdot \P(Y = y)$ by $g(y) \cdot \P(Y = y)$ in the decomposition. Here, the joint probability $\alpha(y)=\P(Y=y, R=1)$ is identifiable from observed data. However, the conditional response probability $\pi(y) = \P(R=1\mid Y=y)$ is generally unidentifiable when missingness depends on the outcome itself. As a result, the mean $\theta$ cannot be point-identified without further assumptions.
We characterize the sharp identification region for $\theta$ by considering all possible values of $\pi(y) \in (0,1]$. Note that the only constraint from observational data is that the probabilities $\P(Y=y)$ must sum to one. Thus, we can define the feasible set for $\pi(y)$ as: \[ \Pi= \left\{ (\pi(1), \dots, \pi(M)): \sum_{y=1}^M \alpha(y) / \pi(y) = 1, \pi(y)\in(0,1]\right\} \,. \]
This induces the identification set for the mean outcome: \[\Theta =\left\{ \sum_{y=1}^M y \cdot \alpha(y)/\pi(y): (\pi(y))_{y\in[M]}\in\Pi\right \}.\]
To simplify notations, let $w(y)=1/\pi(y)-1$. Under this change of variables, the feasible region becomes a polyhedron in $w(y)$, and the mapping from $w(y)$ to $\theta$ is linear. Thus, the identification region $\Theta$ is a closed interval, and its endpoints can be computed by solving the following pair of linear programs:
Therefore, the identified region for the mean outcome is given by $\Theta = [\theta_{\min}, \theta_{\max}]$. In the proposition below, we show that both $\theta_{\min}$ and $\theta_{\max}$ can be solved analytically.
Here we implicitly let $\P(R=1)\E[Y\mid R=1]$ be zero if\; $\P(R=1)$ is zero and hence $\E[Y\mid R=1]$ is undefined. These expressions are attained by setting the weights $w(y)$ to their minimum allowable value $w(y) = 0$ for all but one outcome level. To achieve $\theta_{\min}$, we set $w(y) = 0$ for $y = 2, \ldots, M$ and assign the remaining mass to $y = 1$, the smallest outcome. Conversely, to achieve $\theta_{\max}$, we set $w(y) = 0$ for $y = 1, \ldots, M-1$ and concentrate the remaining weight on $y = M$, the largest outcome. This corresponds to placing as much probability mass as possible on the lowest or highest feasible rating levels, subject to the constraint induced by the observed joint distribution $\alpha(y) = \P(Y = y, R = 1)$.
Notably, the width of the identification region is $\theta_{\max} - \theta_{\min} = (M-1) \cdot \P(R = 0)$, which scales linearly with the probability of missingness. When $\P(R = 0) = 0$, i.e., outcomes are fully observed, the bounds collapse to a point and $\theta$ is point-identified. In contrast, when $\P(R = 0) = 1$, the bounds are equal to the full support range, $[1, M]$, which is uninformative. Without additional information and structural assumptions, the identification region in Proposition (ref) is the best one can hope for. They are sharp bounds for identification in the sense that any other valid identification region from observed data will contain $\Theta$ as a subset.
In many applications, additional covariates are available. For example, online platforms and social science studies often collect structured demographic information such as age, race, or location, and a common practice is to stratify analyses across such covariate groups to reduce heterogeneity and improve precision. Motivated by this, we study how the inclusion of covariates affects the identification region $\Theta$. To build intuition, we focus on the setting where covariates take on finitely many values and can be stratified into discrete groups, although the same ideas can extend to continuous covariates. Let $X \in \mathcal{X} = \{x_1, \ldots, x_K\}$ and consider stratifying the analysis by each covariate level. The mean outcome can be decomposed as \[\theta = \sum_{x \in \mathcal{X}} \sum_{y=1}^M y \cdot \frac{\P(Y=y, R=1\mid X = x)}{\P(R=1\mid Y=y, X = x)} \P(X = x) \,. \] Under the same procedure, we can derive a sharp identification region for the target variable $\theta$.
Interestingly, these bounds coincide exactly with those obtained from the unstratified formulation. Hence, stratifying on covariates does not tighten the identification region. Intuitively, this is because the extreme mass allocations (to outcome $y = 1$ or $y = M$) are feasible within each stratum independently and carry through under marginalization. This motivates the necessity of extra variables with richer structural properties for sharper identification. In the next section, we will consider additional predictions as a covariate that satisfy a conditional independence structure. Under that scenario, we will observe a sharpened identification region.
In many modern settings, we have access to rich unstructured data, such as chat transcripts, reviews, or other interaction logs, that plausibly encode information about the missing outcomes. Directly incorporating such high-dimensional inputs into stratified analyses is typically infeasible. Instead, we propose to leverage pretrained foundation models to extract low-dimensional predictive signals from these data. In particular, we feed them into a pretrained model (e.g., an LLM) to produce an output $F \in \mathcal{F}$ that serves as a proxy for how a human would rate the interaction, thereby injecting additional information about the latent outcome distribution. In practice, we can prompt the model to return a discrete rating (e.g., “choose one of $\{1, \ldots, M\}$”), so $\mathcal{F}$ is finite and, without loss of generality, we can take $\mathcal{F} = [M]$ to match the support of $Y$. We assume $F$ is always observed along with additional covariates $X$. Thus, the final dataset takes the form $\{F_i, X_i, R_i, R_iY_i\}_{i=1}^n$.
We impose the following assumption on the relationship between the model output and the missingness mechanism:
We call $F$ a {Weak Shadow Variable} if it satisfies Assumption (ref) and its connection to classic shadow variable will be clear in the next section. This is a relatively weak assumption that does not require the predictive model to be accurate. It states that, once we condition on the actual outcome and covariates, the model output provides no additional information about whether the outcome is observed. This condition holds in most applications, as users' decisions to respond typically do not depend on an external predictor that they never observe. As a sanity check, if the prediction is perfect, i.e., $F = Y$, then Assumption (ref) is trivially satisfied. A causal diagram for the relationship can be found in Figure (ref).
Assumption (ref) is closely related to the shadow variable/auxiliary variable framework, which has been proposed as an alternative to instrumental variable approaches in the literature on nonrandom missing data miao2016varieties, miao2024identification. In this section, we clarify our connections and how our definition for weak shadow variable generalizes the traditional approach.
The shadow variable framework provides a pathway to handle MNAR data and has been widely adopted in empirical studies ibrahim2001using, kott2018calibration. However, the existence of a shadow variable alone is not sufficient for point identification of the estimand $\theta$. In addition to the two conditions in Definition (ref), prior work requires completeness of the conditional distribution $\P(Y \mid R=1, X, F)$ as another constraint to ensure point identification miao2016varieties, miao2024identification.
The completeness condition is a strengthened version of the dependence constraint in Definition (ref)---not only must the shadow variable $F$ and the outcome $Y$ be dependent, but the variation in $F$ must be sufficient to explain the variation in $Y$. Under our discrete outcome setup, the completeness condition admits a cleaner interpretation in terms of matrix algebra. Define the joint distribution matrix $H = [P(F=f, Y=y)]_{f \in \mathcal{F}, y \in [M]} \in \mathbb{R}^{|\mathcal{F}| \times M}$ and the conditional distribution matrix $B = [P(Y=y \mid F=f, R=1)]_{f \in \mathcal{F}, y \in [M]} \in \mathbb{R}^{|\mathcal{F}| \times M}$. Note that $H$ is independent of the missingness scheme but $B$ is related to the distribution of $R$. Completeness then corresponds directly to a full-rank condition on the matrix $B$, which is further equivalent to a rank condition on $H$ under mild assumptions.
Proposition (ref) translates the abstract completeness condition into a concrete matrix rank condition, which is much simpler to interpret and verify in practice. The condition number $\kappa(H)$ quantifies the degree of invertibility of the joint distribution matrix: when $\kappa(H)$ is close to one, the matrix is well-conditioned and completeness holds strongly; when $\kappa(H)$ is large, the matrix is nearly rank-deficient and completeness barely holds. We will use $H$ rather than $B$ to evaluate the prediction in later experiment section as it does not depend on the missingness scheme.
Point identification using shadow variables relies on solving a Fredholm integral equation, whose numerical stability depends critically on the condition number of matrix $B$. In modern applications where predictive models are generic rather than problem-specific---as is the case with large language models---the predictions may not fully capture variation in the outcome, resulting in ill-conditioned matrices $H$ and, by Proposition (ref), ill-conditioned matrices $B$. As we demonstrate in Section (ref), LLM-generated predictions for real customer service data yield condition numbers on the order of $10^3$ and minimum singular values around $10^{-4}$, rendering point identification numerically unstable.
Our definition of a weak shadow variable departs from the classical shadow variable literature in an important way: we fully relax both the dependence and completeness conditions. We show that neither condition is necessary to identify bounds on $\theta$. Moreover, the proposed identification bounds adapt naturally to the strength of correlation between $Y$ and $F$ conditional on $X$ and $R=1$. Stronger correlation leads to tighter identification regions, and completeness recovers point identification as a limiting case. Another advantage of the partial identification approach is computational: rather than solving potentially ill-conditioned integral equations, we solve a linear program with relaxed constraints, which is numerically more stable. This perspective allows practitioners to leverage modern predictive models as auxiliary tools for tightening identification regions without imposing stringent requirements on predictive accuracy.
In this section, we proceed to develop the sharp identification region for $\theta$ with weak shadow variables. Here our partial identification strategy is based on stratification $\theta = \sum_{x\in \mathcal{X}} \theta_x \cdot \P(X = x)$, where $\theta_x = \mathbb{E}[Y \mid X = x]$ is the conditional mean within stratum $x$. We will first provide an identification region for every $x \in \mathcal{X}$ and then aggregate them together to obtain an identification region for $\theta$. For each $x \in \mathcal{X}$, we have the following decomposition
The second equality follows from the chain rule of probability. Under Assumption (ref), i.e., $F\indep R\mid Y$, we have $ \P(R = 1 \mid Y = y, X = x) = \P(R = 1 \mid F = f, Y = y, X = x)$ for all $f \in \mathcal{F}$, so the denominator does not depend on $f$. We therefore write it as $\pi_x(y)$ for notation simplicity.
The quantity $\alpha_x(f, y)$ is identifiable from observed data, but $\pi_x(y)$ remains unidentifiable. We therefore propose to identify a set of feasible values for $\pi_x(y)$. Here, we leverage the following identity:
where the second equality uses Assumption (ref). We let $\beta_x(f)=\P(R=0, F=f\mid X=x)$ for notation simplicity. Note that $\beta_x(f)$ is identifiable from observed data. Thus, for each $x \in \mathcal{X}$, we can similarly write the feasible region for $\pi_x(y)$ as \[ \Pi_x= \left\{ (\pi_x(1), \dots, \pi_x(M)): \sum_{y=1}^M\alpha_x(f,y)\left(\frac{1}{\pi_x(y)}-1\right)=\beta_x(f)\,,~~~\;\forall f \in \mathcal{F}, \pi_x(y)\in(0,1] \right\} \] and the identification region for $\theta_x$ becomes \[ \Theta_x=\left\{ \sum_{f \in \mathcal{F}} \sum_{y=1}^M y \frac{\alpha_x(f,y)}{\pi_x(y)}: (\pi_x(y))_{y\in[M]}\in\Pi_x \right\} \,. \]
Letting $w_x(y)=1/\pi_x(y)-1$, we obtain a linear representation of the objective and constraints in terms of $w_x(y)$. The feasible set $\Pi_x$ remains convex, and so the identification region $\Theta_x$ is a closed interval. The endpoints are given by the solution to the following pair of linear programs:
where $A_x=[\alpha_x(f,y)]_{f,y}\in [0,1]^{|\mathcal{F}|\times M}$, $\mathbf{w}_x=(w_x(1), \dots, w_x(M))^\top$, and $\bm{\beta}_x=(\beta_x(1), \dots,\beta_x(|\mathcal{F}|))^\top$, $D=\mbox{diag}\{1,2,\ldots,M\}$, and $\bm{1}$ and $\bm{0}$ are vectors of all ones and zeros, respectively. The constraints here restrict $\pi_x(y)$ to lie in $\Pi_x$ where the objective function is the definition of $\theta\in\Theta_x$. Aggregating over covariate strata, we obtain the sharp identification region for the estimand $\theta$.
The identification of $\theta$ depends on the identification of each stratum mean $\theta_x$, which in turn is governed by the linear system $A_x$. If $A_x$ has full column rank, i.e., $\operatorname{rank}(A_x)=M$, then it corresponds to the completeness condition in Definition (ref). In this scenario, the linear constraint system has a unique solution, and we achieve point identification: $\theta_{x,\min} = \theta_{x,\max}$. More generally, if some rows of $A_x$ are linearly dependent (e.g., $F \indep Y \mid X, R = 1$ for every $F \in \mathcal{F}'$ for some subset $\mathcal{F}'\subset\mathcal{F}$) or if the columns are dependent (e.g., $F \indep Y \mid X, R = 1$ for a subset of outcomes $Y \in \mathcal{Y}' \subset \{1, \ldots, M\}$), then the feasible region contains multiple solutions and $\mathbf{w}_x$ is only partially identified. In this sense, our formulation generalizes the classical shadow variable approach: it allows violations of the completeness condition of the shadow variable definition and quantitatively characterizes how the strength of association between $F$ and $Y$ impacts the width of the identification region $\Theta_x$. Again, the identification in Proposition (ref) is sharp in the sense that any other valid identification region based on the observed data will contain $\left[\theta_{\min}^{\mathrm{shad}}, \theta_{\max}^{\mathrm{shad}}\right]$ as a subset.
Lastly, we compare the above identification region with the one defined in Equation (ref) where the shadow variable is not available to understand the effect of the additional prediction $F$. Note that the formulation in linear program (ref) is closely related to linear program (ref), where the constraint for the lower bound in linear program (ref) can be written as a single aggregated constraint $\bm{1}^\top A_x \mathbf{w}_x = \bm{1}^\top \bm{\beta}_x $. Thus, we can use the techniques in aggregation bounds zipkin1980bounds, litvinchev2013aggregation to analyze their differences.
Theorem (ref) gives a characterization of the difference between the identification bound with and without a shadow variable. The proof is provided in Appendix (ref). As a special case, we have the identification with shadow variable is never worse than that without it, i.e., $\theta_{\min}\leq\theta_{\min}^{\mathrm{shad}}\leq\theta_{\max}^{\mathrm{shad}}\leq\theta_{\max} $. Moreover, the amount of improvement depends on the missingness ratio (represented be $\bm{1}^\top\bm{\beta}_x$) and the worst misalignment between missingness data distribution and observed data distribution (characterized by $\left\|\frac{\bm{\beta}_x}{\bm{1}^\top\bm{\beta}_x}-\frac{\bm{a}_{x,M}}{\bm{1}^\top\bm{a}_{x,M}}\right\|_1$). Thus, the shadow variable is especially useful when missingness level is high and when the missingness scheme does not align well with observed scheme, in which case the data is far from MAR.
In practice, when solving the linear programs in Equation (ref), we do not have access to the true probability quantities $A_x = [\alpha_x(f,y)]_{f,y}$ and probability vector $\bm{\beta}_x = [\beta_x(f)]_f$. Instead, we estimate them from observed data and then solve the optimization problem using the estimates to obtain $\hat \theta_{x,\min}$ and $\hat \theta_{x,\max}$. This may cause problems as the system of linear equations may be over-identified, resulting in unstable finite sample performance. In this section, we propose a set expansion estimator that ensures the solution from the linear program in finite sample is a valid estimation for the estimand.
Let $\mathcal{D} = \{(X_i, R_i, R_i Y_i, F_i)\}_{i=1}^n$ denote a random sample of size $n$, and let $n_x$ denote the number of units in stratum $x$, i.e., $n_x = \sum_{i=1}^n \bm{1}(X_i = x)$. We assume that the estimators for $\alpha_x(f,y)$ and $\beta_x(f)$ converge at the $\sqrt{n_x}$ rate, as formalized in the assumption below.
This assumption is mild and holds under standard empirical or maximum likelihood estimators, given the finite support of both $f$ and $y$.
Given the estimators $\hat A_x$ and $ \hat{\bm\beta}_x$, a natural approach is to plug them into the LP formulation and solve the empirical analog of (ref). However, this direct plug-in method can lead to infeasibility or unstable solutions. This is because the matrix $A_x$ may have linearly dependent rows (depending on the informativeness of $F$), so even small estimation errors may perturb the feasible region in a way that makes the empirical LP infeasible, even though the population LP has a valid solution. To address this issue and ensure well-posed optimization, we introduce a bounded constraint on the solution space:
This boundedness assumption equivalently implies a lower bound on the conditional observation probability: $\pi_x(y) = \P(R = 1 \mid Y = y, X = x) \geq 1/(C+1)$. Such positivity conditions are standard in the missing data literature. Without such a lower bound, estimation of the outcome distribution for rarely observed strata becomes ill-posed, as there would be insufficient information in the observed data to learn about the missing values.
Next, we propose to use a set expansion approach, inspired by ChernozhukovHongTamer2007, to ensure that the empirical analog of (ref) remains feasible even when using the estimators $\hat A_x$ and $ \hat{\bm\beta}_x$. Specifically, we quantify the violation of estimated linear constraints within the bounded solution space by computing the following diagnostic: \[\hat m_x=\min_{0\leq\mathbf{w}_x\leq C}\|\hat{A}_x \mathbf{w}_x-\hat{\bm{\beta}}_x\|_{\infty}\,.\]
This quantity captures the minimax deviation from the constraint $\hat{A}_x \mathbf{w}_x=\hat{\bm{\beta}}_x$ over all candidate solutions within the bounded region $0\leq\mathbf{w}_x\leq C$. Intuitively, $\hat m_x$ measures the degree of infeasibility introduced by the estimation error. If the empirical linear program is exactly feasible, then $\hat m_x =0$; otherwise, $\hat m_x > 0$ reflects the minimal slack needed to make the LP feasible again.
Using $\hat m_x$ and a user-specified tolerance sequence $\kappa_{n_x}>0$, we define the set expansion estimator of $\theta_{x,\min}$ as the solution to the following relaxed linear program:
Here, adding $\hat{m}_x$ ensures feasibility, and the set expansion margin ${\kappa_{n_x}}/{\sqrt{n_x}}$ accounts for sampling variability, leading to correct coverage of the target value. The estimator $\hat \theta_{x,\max}$ is defined analogously by replacing the minimization with a maximization in the linear program (ref).
We show that the set expansion estimator is consistent for the true bounds $\theta_{x,\min}$ and $\theta_{x,\max}$.
The convergence rate depends on the set expansion margin $\kappa_{n_x}/\sqrt{n_x}$. Following ChernozhukovHongTamer2007, a typical choice is $\kappa_{n_x} = \log n_x$ or $\log\log n_x$. We recommend choosing a slower choice for $\kappa_{n_x}$ as we have added $\hat{m}_x$ in the constraint for feasibility, which is different from ChernozhukovHongTamer2007. However, if stronger structural conditions hold, such as full column rank of $A_x$, then a smaller slack suffices. This leads to the improved convergence rates, as formalized below.
In this setting, the variable $F$ is sufficiently informative to point identify $\theta_x$. Our estimator achieves the same convergence rate as those derived under shadow variable conditions for point identification of $\theta$ miao2024identification. Hence, our framework generalizes existing results by accommodating both point-identified and partially identified regimes and showing the convergence rates in each regime. The proof for both theorems relies on the convergence of the feasible region of the linear program (ref) and the continuity of the objective of the linear programs with respect to the feasible region. The proof of both theorems is provided in Appendix (ref).
We now extend the framework to causal inference, where a platform conducts a randomized experiment to evaluate a new service feature, and customer ratings may be MNAR in both arms. The goal is to estimate, or bound, the average treatment effect. Partial identification is especially useful here, since decision-making often requires only the sign of the treatment effect rather than a point estimate. For notational simplicity, we suppress $X$ throughout the exposition; all results extend directly to the covariate-adjusted case by conditioning on $X$ and then aggregating over the distribution of $X$.
Consider a randomized experiment where each unit $i$ has potential outcomes $\{Y_i(1), Y_i(0), R_i(1), \allowbreak R_i(0)\}$ drawn i.i.d.\ from a distribution $\mathcal{D}$, with $Y_i(d) \in [M]$ denoting the discrete rating outcome and $R_i(d) \in \{0,1\}$ the observation indicator under treatment $d \in \{0,1\}$. Under treatment assignment $D_i$, the observed data are $O_i = (D_i, R_i(D_i), R_i(D_i) Y_i(D_i))$, and randomization ensures $D_i \indep (Y_i(1), Y_i(0), R_i(1), R_i(0))$. Our target is the average treatment effect \[ \tau = \E[Y(1) - Y(0)] \,. \]
We can bound $\tau$ by extending the LP-based identification strategy of Section (ref) to each arm. Specifically, we allow missingness to depend on the outcome within each arm: $\pi_d(y) = \P(R(d)=1 \mid Y(d)=y)$ remains unidentified under MNAR. Let $\alpha_d(y) = \P(Y(d)=y, R(d)=1)$ be the probability of observing rating $y$ in arm $d$, which is identifiable from the observed data.
Using the decomposition $\P(Y(d)=y) = \alpha_d(y)/\pi_d(y)$, together with the constraint $\sum_{y=1}^M \P(Y(d) = y) = 1$, the feasible region for each arm is $$\Pi_d = \left\{ (\pi_d(1), \ldots, \pi_d(M)): \sum_{y=1}^M \alpha_d(y)/\pi_d(y) = 1,\; \pi_d(y) \in (0,1] \right\}\,.$$
The identification region for $\tau$, denoted by $\mathcal{T}$, is obtained by optimizing over all feasible missingness mechanisms in both arms. Letting $w_d(y) = 1/\pi_d(y) - 1$ and $\mathcal{T}=[\tau_{\min}, \tau_{\max}]$. Then $\tau_{\min}$ and $\tau_{\max}$ are the optimal values of the following pair of linear programs:
Here, “${\max/\min}$” refers to computing the upper and lower bounds, with “$\tau_{\max/\min}$” denoting the corresponding identified bounds, respectively.
Although $\tau$ is not point-identified, reliable decisions may still be possible if the entire identification region lies on one side of zero. This allows the platform to determine the sign of the treatment effect and make deployment decisions without requiring point identification. We provide two sufficient conditions under which the treatment effect is guaranteed to be nonnegative. Notably, both conditions can be tested using observed data.
This condition has an intuitive interpretation: the left-hand side represents the difference in observed outcome means weighted by the joint observation probabilities, while the right-hand side captures the worst-case scenario where all missing control outcomes equal the maximum score $M$. When the observed advantage of treatment exceeds this worst-case penalty, we can conclude a positive treatment effect regardless of the true missingness mechanism.
This condition holds when the missingness mechanism is identical across treatment arms---a plausible assumption when treatment does not affect the propensity to respond. Under this constraint, if the treatment shifts probability mass from lower to higher ratings in a near monotone fashion, then the treatment effect is guaranteed to be positive. This is a much weaker condition than Proposition (ref) as we do not require the magnitude of the differences. The proofs of both propositions are provided in Appendix (ref).
We now introduce a weak shadow variable $F_i(d) \in \mathcal{F}$ representing, for instance, an LLM-generated prediction of the customer's satisfaction based on the conversation transcript. The shadow variable is fully observed for all units regardless of whether they provide ratings. Under this scenario, the observed data become $O_i = (D_i, R_i(D_i), F_i(D_i), R_i(D_i) Y_i(D_i))$.
As in Assumption (ref), we assume the shadow variable satisfies conditional independence with respect to missingness:
This is the causal analogue of Assumption (ref): conditional on the true outcome, the externally generated prediction carries no additional information about the decision to respond. Under Assumption (ref), the identifiable quantities expand to include $\alpha_d(f,y) = \P(R(d)=1, F(d)=f, Y(d)=y)$ and $\beta_d(f) = \P(R(d)=0, F(d)=f)$. Following the same derivation as in Section (ref), we can express the constraints on the missingness mechanism as: \[\beta_d(f) = \sum_{y=1}^M \alpha_d(f,y) \left(\frac{1}{\pi_d(y)} - 1\right), \quad \forall f \in [M].\]
Letting $w_d(y) = 1/\pi_d(y) - 1$ and defining matrices $A_d = [\alpha_d(f,y)]_{f,y} \in \R^{|\mathcal{F}| \times M}$, vectors $\mathbf{w}_d = (w_d(1), \ldots, w_d(M))^\top$, $\bm{\beta}_d = (\beta_d(1), \ldots, \beta_d(M))^\top$, and $D = \mathrm{diag}\{1, \ldots, M\}$, the identification region is characterized by:
The key difference from the formulation without shadow variables is that we now have $|\mathcal{F}|$ linear constraints per arm (one for each value of $F$), compared to a single aggregated constraint. This additional structure tightens the feasible region and consequently the identification region.
The set expansion estimator from Section (ref) extends directly to this setting. The following proposition shows that the resulting estimators attain the same consistency rates as Theorems (ref) and (ref), with the effective sample size $\underline{n} = \min(n_0, n_1)$, where $n_0$ and $n_1$ are the number of control and treated observations, respectively.
We evaluate the proposed identification bounds and set-expansion estimator on both simulated data and real-world customer service dialogue data. In both settings, the goal is to estimate the mean of an outcome variable that is MNAR, with an auxiliary prediction serving as a shadow variable. In the first simulation data, the prediction is synthetically generated, while in the second experiment, the prediction is generated through a Large Language Model (LLM) by prompting the actual dialogue.
We compare against several canonical baseline methods representing three distinct approaches: naive handling of missing data, classical MNAR models with parametric assumptions, and more recent prediction-powered inference that assumes MAR. Let $N$ denote the total number of units, $N_1 = \sum_{i=1}^N R_i$ the number of observed outcomes, and $N_0 = N - N_1$ the number of missing outcomes. The baseline estimators are defined as follows:
We simulate a platform rating setting with $M = 5$ and true mean $\mu^* = 3.0$. The outcome distribution is a discretized normal centered at $\mu^*$. For the construction of shadow variables, we examine two configurations of $\mathbb{P}(F \mid Y)$ that yield qualitatively different identification regimes:
For the missingness mechanism, we consider a MNAR pattern motivated by empirical evidence that customers with more extreme satisfaction levels (either lowest or highest) are more likely to leave reviews on service platforms, e.g., hu2017self. Specifically, we set the response probabilities to $\pi(1) = 0.30$, $\pi(2) = 0.10$, $\pi(3) = 0.05$, $\pi(4) = 0.70$, and $\pi(5) = 0.95$. The resulting overall response rate is approximately 36.6%, and complete-case analysis yields substantial upward bias as response rates are higher for higher scores.
The identification region for the mean outcome in both cases without shadow variable is $[2.15, 4.64]$ as they enjoy the same data distribution without shadow variables. Under the point identification configuration, the oracle identification region with the shadow variable collapses to the singleton $\{3.0\}$, while under the set identification configuration, the shadow variable tightens the bounds to $[2.78, 3.20]$, representing a reduction in width from 2.49 to 0.42. This demonstrates the strong identification power of the shadow variable under this setup.
We compare with all five baseline methods mentioned earlier with sample sizes ranging from 1,000 to 10,000. We average results over 100 runs and present the point estimators along with their confidence intervals in Figure (ref). We observe that in both cases, the set expansion estimator successfully covers the true value, while the other estimators all hold a consistent upward bias. This is because the prediction holds a general upward bias that is passed to the final estimator due to naive handling of the MNAR pattern (CCA, PPI, and NI) or misspecified parametric model (Heckman and PM). Moreover, the confidence intervals for those baseline estimators converge outside of the identification region, meaning a consistently higher estimate. Additionally, in the set-identified case, both the upper bound and lower bound from the set expansion estimator converge to their oracle value, confirming the theoretical guarantee of Theorems (ref) and (ref).
Next, we validate our estimator using a semi-synthetic experiment based on real dialogue satisfaction data from e-commerce customer service interactions. We simulate for a wide range of missingness patterns to test the robustness and efficiency of the proposed approach.
We use the public User Satisfaction Simulation (USS) dataset sun2021simulating, a collection of human-annotated dialogues compiled from five public dialogue corpora. In our experiment, we focus on one of the corpora called JD.com, which is the second largest online retailer in China. The dataset consists of 3,300 Chinese customer service dialogues, each containing multiple conversation turns between users and automated systems. Two examples of the dialogues can be found in Appendix (ref). Independent human annotators provide overall satisfaction ratings on a 1--5 scale (1 = very dissatisfied, 5 = very satisfied) for each dialogue. Complete distribution can be found in Table (ref), which shows a concentrated distribution with medium scores (3 or 4).
Next, we generate the LLM-predicted satisfaction rating by prompting GPT-4o-mini with the dialogue transcript and asking for a satisfaction rating. While we use GPT-4o-mini in our main experiments, the general insights carry over to other pretrained foundation models. We consider four prompting strategies:
Table (ref) summarizes the predictive accuracy across four prompting strategies. The Mean Absolute Errors (MAEs) range from 0.476 to 0.607, indicating modest prediction accuracy. We also report the minimal singular value of the joint distribution matrix $H := (\P(F_i=f, Y_i=y))_{f,y \in [5]}$, which quantifies the completeness condition required for point identification. The minimal singular values are on the order of $10^{-4}$ to $10^{-5}$ and the condition numbers are at the scale of $10^4$, indicating that the completeness condition barely holds and point identification is fragile.
Since the USS dataset contains complete outcome data, we simulate MNAR missingness by selectively masking outcomes according to outcome-dependent response probabilities $\pi(y) = \P(R = 1 \mid Y = y)$. This semi-synthetic design allows us to evaluate estimator performance against known ground truth. We generate 1,000 random missingness patterns where each $\pi(y)$ is drawn independently from $\text{Uniform}(0.1, 0.9)$, and additionally examine three representative patterns reflecting realistic nonresponse mechanisms on service platforms:
The exact missingness probability can be found in Figure (ref).
We compare two interval estimators: Set Expansion (with $\kappa_n = 0.5$ and $C = 50$) and Aggregated LP (which does not use the shadow variable), and six point estimators with five mentioned above and an additional LLM Raw estimator that directly outputs the average of the LLM prediction. All point estimators assume some form of MAR, either unconditionally or conditional on $F$, which fails under our MNAR data generating process. Note that when calculating MAE, we use the midpoint of the set expansion output interval as the point estimate.
Figure (ref) displays results under the three representative missingness patterns using CoT prompting with missingness probability on the top left of each figure. The true population mean is $\mu = 3.73$. The Set Expansion interval (blue region) covers the true mean in all three cases, while the Aggregated LP interval (green region) is 3--5 times wider, confirming that the shadow variable substantially improves estimation. All point estimators exhibit systematic bias under at least one pattern: CCA overestimates when lower scores are missing and underestimates when higher scores are missing, while NI and PPI inherit bias from the violated MAR assumption. The bias is largest when the missingness pattern opposes the LLM prediction bias. In our data, LLM predictions are biased downward ($\bar{F} < \bar{Y}$), so Panel (a)---where higher scores are selectively missing---produces the largest downward bias in point estimators. Conversely, when missingness and prediction bias align as in Panel (c), point estimators have smaller errors. Crucially, the proposed Set Expansion estimator correctly covers the true mean regardless of the missingness pattern.
Table (ref) summarizes performance across all four prompting strategies over 1,000 random missingness patterns. The Set Expansion estimator achieves the lowest MAE in all settings, outperforming the best point estimator by a factor of 2--3 while maintaining coverage above 98%. The fine-tuned model shows higher MAE (0.045) than CoT (0.031), reflecting both the smaller held-out sample size and out-of-sample generalization challenges. Naive Imputation exhibits highly variable performance---excellent with few-shot prompting (MAE = 0.059) but poor otherwise---highlighting its sensitivity to LLM calibration. The Heckman two-step estimator performs poorly (MAE $\approx$ 1.0) because it requires valid exclusion restrictions---variables affecting selection but not the outcome---which are unavailable when only $F$ is observed. Pattern-Mixture, which simply assumes MAR conditional on $F$, performs substantially better (MAE $\approx$ 0.09). The average interval width for Set Expansion ranges from 0.25 to 0.37 across prompting strategies, compared to 1.46 for Aggregated LP, representing 75--83% width reduction.
These results validate three key claims: (i) the shadow variable framework reduces interval width by over 80% compared to methods without auxiliary predictions; (ii) the set expansion estimator maintains valid coverage across diverse MNAR patterns where point estimators fail; and (iii) the method is robust to LLM prediction quality, with even zero-shot prompting yielding substantial improvements.
We study the problem of estimating population quantities when outcomes are missing not at random, a pervasive challenge on digital platforms and in social surveys. Rather than imposing strong parametric assumptions or seeking point identification under potentially fragile conditions, we adopt a partial identification perspective and show that sharp bounds on the mean can be computed via a pair of linear programs. Our key insight is that predictions from pretrained models---including large language models---can be incorporated as weak shadow variables to tighten these bounds. Unlike classical shadow-variable approaches, our framework does not require completeness or strong predictive accuracy; it extracts useful information from imperfect predictions while providing valid coverage guarantees through a set-expansion estimator. In simulations and semi-synthetic experiments on customer-service dialogues, even simple LLM predictions reduce identification intervals by 75--83% and maintain coverage above 98% across diverse missingness patterns.
Our work opens several promising directions for future research. One is to leverage richer model outputs, e.g., multiple prompts and ensembles, to construct auxiliary signals that tighten identification regions. Another is to study data collection and incentive design, e.g., how to design elicitation mechanisms to minimize the bound width under budget constraints and strategic behavior.
\setcounter{page}{1}