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.
45,657 characters · 12 sections · 46 citation commands
Automatic debiased machine learning and sensitivity analysis for sample selection models
\thispagestyle{empty}
In many empirical studies, researchers face the challenge that outcomes are only observed for a subset of the sample population. Returns to education studies observe wages only for employed individuals. Job training evaluations miss earnings data for unemployed participants. Clinical trials lose patients before outcome measurement and also educational interventions suffer when students do not take standardized tests. This phenomenon, commonly referred to as sample selection or outcome attrition, complicates the estimation of causal effects heckman1976common,heckman1979sample,hausman1979attrition,little1995modeling. The problem becomes even more complex when treatment assignment is itself non-random. In such cases, researchers confront what Bia2024 describe as the “double selection problem”, involving both selection into treatment and selection into outcome observability. Standard methods for confounding adjustment, such as regression or propensity score weighting, fail when outcomes are selectively missing. Even inverse probability weighting, which addresses treatment selection, requires modification to handle missing outcomes robins1994estimation,hernan2004structural. The machine learning literature offers powerful tools for high-dimensional covariate adjustment, yet it also introduces new challenges. In particular, regularization inherent in machine learning estimators can induce bias that invalidates standard inference procedures unless appropriate orthogonality conditions are imposed chernozhukov2018double. Bia2024 address this issue by deriving a Neyman-orthogonal score function for treatment effect estimation in the presence of sample selection. dolgikh2025double also propose double machine learning estimators for treatment effect estimation in the multivariate sample selection model with ordinal selection equations. An alternative approach is provided by the Riesz representation theorem. Instead of relying on Neyman-orthogonal score functions, target parameters can be characterized through unique weighting functions called Riesz representers chernozhukov2022debiased. The Riesz framework offers several advantages: it avoids unstable propensity score inversions, enables direct estimation via variational or adversarial methods, and naturally accommodates sensitivity analysis. Existing applications include the estimation of the Average Treatment Effect (ATE) and other policy-relevant causal parameters in settings without selection chernozhukov2022debiased and the estimation of the Average Treatment Effect on the Treated (ATT) in Difference-in-Differences models bach2025sensitivityanalysistreatmenteffects.
This paper extends the Riesz representation methods to sample selection models. We show that the ATE identified by Bia2024 via efficient scores also admits identification through a Riesz representation. The corresponding representer takes the form of inverse probability weights that adjust simultaneously for treatment assignment and sample selection. The Riesz representer framework is particularly useful for analyzing bias induced by unobserved selection confounding. It yields an interpretable decomposition of the omitted variable bias. Building on cinelli2020making and chernozhukov2022long, we express the bias as the product of three terms: (i) a scale factor identifiable from observed data, (ii) the strength of confounding in the outcome equation, and (iii) the strength of confounding in the selection equation. This decomposition delivers sharp bounds on the magnitude of bias without requiring the specification of the full joint distribution of unobservables. A key insight is that observed covariates provide natural benchmarks for calibrating these sensitivity parameters imbens2003sensitivity,altonji2005selection,oster2019unobservable. In our simulation study, we investigate the finite-sample behavior of the proposed ForestRiesz estimator and find that it performs well in finite samples when estimating the ATE. As an empirical contribution, we study the gender wage gap in the U.\,S. using data from the American Community Survey. We find that our ForestRiesz approach yields larger treatment effect estimates than a standard double machine learning approach which does not account for sample selection. This suggests that ignoring sample selection leads to an underestimation of the gender wage gap, as wage reporting behavior differs systematically between female and male respondents.
Estimation of treatment effects is fundamental to empirical research in economics, medicine, and the social sciences. This section introduces the Average Treatment Effect (ATE) within the potential outcomes framework and examines the problem of sample selection, which occurs when outcome data are missing for some units in the analysis.
To formally define causal effects, we rely on the Potential Outcomes framework, often associated with rubin1974estimating, rubin1977assignment. Let $D$ be a variable that represents the treatment status assigned to an individual unit $i$. For clarity, we consider a binary treatment where $D_i=1$ if unit $i$ receives the treatment and $D_i=0$ if unit $i$ receives the control, though the framework readily extends to multiple discrete treatments $d \in \{0, 1, ..., Q\}$. For each unit \(i\), we define two potential outcomes: \(Y_i(1)\) is the outcome that unit \(i\) would have experienced under treatment (\(D_i = 1\)), while \(Y_i(0)\) is the outcome it would have experienced under control (\(D_i = 0\)). We assume SUTVA: unit $i$'s potential outcomes are unaffected by other units' treatment assignments, and $Y_i(1)$ and $Y_i(0)$ are well-defined for each unit rubin1980randomization. Given the impossibility of observing individual treatment effects directly for a single unit, empirical research typically focuses on estimating average causal effects across a population or subpopulation. The most common target parameter is the Average Treatment Effect (ATE) for the entire population: \[ \text{ATE} = \mathbb{E}[Y(1) - Y(0)] = \mathbb{E}[Y(1)] - \mathbb{E}[Y(0)], \] where the expectation $\mathbb{E}[\cdot]$ is taken over the distribution of units in the population of interest.
In many practical applications, the outcome variable $Y$ is not observed for all units in the sample. This issue is known as sample selection, outcome attrition, or nonresponse. Let $S$ be a binary indicator variable such that $S_i=1$ if the outcome $Y_i$ is observed for unit $i$, and $S_i=0$ otherwise. If the mechanism determining whether the outcome is observed ($S=1$) is related to the potential outcomes $Y(d)$ themselves, even after conditioning on treatment status $D$ and covariates $X$, then the subsample for whom we observe the outcome is no longer representative of the full population concerning the outcome process. Simply performing an analysis on the selected sample without accounting for the selection mechanism can introduce sample selection bias. When non-random treatment assignment (violating $Y(d) \perp D$) occurs simultaneously with non-random sample selection (violating $Y(d) \perp S$), researchers face a so-called double selection problem, as emphasized by Bia2024. In this situation, valid estimation requires assumptions addressing both sources of potential bias. The first assumption invokes conditional independence of the treatment:
By Assumption (ref), no unobservables jointly affect the treatment and the potential outcomes conditional on covariates $X$. Analogous to how Assumption (ref) addresses confounding in treatment assignment, specific assumptions are required to handle sample selection.
A common starting point is another selection-on-observables assumption, but applied to the selection process $S$. This is often termed the Missing-At-Random (MAR) assumption rubin1976inference or conditional independence of selection. In the context of treatment evaluation, it means that conditional on observed variables (importantly, treatment status $D$ and covariates $X$), the selection indicator $S$ is independent of the potential outcomes:
This conditional independence of selection assumption states that, within groups defined by a specific treatment status $d$ and covariate values $x$, whether an outcome $Y(d)$ is observed ($S=1$) or missing ($S=0$) does not depend on the potential outcome's value itself. Selection is allowed to depend on treatment $D$ and covariates $X$, but not on any unobserved factors related to $Y(d)$ once $D$ and $X$ are accounted for. However, this might be violated in many real-world scenarios. Selection could depend on unobserved factors (denoted $A$) that also influence the potential outcome, even after conditioning on $D$ and $X$. This is known as non-ignorable nonresponse or selection based on unobservables. In the following, we consider a much weaker condition, i.\,e., selection independence only holds if we could condition on the additional unobserved factors $A$. This leads to the following assumption:
Under Assumption (ref), selection is independent of potential outcomes once we account for treatment status, observed covariates, and the unobserved selection confounding factors $A$. While we cannot observe $A$ directly, we can apply the framework of chernozhukov2022long to provide sharp bounds on the size of the omitted variable bias that results from not observing $A$. When introducing unobserved confounders $A$ in the selection process $S$, we need to specify how this affects the treatment assignment as well. We impose the following assumption:
It is worth noting that this assumption is implied by $X$ being the only causal parent of $D$, e.\,g., in a stratified randomized control trial. The analysis allowing unobserved confounders $A$ in both the sample selection and the treatment assignment will be pursued in an extended version of this work.
We also make the following assumption. Let $p_d(X):=\mathbb{P}(D=d\mid X)$ for $d\in\{0,1\}$ and let $\pi_0(d,X,A):=\mathbb{P}(S=1\mid D=d,X,A)$.
We refer to Equation (ref) as a weak overlap condition since it requires only integrability of inverse propensities (rather than uniform lower bounds). The first part of the assumption is a conventional common support condition, which ensures that treatment assignment is non-degenerate and the probability of selection is always non-zero for each conditioning value.
Further, we denote the conditional mean outcome by $\mu_d(X) = \mathbb{E}[Y|D=d, S=1, X] $. Under Assumption (ref), Assumption (ref), and conditional independence of selection in Equation (ref), the ATE is identified by: $$ \theta_0 = \mathbb{E}[\phi_1 - \phi_0] $$ with
being the efficient score function derived by Bia2024. Hence, the ATE is identified using outcomes $Y$ from the selected sample ($S=1$) and selection indicators $S$ for all units. Intuitively, identification involves modeling the conditional outcome mean within the selected sample, $\mathbb{E}[Y|D=d, S=1, X]$, and then appropriately adjusting or re-weighting based on estimates of the treatment propensity score $p_d(X) = \mathbb{P}(D=d|X)$ and the selection propensity score $\pi_s(d, X) = \mathbb{P}(S=1| D=d, X)$.
Many empirical problems now involve rich covariates. Machine-learning methods like Lasso, random forests, and neural networks can estimate nuisance functions such as conditional means and propensities in these settings. They achieve good prediction through regularization and model selection. These devices, however, typically introduce bias. If we plug a regularized estimate $\hat g$ into a target functional, the resulting estimator can inherit non-negligible bias and invalidate $\sqrt{n}$-consistent inference.
Debiased machine learning (DML) addresses this problem by using Neyman-orthogonal scores Levit1975,IbragimovHasminskii1979,chernozhukov2018double. In this framework, a score $\psi(W,\theta,g)$ identifies $\theta_0$ via the moment condition $$\mathbb{E}[\psi(W,\theta_0,g_0)]=0$$ and is constructed so that small errors in $g$ have only a second-order effect on the moment, where $W$ denotes the data. Bia2024 derive a Neyman-orthogonal score for high-dimensional sample selection models (see Equation (ref)) and use cross-fitting to obtain valid inference.
A complementary approach uses the Riesz representer. For many parameters of interest (including the ATE), we can write \[ \theta_0=\mathbb{E}[m(W,g_0)], \] where the map $g\mapsto \mathbb{E}[m(W,g)]$ is linear and continuous on a suitable function class. The Riesz Representation Theorem then yields a unique function $\alpha_0$, called the Riesz representer, such that \[ \mathbb{E}[m(W,g)] = \mathbb{E}[\alpha_0(Z)\,g(Z)] \] for all admissible functions $g$, where $Z$ collects the arguments of $g$ and $\alpha_0$.
The condition $\mathbb{E}[\alpha_0(Z)^2]<\infty$ is closely linked to $\theta_0$ having a finite semiparametric efficiency bound newey1994asymptotic, hirshberg2021augmented, chernozhukov2022automatic. In our setting, the efficient score of Bia2024 admits an analogous Riesz representation that combines treatment and selection propensity weights. The representer also leads to a generic orthogonal score. For target parameters of the form $\theta_0=\mathbb{E}[m(W,g_0)]$ with $g_0(Z)=\mathbb{E}[Y\mid Z]$, consider
where $g$ and $\alpha$ approximate $g_0$ and $\alpha_0$. It is worth noting that in our sample selection model, $Y$ is only observed when $S=1$, so we basically consider $Y=SY$. As shown by chernozhukov2022automatic, evaluating at the true $\theta_0$ yields \[ \mathbb{E}[\psi(W,\theta_0,g,\alpha)] = -\,\mathbb{E}\bigl[(\alpha(Z)-\alpha_0(Z))(g(Z)-g_0(Z))\bigr]. \] Thus, the score is doubly robust: its expectation is zero if either $g=g_0$ or $\alpha=\alpha_0$, and errors enter only through their product. Combined with cross-fitting, this property delivers $\sqrt{n}$-consistent inference with flexible first stages chernozhukov2018double.
The Riesz formulation also plays an important role for estimation and sensitivity analysis. It casts the problem as learning a weighting function $\alpha_0$ jointly with $g_0$, which aligns well with variational, adversarial, and forest-based methods and can improve numerical stability and transparency. By learning the Riesz representer directly rather than relying on plug-in inverse probability weights, this approach can reduce instability when estimated propensities or selection probabilities are small. Most crucially for our case, the same representer-based structure naturally supports the sensitivity analysis in Section (ref).
Our goal is to identify the Average Treatment Effect (ATE), $\theta_0 = E[Y(1) - Y(0)]$, in the sample selection model described in Section (ref), where non-random treatment assignment occurs simultaneously with non-random sample selection. Under Assumptions (ref)--(ref), the ATE admits the following representation in the long model (i.\,e., in a hypothetical setting where the latent factors $A$ were observed):
where $W:=(Y,D,S,X,A)$ is the so-called long data vector and $g_0(d,x,a):=\mathbb{E}[Y\mid D=d,S=1,X=x,A=a]$ is the long regression. Since $A$ is not observed in practice, we are only able to identify the so-called “short" parameter $$\theta_s = \mathbb{E}[m(W_s, g_s)]= \mathbb{E}[g_s(1, X) - g_s(0, X)]$$ from the observed short data vector $W_s:=(Y, D, S, X)$, where $g_s(d,X) = \mathbb{E}[Y| D=d, S=1, X]$ is the short regression. Since both parameters have a representation of the form $\theta=\mathbb{E}[m(W, g)]$, the Riesz Representation Theorem guarantees the existence of a Riesz representer $\alpha$, such that $\theta=\mathbb{E}[\alpha(Z) g(Z)]$. The following main theorem of this paper, provides the explicit form of the Riesz representer in sample selection models.
The formal proof is given in Appendix (ref). Intuitively, the Riesz representer reweights the data to mirror what we would see in a randomized experiment. Weighting by $1/\mathbb{P}(D=d \mid X)$ increases the influence of units with observed characteristics $X$ that are unlikely to receive treatment $d$. This reweighting aligns the distribution of observed confounders, mimicking the balance achieved through random assignment. Introducing sample selection creates an additional challenge: outcomes are observed only when $S=1$. To correct for this, we apply a second set of inverse-probability weights based on the likelihood of selection. Since we do not rely on conditional independence of selection in Equation (ref) but rather on Assumption (ref), conditional independence of selection holds only after controlling for the unobserved variables $A$. As a result, our correction for selection must also account for these unobservables. Weighting by $1/\pi_0(d,X,A)$ in the long parameter, or $1/\pi_s(d,X)$ in the short parameter, gives more weight to units that were less likely to be selected into the observed sample, thereby restoring representativeness relative to the full population. The distinction between long and short parameters reflects whether the weighting scheme accounts for the unobserved confounders $A$ in the selection process $S$ or not.
With observed data we are only able to identify the short parameter, although we are interested in the long parameter $\theta_0$. The Riesz representer theorem gives us a direct formula for the omitted variable bias arising from not controlling for $A$ in the selection into observability. Following chernozhukov2022long, the difference between the long parameter $\theta_0$ and the short parameter $\theta_s$ is given by $$\theta_0 - \theta_s = \mathbb{E}[(g_0 - g_s)(\alpha_0 - \alpha_s)],$$ which can be interpreted as the covariance between the error parts of $g$ and $\alpha$. Therefore, the (squared) bias is bounded by $$|\theta_0 - \theta_s|^2 = \rho^2B^2 \leq B^2,$$ where $B^2:= \mathbb{E}[(g_0-g_s)^2]\mathbb{E}[(\alpha_0-\alpha_s)^2]$ and $\rho^2 := \text{Cor}^2(g_0-g_s,\alpha_0-\alpha_s)$. Furthermore, this squared bias bound \( B^2 \) has an intuitive decomposition that helps to understand the role of confounding in sample selection models. The squared bias bound \( B^2 \) can be decomposed as \[ B^2 = \widetilde{S}^2 C^2_Y C^2_S, \] where $\widetilde{S}^2:= \mathbb{E}[(Y - g_s)^2] \mathbb{E}[\alpha_s^2]$, $C^2_Y := \frac{\mathbb{E}[(g_0 - g_s)^2]}{\mathbb{E}[(Y - g_s)^2]}$ and $C^2_S := \frac{\mathbb{E}[(\alpha_0 - \alpha_s)^2]}{\mathbb{E}[\alpha_s^2]}$. Therefore, the bound \( B^2 \) is the product of \( \widetilde{S}^2 \), a scaling factor identifiable from observed data, \( C^2_Y \) that measures confounding strength in the outcome equation and \( C^2_S \) that measures confounding strength in the selection equation. For \( C^2_Y \) and \( C^2_S \) researchers need to make informed assumptions about the impact of unobserved confounding. More formally, it holds that
which measures the proportion of residual outcome variation (variation not explained by observed variables) that can be explained by the latent confounders $A$. It is by definition $\eta^2_{Y \sim A \mid D,X,S=1}$, the partial $R^2$ of $Y$ on the confounder $A$, after adjusting for $D$ and $X$, conditional on $S=1$. Further, it holds \[ C_S^2=\frac{\mathbb{E}[\alpha_0^2]-\mathbb{E}[\alpha_s^2]}{\mathbb{E}[\alpha_s^2]} =\frac{1-R^2_{\alpha_0\sim\alpha_s}}{R^2_{\alpha_0\sim\alpha_s}}, \quad \text{with} \quad R^2_{\alpha_0\sim\alpha_s}=\frac{\mathbb{E}[\alpha_s^2]}{\mathbb{E}[\alpha_0^2]}. \] It is worth noting that \( R^2_{\alpha_0 \sim \alpha_s} = \mathbb{E}[\alpha_s^2]/\mathbb{E}[\alpha_0^2] \) measures how much variation in the true Riesz representer \( \alpha_0 \) is explained by the short Riesz representer \( \alpha_s \). Therefore, \( 1 - R^2_{\alpha_0 \sim \alpha_s} \) (bounded between 0 and 1) measures the proportion of variation in \( \alpha_0 \) that is explained by the omitted confounder \( A \). While this parameter also admits an interpretation as a gain in precision, we find it more informative to use the following quasi-Gaussian approach for interpretation: \\ Quasi-Gaussian Selection Sensitivity. In practical applications, it might be difficult to think of plausible values for \( 1 - R^2_{\alpha_0 \sim \alpha_s}\), a technical and likely unfamiliar parameter. Instead, we find it useful to represent the selection indicator $S$ in a form of a latent index $S^*$ with Gaussian shocks crossing a threshold: Let $S=\mathbf{1}\{S^*>0\}$ with $$ S^* = h(D,X) - U \quad \text{and} \quad U \mid D,X \sim N(0,1).$$ This representation does not entail loss of generality. We can then model confounding as follows: $$ U = \mu_S A + \sqrt{1-\mu_S^2}\,\varepsilon_S, \quad \text{with} \quad A,\varepsilon_S \stackrel{\mathrm{i.i.d.}}{\sim} N(0,1), $$ independent of $(D,X)$. Thus, $\mu^2_S$ is the $R^2$ in the regression of the Gaussian shock $U$ on the latent confounder $A$. By definition, it is also equal to $\eta^2_{S^* \sim A | D, X}$, the nonparametric partial $R^2$ in the regression of the latent index $S^*$ on $A$, after nonparametrically partialling out $(D,X)$. It is therefore easy to interpret. We can also map $\mu^2_S$ to the technical sensitivity parameter as follows. We compute the short selection probability $\pi_s(d,x)=\mathbb{P}(S=1\mid D=d,X=x)=\Phi(h(d,x))$, so $h(d,x)=\Phi^{-1}(\pi_s(d,x))$ is identified from the short model, and the long probability is $$\pi_0(d,x,a)=\mathbb{P}(S=1\mid D=d,X=x,A=a)=\Phi\left((h(d,x)-\mu_S a)/\sqrt{1-\mu_S^2}\right).$$ We show in Appendix (ref) that $\mathbb{E}[\alpha_0^2]$ and $\mathbb{E}[\alpha_s^2]$ can be expressed in terms of these probabilities and can therefore be seen as functions of $\mu^2_S$. We then derive the maps from the interpretable to the technical sensitivity parameters: $\mu^2_S \mapsto 1-R^2_{\alpha_0\sim\alpha_s} (\mu^2_S).$ This yields a one-parameter, probit-scale calibration of selection confounding that is directly compatible with the Riesz-based bias bounds. Note that this does not impose any assumptions on the data, but is rather an interpretation device. While sensitivity analysis maps assumptions about the unobserved confounder \(A\) (which might affect both outcome \(Y\) and selection \(S\)) to potential bias in the ATE estimate \(\theta_s\), it does not tell us how plausible those assumptions are. Researchers must therefore make informed judgments about the two partial $R^2$ measures that capture how strongly $A$ predicts the outcome $Y$ and the selection index $S^*$. This task can be aided by a benchmarking approach, following imbens2003sensitivity, altonji2005selection, oster2019unobservable, cinelli2020making, and chernozhukov2022long, which uses the observed influence of specific covariates \(X_j\) as a reference point for the potential influence of an unobserved confounder $A$. We outline this approach in Appendix (ref).
Since \(\alpha_0\) is generally unknown, constructing a feasible estimator based on the orthogonal score (Equation (ref)) in the DML framework requires an estimate \(\hat{\alpha}\). The traditional method for obtaining \(\hat{\alpha}\) is a plug-in approach. While conceptually straightforward, this plug-in approach for estimating the Riesz representer suffers from several drawbacks, particularly in high-dimensional or complex settings. Deriving the analytical form of \(\alpha_0\) can be mathematically challenging or even intractable for more complex parameters of interest beyond the standard ATE. The formula for \(\alpha_0\) also frequently involves division by estimated probabilities or densities (see, e.\,g., $\hat{p}_d(X)$ and $\hat\pi_s(d,X)$ in Theorem (ref)). If these estimated quantities are close to zero, the resulting \(\hat{\alpha}\) can become extremely large. This occurs when the common support assumption (positivity) is empirically violated in the sample. Such large values can lead to unstable estimates of the target parameter \(\theta_0\). Recognizing the limitations of the plug-in method, recent research has focused on methods that estimate the Riesz representer $\alpha_0$ directly, without needing its explicit analytical formula or relying on potentially unstable inverse weighting schemes. Two prominent direct approaches are variational methods (Riesz Regression) and adversarial (minimax) methods chernozhukov2020adversarial, chernozhukov2022automatic, chernozhukov2022riesznet. In this paper, we rely on the ForestRiesz, also developed by chernozhukov2022riesznet, that adapts the random forest methodology to estimate the Riesz representer. Within this framework, the Riesz representer is modeled as locally linear with respect to a pre-specified feature map $ a(Z) = \langle r(D, X, S), \beta(X) \rangle $, where $Z=(D,X,S)$, $r(D, X, S)$ represents a smooth feature map (e.\,g., a polynomial series) and $\beta(X)$ denotes local coefficients that vary with covariates $X$. The algorithm constrains splits to covariates $X$ exclusively to preserve sufficient variation in the treatment variable $D$ within each node. chernozhukov2022riesznet show that this problem falls in the class of problems defined via solutions to moment equations $m(\cdot)=0$. Therefore, we can apply the framework of Generalized Random Forests of 10.1214/18-AOS1709 to solve this local moment problem via random forests. For each node in the forest, the algorithm computes a Jacobian matrix and a local moment vector $$J(\text{node}) = \frac{1}{|\text{node}|}\sum_{i \in \text{node}} r(Z_i)r(Z_i)^\top \quad \text{and} \quad M(\text{node}) = \frac{1}{|\text{node}|}\sum_{i \in \text{node}} m(W_i;r).$$ The optimal coefficient vector within each node is given by $\beta(\text{node}) = J(\text{node})^{-1} M(\text{node}).$ ForestRiesz grows the forest by recursively splitting nodes based solely on the covariates \(X\). For each candidate split, the two resulting child nodes are evaluated by computing their respective \(J\) and \(M\). The splitting rule seeks to maximize the stability-adjusted signal by minimizing the aggregate local Riesz loss: \[ -\sum_{\text{child} \in \{1,2\}} |\text{child}| \; \beta(\text{child})^\top J(\text{child})\,\beta(\text{child}). \] This criterion favors splits that yield child nodes where the local moment \(M\) is both strong and well-supported by a diverse (i.\,e., well-spread) feature set, while penalizing splits that produce nodes with nearly singular \(J\). ForestRiesz incorporates multitasking capabilities, wherein the forest simultaneously learns the regression function $\hat{g}$ and the Riesz representer $\hat\alpha$ by augmenting the node-splitting criteria with regression-based objectives. The final estimate is given by
or, better yet, its cross-fitted form to avoid overfitting, leveraging Equation (ref) as proposed in chernozhukov2022riesznet, where $\mathbb{E}_n$ denotes the sample mean.
The finite-sample properties of the proposed ForestRiesz (FR) estimator are assessed with a simulation study. The data-generative process (DGP) follows the conditional missing-at-random (MAR) design outlined in Appendix E of Bia2024, with pre-treatment covariates $X$, a selection and treatment indicator $S, \thickspace D \in\{0,1\}$, error terms $u,\thickspace v,$ and $w$, and an outcome variable $Y$, that is only observed if $S = 1$:\\ $$ Y_i = \thickspace \theta_0 D_i \thickspace + X_i'\beta_0 + u_i, \quad S_i = \thickspace \mathbf{1}\{D_i \thickspace + \thickspace X_i'\beta_0 \thickspace + \thickspace v_i \thickspace > 0 \}, \quad D_i = \mathbf{1}\{X_i'\beta_0 \thickspace + \thickspace w_i > 0\}, $$\\[-1.16 ex] with $X_i \sim N(0,\sigma^2_X)$, $(u_i,v_i) \sim N(0,\sigma^2_{u,v})$, and $w_i \sim N(0,1)$. For MAR to hold, $\sigma^2_{u,v}$ is specified as an identity matrix, implying that conditional on the treatment indicator and covariates none of the unobservables jointly affect the selection and outcome equation.
In the DGP, we set the true ATE to $\theta_0 = 1$. To benchmark the performance of the ForestRiesz, we compare it to an interactive regression model (IRM) chernozhukov2018double, which does not adjust for the sample selection mechanism of the DGP, and to the sample selection model (SSM) by Bia2024, which uses efficient Neyman-orthogonal score functions within the DML framework to address sample selection. The benchmark estimators are implemented via the doubleML package Bach_DoubleML_-_Double, using random forests\,\footnote{For the exact specification of hyperparameters of the random forests and the DML parameters see Appendix (ref).} for estimating the nuisance functions and three-fold cross-fitting to prevent overfitting bias.
For the sample sizes $N \in \{1000,\thickspace 4000,\thickspace 16000\}$, Table (ref) reports each estimator's average results across $200$ Monte Carlo iterations. For each estimator and sample size, it presents the estimate (ATE), the standard error (SE), and the corresponding bias (MAE). Across all sample sizes, the IRM model underestimates $\theta_0$, since it does not account for sample selection. By contrast, both SSM and FR converge to the true $\theta_0=1$ when the number of observations increases. Moreover, as standard errors scale with $1/\sqrt{N}$, quadrupling the sample size reduces the standard errors of all estimators by approximately one half.
A more detailed comparison of the SSM and FR simulation results suggests a different bias-variance trade-off. Across all sample sizes, SSM yields smaller standard errors, whereas FR results indicate a faster decline in bias as the sample size increases. It is worth noting that the FR model is used without any tuning, while for the SSM we explored different random forest depths to improve propensity scores estimation and reduce bias. To complement the previously described considerations, Figure (ref) in Appendix (ref) presents the distribution of the ATE estimates across all Monte Carlo iterations. Furthermore, Appendix (ref) presents additional results for the SSM estimator, showing that under the Lasso specifications used to learn the nuisance parameters in the score of Bia2024, the SSM bias declines as expected given the linearity of the DGP.
These considerations highlight the importance of the choice of machine learning methods and hyperparameter tuning in the SSM approach, and more generally within the DML framework (pmlr-v236-bach24a), and demonstrate that the FR approach is considerably more robust. A more detailed empirical comparison between the DML-based methods and the Riesz representer approach is left for future research.
As an empirical application, we apply our method to estimate the gender wage gap in the U.\,S. We use data from the 2016 American Community Survey (ACS), which provides a representative 1\,% sample of the U.\,S. population under mandatory participation. Since some respondents do not report their wages, even though they are employed, any gender wage gap analysis based on the ACS data is subject to a sample selection problem. The dataset contains 158 variables for socio-economic characteristics at the individual and the household level, for example referring to education, industry, and occupation. We follow the study of bach2024heterogeneity and focus on two sub-populations in the ACS: respondents with a high school degree and those with a college degree. Our treatment variable $D$ is the gender of a respondent, with $D=1$ indicating a female respondent. Our outcome variable $Y$ denotes (log) weekly wages (in USD) and the indicator $S$ indicates whether $Y$ is observed (i.\,e., the respondent has reported her wage). In the high school sub-population, we have $372\ 728$ respondents and in the college sub-population $297\ 178$ individuals.
In order to estimate the gender wage gap, we apply the proposed ForestRiesz, where one fits a random forest that jointly learns the Riesz representer $\alpha$ and the regression function $g$ in one step as described in Section (ref). To demonstrate the relevance of our Riesz representer approach in sample selection models, we compare our estimation results with those obtained from the interactive regression model (IRM) and the SSM approach, both implemented using the doubleML package Bach_DoubleML_-_Double, as in the simulation study. We apply the three estimators to the high school and college subsamples and report point estimates, standard errors, and p-values. Table (ref) presents the estimation results for the college and the high school subsamples. For all three regression models, we find a significant gender wage gap in both subsamples, with a larger gap in the high school subsample than in the college subsample, in line with previous findings in bach2024heterogeneity. Since the estimated wage gap is approximately $3$ percentage points larger using the Riesz representer approach compared to IRM, our results suggest that we underestimate the gender wage gap when not controlling for non-reporting respondents. Applying a logit model to the reporting indicator $S$, we find that never-married female workers with a high university degree (professional degree) have a higher probability of reporting their income than their male counterparts, and that the relationship between experience and reporting also differs between men and women (see Table (ref) in Appendix (ref)). Because these covariates are also among the strongest predictors of wages (see Table (ref) in Appendix (ref)), estimates of the gender wage gap are subject to selection bias if these patterns are ignored. While the IRM model does not address this issue, both the ForestRiesz (FR) and the SSM approach correct for it by reweighting respondents with a lower probability of wage reporting.
Next, we conduct a sensitivity analysis to assess the robustness of our estimated treatment effects to unobserved confounding. Using observed covariates $X_j$ as benchmarks, this approach evaluates how influential an unobserved confounder $A$ would need to be to overturn our main findings. We perform this analysis for all covariates and report results for the six most influential covariate groups in the college subsample in Table (ref) in Appendix (ref). For each group $j$, the table reports the share of additional outcome variation $G_{Y,j}$, selection variation $G_{S,j}$, and their alignment measure $\rho_j$, detailed in Appendix (ref). Overall, the results indicate that the estimated gender wage gap in the college subsample is highly robust. Omitting the most influential covariate group, marital status, changes the ATE estimate by only $0.55$ percentage points. Notably, although education explains the largest share of variation in wages and in the Riesz representer (high $G_Y$ and $G_S$), it has virtually no effect on the estimated gender wage gap (low $\Delta \theta$), reflecting the weak correlation between the residual component of the outcome and Riesz representer models (small $|\rho|$). We further assess robustness through sensitivity analyses based on these benchmarks. First, we construct confidence intervals that account for unobserved confounding as strong as the marital status covariate. Figure (ref) in Appendix (ref) shows that even under this conservative scenario, the estimated ATE remains statistically significant. Second, we examine the magnitude of unobserved confounding required to overturn our conclusions. Figure (ref) in Appendix (ref) illustrates the potential bias as a function of $C^2_Y = \eta^2_{Y \sim A|D,X,S=1}$ and $\eta^2_{S^*\sim~ A |D, X}$, assuming the worst-case alignment ($\rho=1$). The robustness value (RV) for the college subsample is $0.063$, implying that an unobserved confounder would need to explain at least $6.3\,\%$ of both residual outcome and selection variation to nullify the estimated effect. This is substantially more than any observed covariate in our data can explain.
One main contribution of the paper is a bounds analysis for treatment effects when the traditional sample-selection model’s conditional missing-at-random (MAR) assumption fails. Although MAR is widely used, it is often hard to defend in applications. We relax MAR by introducing a latent confounder that affects selection and then derive the Riesz representer for the average treatment effect (ATE), which combines treatment-propensity weighting with selection-probability weighting. Using the resulting Riesz representers for the short and long models, we decompose the omitted-variable bias into three interpretable components. This decomposition yields sharp, distribution-free bounds on the magnitude of bias and provides a practical sensitivity-analysis toolkit for violations of MAR.
A second contribution is to adapt the ForestRiesz method of chernozhukov2022riesznet to treatment-effect estimation under sample selection. This automatic debiased machine learning approach jointly learns the outcome regression and the Riesz representer, avoiding the numerical instability of plug-in estimators that require direct inversion of estimated probabilities. Our simulations highlight the advantages of the ForestRiesz framework over more standard doubly robust plug-in approaches in finite samples. We illustrate the practical benefits of the method in an application to the U.\,S. gender wage gap using the American Community Survey. We find that ignoring sample selection leads to an underestimation of the wage gap, driven by systematic gender differences in wage reporting. A benchmarking-based sensitivity analysis indicates that this conclusion is robust.
Overall, our results highlight the importance of explicitly accounting for sample selection, particularly in survey-based studies, and demonstrate that the ForestRiesz estimator offers a robust, interpretable, and computationally attractive approach for causal inference in the presence of selective outcome observability.