EconBase
← Back to paper

Model Averaging and Double Machine Learning

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.

104,283 characters · 10 sections · 91 citation commands

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

Model Averaging and Double Machine Learning

\thispagestyle{empty}

abstract\singlespacing This paper discusses pairing double/debiased machine learning (DDML) with stacking, a model averaging method for combining multiple candidate learners, to estimate structural parameters. In addition to conventional stacking, we consider two stacking variants available for DDML: short-stacking exploits the cross-fitting step of DDML to substantially reduce the computational burden and pooled stacking enforces common stacking weights over cross-fitting folds. Using calibrated simulation studies and two applications estimating gender gaps in citations and wages, we show that DDML with stacking is more robust to partially unknown functional forms than common alternative approaches based on single pre-selected learners. We provide Stata and R software implementing our proposals. \noindentKeywords: causal inference, partially linear model, high-dimensional models, super learners, nonparametric estimation \\ JEL: C21, C26, C52, C55, J01, J08

\setcounter{page}{1}

Introduction

Motivated by their robustness to partially unknown functional forms, supervised machine learning estimators are increasingly leveraged for causal inference. For example, lasso-based approaches such as the post-double-selection lasso (PDS lasso) of Belloni2014a have become popular estimators of causal effects under conditional unconfoundedness in applied economics gilchrist2016something,dhar2022reshaping. Yet, a recent literature also raises practical concerns about the use of machine learning for causal inference. wuthrich2021 find that lasso often fails to select relevant confounders in small samples while inference based on linear regression performs relatively well. giannone2021economic and kolesar2023 argue that the sparsity assumption, on which the lasso fundamentally relies, is frequently not plausible in economic data sets. angrist2022 show that conditioning on confounders using random forests may yield spurious results in IV regressions.\footnote{See also Angrist2022b for additional discussion.} In an application to the evaluation of active labor market programs, Goller2020 find that random forests are not suitable for the estimation of propensity scores. A key characteristic shared by many of these studies using machine learning for causal inference is the focus on a single pre-selected machine learner.

This paper revisits the application of machine learning for causal inference in light of this recent literature. In particular, we highlight the benefits of pairing double/debiased machine learning (DDML) estimators of Chernozhukov2018 with stacking wolpert1996, Breiman1996a, laan2007. DDML can leverage generic machine learners meeting mild convergence rate requirements for the estimation of common (causal) parameters. Stacking is a form of model averaging that allows selecting among or combining multiple candidate machine learners relying on different regularization assumptions rather than requiring an {\it ad hoc} choice between them. Based on a diverse set of applications and calibrated simulation studies, we show that the synthesis of stacking and DDML improves the robustness of estimates of target parameters to the underlying structure of the data, and illustrate the finite sample performance of stacking-based DDML estimators. The results suggest that stacking with a rich set of candidate estimators can address some of the shortcomings highlighted in the recent literature on causal inference with single pre-selected machine learners.

We further consider two alternate ways of combining stacking and DDML aimed at improving practical feasibility and stability in finite samples: {\it Short-stacking} leverages the cross-fitting step of DDML to reduce the computational burden of stacking substantially. {\it Pooled stacking} decreases the variance of stacking-based learners across the DDML cross-fitting folds. Both approaches facilitate interpretability compared to conventional stacking by enforcing common stacking weights. We complement the paper with software packages for Stata and R that implement the proposed approaches Ahrens2023_ddml,ddml_R.

Model averaging techniques have a long tradition in economics and statistics. In the time-series literature, the idea of using `optimal' weights to combine forecasts goes back to the 1960s CraneCrotty1967,BatesGranger1969. Loss-minimizing combinations of a pre-specified set of estimators were introduced under the term stacking to the statistics literature by Wolpert1992 and Breiman1996a and generalized by laan2007. Stacking fits a final (typically parametric) learner on a set of cross-validated predicted values derived from distinct candidate learners. A popular choice is to constrain the weights attached to each candidate learner to be non-negative and sum to one. The benefits of combining multiple estimators into a `super learner' via stacking to improve robustness to the structure of the underlying data-generating process are well-known in the econometrics and statistics literature. Under appropriate restrictions on the data generating process and loss-function, laan2003 show asymptotic equivalence between stacking and the best-performing candidate learner.\footnote{Hastie2009 and van2011targeted provide textbook treatments of stacking and super learning. See also bhansen2012jma for discussion of jackknife (leave-one-out) stacking. clydec2013 and le2017 provide a Bayesian interpretation of stacking in the setting where the true model is not among the candidates.}

Model averaging methods and stacking are widely used in time-series forecasting and macro-econometrics steel2020,Wang2023review. Yet, despite its theoretical appeal, stacking has hitherto been rarely used for the estimation of causal effects in economics or other social sciences. An important exception is van2011targeted, who recommend stacking in the context of Targeted Maximum Likelihood. Instead, estimators are often based on parametric (frequently linear) specifications or single pre-selected machine learners. This can have severe consequences for the properties of causal effect estimators if the given choice is ill-suited for the application at hand. A simple example is shown in Figure (ref) which compares the performance of DDML using either cross-validated (CV) lasso or a feed-forward neural network to estimate a partially linear model across two different data-generating processes. The results show that the bias associated with each learner strongly depends on the structure of the data. Since true functional forms are often unknown in the social sciences, indiscriminate choices of machine learners in practice can thus result in poor estimates. DDML with stacking is a practical solution to this problem. As the example showcases, DDML using stacking is associated with low bias when considering a rich set of candidate learners that are individually most suitable to different structures of the data.

figure[figure omitted — 1,715 chars of source]

We conduct simulation studies calibrated to real economic datasets to demonstrate that stacking approaches can safeguard against ill-chosen or poorly tuned estimators in practical settings. Throughout, stacking estimators are associated with relatively low bias regardless of the simulated data-generating process, strongly contrasting the data-dependent performance of the causal effect estimators based on single pre-selected learners. The proposed stacking approaches thus appear relevant in the ubiquitous scenario where there is uncertainty about the set of control variables, correct functional form or the appropriate regularization assumption.

By revisiting the simulation design of wuthrich2021, we further show that stacking can outperform linear regression for even small sample sizes. We argue that the poor small sample performance of lasso-based approaches is partially driven by the choice of covariate transformations and illustrate how stacking can accommodate a richer set of specifications, including competing parametric models. We also find that short-stacking and pooled stacking may outperform DDML paired with conventional stacking in small to moderate sample sizes. Paired with its lower computational cost, this finding suggests that short-stacking may be an attractive baseline approach to select and combine competing reduced form specifications.

Finally, we demonstrate the value of pairing of DDML with stacking with two applications. First, we examine gender gaps in citations of articles published in top-30 economic journals from 1983 to 2020, and assess how the difference in citations change when conditioning on content and quality proxied by the abstract text. Estimating these conditional differences is a challenging statistical problem due to the non-standard nature of text data, which is increasingly encountered in economic applications ash2023a,chen2023gender,eberhardt2022. Second, we revisit a UK sample of the OECD Skill Survey to estimate semiparametric Kitagawa-Oxaca-Binder estimates of the unexplained gender wage gap. Both applications highlight that estimators of structural parameters based on single learners can be highly sensitive to the underlying structure of the data and/or poor tuning. The applications further demonstrate that DDML with stacking is a simple and practical solution to resolve the difficult problem of choosing a particular candidate learner in practice. Further, we observe that the optimal stacking weights often vary across reduced-form equations --- meaning that different conditional expectation functions in the same data set are best estimated using different learners. This behavior sharply contrasts with common estimation approaches, such as OLS and PDS lasso, that impose the same form for each conditional expectation function.

The remainder of the paper is organized as follows: Section (ref) provides a brief review of DDML. Section (ref) discusses DDML with stacking, short-stacking, and pooled stacking. Section (ref) presents our calibrated simulation studies. Section (ref) discusses the applications, and Section (ref) concludes.

Double/Debiased Machine Learning

This section outlines double/debiased machine learning as discussed in Chernozhukov2018. Throughout, we focus on the partially linear model as a natural extension of commonly applied linear regression methods. Despite its simplicity, the partially linear model illustrates practical challenges in the application of DDML that can be addressed by stacking. We highlight, however, that our discussion also applies to the wide range of models outlined in Chernozhukov2018 and more generally to estimation of low-dimensional structural parameters in the presence of high-dimensional nuisance functions.\footnote{A key example not explicitly discussed in Chernozhukov2018 is doubly-robust estimation of difference-in-difference parameters with staggered treatment assignment as in callaway2021difference and chang2020double. In settings with conditional parallel trends assumptions, high-dimensional nuisance functions arise in the estimation of group-time specific average treatment effect on the treated. The pairing of DDML and stacking, as proposed in this paper, also directly applies to the estimator of callaway2021difference under a conditional unconfoundedness assumption.} Stacking could also be applied in settings where DDML is applied with nonparametric targets as in ColangeloLee.

The partially linear model is defined by a random vector $(Y, D, X^\top, U)$ with joint distribution characterized by

align[align omitted — 76 chars of source]

where $Y$ is the outcome, $D$ is the scalar variable of interest, and $X$ is a vector of control variables. The parameter of interest $\theta_0$ and the unknown nuisance function $g_0$ are such that the corresponding residual $U$ satisfies the conditional orthogonality property $E[Cov(U, D\vert X)]= 0$. These properties are analogous to the orthogonality properties of residuals in multiple linear regression with the key difference here being that $g_0$ need not be linear in the controls.

Albeit a seemingly small change in specification, the partially linear model has several important advantages over linear regression. For discrete $D$, for example, results in angrist1999empirical imply that $\theta_0$ can be interpreted as a positively weighted average of incremental changes in the conditional expectation function $E[Y\vert D = d, X]$. Under appropriate conditional unconfoundedness assumptions, $\theta_0$ thus corresponds to a convex combination of conditional average treatment effects.\footnote{Similarly, for continuous $D$, $\theta_0$ corresponds to a positively weighted average of derivatives of the conditional expectation function $E[Y\vert D = d, X]$ with respect to $d$. Under a conditional unconfoundedness assumption, $\theta_0$ is thus a convex combination of derivatives of the causal response function.} Importantly, these interpretations remain valid even if the additive separability assumption of the partially linear model fails. Linear regression coefficients, in contrast, do not correspond to positively weighted averages of causal effects without imposing strong linearity assumptions that are questionable in real applications.\footnote{In the context of IV estimation where instrument validity relies on observed confounders, blandhol2022tsls emphasize that, in the absence of strong functional form assumptions, two stage least squares does not generally correspond to a convex combination of local average treatment effects (LATE). We note that the IV analogue to the partially linear model does admit a causal interpretation under the LATE assumptions, just as the partially linear model admits a weakly causal interpretation under conditional unconfoundedness.}

The advantages of the partially linear model in the interpretation of its parameter of interest come at the cost of a more challenging estimation problem relative to estimating a model that is linear in a pre-specified set of variables. Estimators for $\theta_0$ are based on the solution to the moment equation

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

given by

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

where $\ell_0(X) \equiv E[Y \vert X]$ and $m_0(X) \equiv E[D \vert X]$ are the conditional expectations of the outcome and variable of interest given the controls, respectively. Since conditional expectation functions are high-dimensional in the absence of strong functional form assumptions, a sample analogue estimator for $\theta_0$ requires nonparametric first-step estimators for the nuisance parameters $\ell_0$ and $m_0$. While nonparametric estimation generally reduces bias compared to linear regression alternatives, the increased variance associated with more flexible functional form estimation introduces additional statistical challenges: To allow for statistical inference on $\theta_0$, the nonparametric estimators need to converge sufficiently quickly to the true conditional expectation functions as the sample size increases.

DDML defines a class of estimators that allows for statistical inference on the parameter of interest $\theta_0$ while only imposing relatively mild convergence requirements on the nonparametric estimators. These mild requirements are central to the wide applicability of DDML as they permit the use of a large variety of machine learners.\footnote{The exact convergence rate requirement for nonparametric estimators depends on the parameter of interest. Chernozhukov2018 name the crude rate requirement of $o(n^{-1/4})$, but provide examples where the rate requirement is considerably weaker. Recent contributions show that these requirements are satisfied by specific instances of machine learners; see, e.g., results for lasso bickel2009simultaneous, Belloni2012, random forests wager2015adaptive, wager2018, Athey2019a, neural networks schmidt2020nonparametric, farrell2021deep, and boosting luo2016high. The exact asymptotic properties of many other machine learners remain an active research area.}

Two key devices permit the mild convergence requirements of DDML: Identification of the parameter of interest based on Neyman-orthogonal moment conditions and estimation using cross-fitting. Neyman-orthogonal moment conditions are insensitive to local perturbations around the true nuisance parameter.\footnote{In the context of the partially linear model, the formal Neyman-orthogonality requirement is

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

for arbitrary measurable functions $\ell$ and $m$, which can easily be verified using properties of the residuals.} Cross-fitting is a sample-splitting approach that addresses the {\it own-observation bias} that arises when the nuisance parameter estimation and the estimation of $\theta_0$ are applied to the same observation. In practice, cross-fitting is implemented by randomly splitting a sample $\{(Y_i, D_i, X^\top_i)\}_{i\in I}$ indexed by $I = \{1, \ldots, n\}$ into $K$ evenly-sized folds, denoted as $I_1,\ldots, I_K$. For each fold $k$, the conditional expectations $\ell_0$ and $m_0$ are estimated using only observations not in the $k$th fold --- i.e., in $I^c_k\equiv I \setminus I_k$ --- resulting in $\hat{\ell}_{I^c_{k}}$ and $\hat{m}_{I^c_{k}}$, respectively, where the subscript ${I^c_{k}}$ indicates the subsample used for estimation. The out-of-sample predictions for an observation $i$ in the $k$th fold are then computed via $\hat{\ell}_{I^c_{k}}(X_i)$ and $\hat{m}_{I^c_{k}}(X_i)$. Repeating this procedure for all $K$ folds then allows for computation of the DDML estimator for $\theta_0$:

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

where $k_i$ denotes the fold of the $i$th observation.

Since the cross-fitting algorithm depends on the randomized fold split, and since some machine learners rely on randomization too, DDML estimates vary with the underlying random-number generator and seed. To reduce dependence on randomization, it is thus worthwhile to repeat the cross-fitting procedure and apply mean or median aggregation over DDML estimates Ahrens2023_ddml. We show in Section (ref) that repeating the cross-fitting procedure is a useful diagnostic tool, allowing to gauge the stability of DDML estimators.

Under the conditions of Chernozhukov2018 --- including, in particular, the convergence requirements on the nonparametric estimators --- $\hat{\theta}_n$ is root-$n$ asymptotically normal around $\theta_0$. As already highlighted by the example in Figure (ref), however, a poorly chosen or poorly tuned machine learner for the estimation of nuisance parameters $\hat{\ell}$ and $\hat{m}$ can have detrimental effects on the properties of $\hat{\theta}_n$. Since no machine learner can be best across all settings, this raises the difficult question of which learner to apply in a particular setting. In the next section, we discuss how DDML can be paired with stacking to provide a practical solution to the choice of learner. We also illustrate how the cross-fitting structure naturally arising in DDML estimators can be leveraged to substantially reduce the computational burden otherwise associated with stacking.

Pairing DDML with Stacking Approaches

This section discusses the estimation of structural parameters by pairing DDML with stacking approaches. After the discussion of DDML with conventional stacking, we introduce two stacking variants that leverage the cross-fitting structure of DDML estimators: short-stacking and pooled stacking. To fix ideas, we focus on the nuisance parameter $\ell_0(X)=E[Y\vert X]$ arising in the partially linear model where we consider an i.i.d.\ sample $\{(Y_i,X_i)\}_{i\in I}$. Further, we consider a rich set of $J$ pre-selected base or candidate learners. The set of learners could include distinct parametric and nonparametric estimators --- e.g., linear or logistic regression, regularized regression such as the lasso, or tree-based methods such as random forests --- as well as the same algorithm with varying (hyper-) tuning parameters or different (basis) expansions of the control variables. It is important to note that the set of candidate learners for stacking can readily incorporate commonly used unregularized learners such as linear or logistic regression; in practice, sometimes the best-performing candidate learner may be one such learner.

figure[figure omitted — 1,098 chars of source]

\paragraph{DDML with conventional stacking.} Combining DDML with conventional stacking involves two layers of re-sampling, as we illustrate in Figure (ref). The {\it cross-fitting layer} divides the sample into $K$ cross-fitting folds, denoted by $I_1,\ldots,I_K$. In each cross-fitting step $k\in\{1,\ldots,K\}$, the stacking learner is trained on the training sample which excludes fold $I_k$ and which we label $T_k\equiv I\setminus I_k$. Fitting the stacking learner, in turn, requires sub-dividing the training sample $T_k$ further into $V$ cross-validation folds. This second sample split constitutes the {\it cross-validation layer}. We denote the cross-validation folds in cross-fitting step $k$ by $T_{k,1},\ldots,T_{k,V}$. Each candidate learner $j\in\{1,\ldots,J\}$ is cross-validated on these folds, yielding cross-validated predicted values for each learner.

The final learner fits the outcome $Y_i$ against the cross-validated predicted values of each candidate learner. The most common choice is to construct a convex combination via constrained least squares (CLS), with weights restricted to be non-negative and summing to one. Specifically, for each $k$, candidate learners are combined to solve \[ \underset{w_{k,1},\ldots,w_{k,J}}{\min} \sum_{i\in T_k} \left( Y_i - \sum_{j=1}^J w_{k,j} \hat{\ell}^{(j)}_{T_{k,v(i)}^c}(X_i) \right)^2 \qquad \textrm{s.t.}\ w_{k,j}\geq 0,\ \sum_{j=1}^J |w_{k,j}|=1.\] Here, $\hat{\ell}^{(j)}_{T_{k,v(i)}^c}(X_i)$ denotes the out-of-sample predicted value for observation $i$, which is calculated from training candidate learner $j$ on sub-sample $T_{k,v(i)}^c \equiv T_k \setminus T_{k,v(i)}$, i.e., all step-$k$ cross-validation folds but fold $(k,v(i))$ which is the fold of the $i$th observation. We call the resulting $\hat{w}_{k,j}$ the {\it stacking weights}. The stacking predictions are obtained as $\sum_j \hat{w}_{k,j} \hat{\ell}_{T_k}^{(j)}(X_i)$ where each learner $j$ is re-fit on $T_k$.

Although various options for combining candidate learners are available, CLS facilitates the interpretation of stacking as a weighted average of candidate learners Hastie2009. Due to this constraint, CLS tends to set some stacking weights to exactly zero. The constraint also regularizes the final estimator, which is important to mitigate issues arising from potential multicollinearity of the candidate learners. An alternative to CLS, which we refer to as {\it single-best learner}, is to impose the constraint that $w_{k,j}\in \{0,1\}$ and $\sum_{j} w_{k,j}=1$, implying that only the candidate learner with lowest cross-validated loss is used as the final estimator. Under appropriate restrictions on the data-generating process and loss function, laan2003 show asymptotic equivalence between stacking and the best-performing candidate learner.\footnote{The {\it scikit-learn} sklearn_api routines {\tt StackingRegressor} and {\tt StackingClassifier} implement stacking for Python. In Stata, stacking regression and classification are available via pystacked, which is a Stata front-end for these Python routines Ahrens2022.}

A drawback of DDML with stacking is its computational complexity. Considering the estimation of a single candidate learner as the unit of complexity (and ignoring the cost of fitting the final learner), DDML with stacking heuristically has a computational cost proportional to $K\times V \times J$. For example, when considering DDML with $K=5$ cross-fitting folds and $J=10$ candidate learners that are combined based on $V=5$ fold cross-validation, more than 250 candidate learners need to be individually estimated. Although DDML with stacking is “embarrassingly parallel” and can thus be expected to decrease in computational time nearly linearly in the number of available computing processes, the increased complexity limits its application to moderately complex applications. Another potential concern (which we investigate in Section (ref)) is that DDML with stacking might not perform well in small samples, given that candidate learners are effectively trained on approximately $\frac{(K-1)(V-1)}{KV}$% of the full sample (see Figure (ref)). These two concerns motivate short-stacking.

\paragraph{DDML with short-stacking.} In the context of DDML, we propose to take a short-cut: Instead of fitting the final learner on the cross-{\it validated} fitted values in each step $k$ of the cross-fitting process, we can directly train the final learner on the cross-{\it fitted} values using the full sample; see Figure (ref). Formally, candidate learners are combined to solve \[ \underset{w_1,\ldots,w_J}{\min} \sum_{i=1}^n \left( Y_i - \sum_{j=1}^J w_j\hat{\ell}^{(j)}_{I_{k(i)}^c}(X_i) \right)^2 \qquad \textrm{s.t.}\ w_j\geq 0,\ \sum_j |w_j|=1\] where $w_j$ are the short-stacking weights. Cross-fitting thus serves a double purpose: First, it avoids the own-observation bias by avoiding overlap between the samples used for estimating high-dimensional nuisance functions and the samples used for estimating structural parameters. Second, it yields out-of-sample predicted values which we leverage for constructing the final stacking learner. As a consequence, the computational cost of DDML with short stacking is heuristically only proportional to $K\times J$ in units of estimated candidate learners. In the example from the previous paragraph, short-stacking thus requires estimating about 200 fewer candidate learners.

We recommend DDML with short-stacking in settings where the number of candidate learners is small relative to the sample size, i.e., $J\ll n$. We believe this setting provides a good approximation to current applications of machine learning in economics and other social sciences where it is rare to consider more than a few candidate learners. If instead the number of considered learners is very large relative to the sample size --- i.e., settings in which inference for standard linear regression on $J$ variables is invalid --- pairing DDML with short-stacking may introduce bias.\footnote{Suppose, for simplicity, we consider ordinary (unconstrained) least squares as the final learner. Heuristically, the regression of $Y_i$ against $J$ sets of cross-fitted predicted values is akin to a conventional least squares regression of $Y_i$ against $J$ observed regressors where good performance would require $J/n\rightarrow 0$, ignoring that the cross-fitted predicted values are estimated. The additional regularization by constrained least squares should further weaken this rate requirement.}

figure[figure omitted — 848 chars of source]

\paragraph{DDML with pooled stacking.} While DDML with conventional stacking has one vector of weights per cross-fitting fold, short-stacking yields a single weight for each learner. A single weight for each learner decreases the variance of the final estimator and facilitates the interpretation of the stacking weights. Another way of achieving common stacking weights is DDML with pooled stacking. Pooled stacking relies on the same two-layer re-sampling strategy as conventional stacking, but combines candidate learners to solve \[ \underset{w_{1},\ldots,w_{J}}{\min} \sum_{i\in I}\sum_{k\neq k(i)} \left( Y_i - \sum_{j=1}^J w_{j} \hat{\ell}^{(j)}_{T_{k,v(i)}^c}(X_i) \right)^2 \qquad \textrm{s.t.}\ w_{j}\geq 0,\ \sum_{j=1}^J |w_{j}|=1.\] That is, pooled stacking collects the cross-validated predicted values that are calculated in each step $k$ of the cross-fitting process for each learner $j$ and estimates the stacking weights based on the pooled data set. We note that the computational costs are approximately the same as for DDML with conventional stacking.

The Practical Benefits of DDML with Stacking: Two Simulation Studies

In this section, we discuss two simulation studies illustrating the advantages of pairing DDML with stacking over alternative approaches based on single pre-selected learners. We begin with a simulation calibrated to household data on wealth and 401k eligibility from the 1991 wave of the Survey of Income and Program Participation (SIPP) in Section (ref). In Section (ref), we revisit the simulation of wuthrich2021 to assess the robustness of DDML with stacking approaches in very small samples.

Simulation calibrated to the SIPP 1991 household data

To assess the performance of DDML with conventional stacking, short-stacking and pooled stacking in a realistic setting, we consider the analysis of 401(k) eligibility and total financial assets in poterba1995 as the basis for an empirically calibrated Monte Carlo simulation. The application has recently been revisited by Belloni2017, Chernozhukov2018, and wuthrich2021 to approximate high-dimensional confounding factors using machine learning. We focus on estimating the partially linear model discussed in the previous section. The outcome is measured as net financial assets, the treatment variable is an indicator for eligibility to the 401(k) pension scheme, and the set of controls includes age, income, education in years, family size, as well as indicators for two-earner status, home ownership, and participation in two alternative pension schemes.

Algorithm\onehalfspacing \caption{Algorithm for the calibrated Monte Carlo simulation} \begin{minipage}{.93\linewidth} \begin{tcolorbox}[standard jigsaw,opacityback=0] Let $\{(y_i, d_i, x_i)\}_{i = 1, \ldots, n}$ denote the observed sample, where $i$ is a household in the 1991 SIPP and $y_i$, $d_i$, and $x_i$ respectively denote net financial assets, an indicator for 401(k) eligibility, and the vector of control variables. \begin{enumerate}[nosep] • Using the full sample, obtain the slope coefficient $\hat{\theta}_{OLS}\approx 5\,896$ from linear regression of $d_i$ against $d_i$, and $x_i$ in the original data. Construct the partial residuals $y_i^{(r)} = y_i - \hat{\theta}_{OLS} d_i, \, \forall i$. • Fit a supervised learning estimator (either linear regression or gradient boosting) to predict $y_i^{(r)}$ with the controls $x_i$. Denote the fitted estimator by $\Tilde{g}$. Similarly, fit a supervised learning estimator to predict $d_i$ with $x_i$ and denote the fitted estimator by $\Tilde{h}$. • Repeat to generate simulated samples of size $n_b$: \begin{enumerate}[nosep] • Sample from the empirical distribution of $x_i$ by bootstrapping $n_b$ observations from the original data. Denote the bootstrapped sample by $\mathcal{D}_b$. • Draw $\nu_i \overset{iid}{\sim} \mathcal{N}(0, \kappa_1)$ and $\varepsilon_i \overset{iid}{\sim} \mathcal{N}(0, \kappa_2)$, where $\kappa_1$ and $\kappa_2$ are simulation hyperparameters. Define \begin{align*} \Tilde{d}^{(b)}_i &= \mathbbm{1}\{\Tilde{h}(x_i) + \nu_i \geq 0.5\} \\ \Tilde{y}^{(b)}_i &= \theta_0 \Tilde{d}^{(b)}_i + \Tilde{g}(x_i) + \varepsilon_i \qquad\forall i \in \mathcal{D}_b \end{align*} where we set $\theta_0=6\,000$ to roughly resemble the magnitude of the regression coefficient of 401(k) eligibility in the full data. \end{enumerate} \end{enumerate} \end{tcolorbox}\end{minipage}\singlespacing\scriptsize \parbox{\linewidth}{{\it Notes:} We set the hyper-parameter $\kappa_1$ and $\kappa_2$ to approximately match variance of 401(k) eligibility and wealth in the data. The values of the simulation hyperparameters $(\kappa_1, \kappa_2)$ differ slightly depending on the supervised learning estimator used to fit the reduced form equations in the data. We take $\kappa_1= 0.35$ in both scenarios but take $\kappa_2 = 55\,500$ when using linear regression and $\kappa_2 = 54\,000$ when using gradient boosting. Differences arise because gradient boosting reduces residual variance in the true data.}

The simulation involves three steps. In the calibration step, we fit two generative models to the $n=9\,915$ households from the 1991 wave of the Survey of Income and Program Participation. The first generative model is fully linear while the second is partially linear, allowing controls to enter non-linearly through gradient-boosted trees fitted to the real data. This approach is aimed at extracting and magnifying the linear or non-linear structures in the empirical conditional distributions, respectively, enabling us to compare the performance of estimators across favorable and unfavorable structures of the data. The generative step then simulates datasets of size $n_b=\{9\,915,\,99\,150\}$ from the respective fully linear model and the partially linear model. Throughout, we set the effect of 401(k) eligibility on total financial wealth to $\theta_0=6\,000$. Finally, in the estimation step, we fit various estimators to bootstrapped samples of the generated datasets and assess their statistical properties. We outline the steps used for constructing the two generative models in more detail in Algorithm (ref).

For each bootstrap sample, we calculate estimates of the effect of 401(k) eligibility on simulated net financial assets. The estimators we consider are linear regression, the post-double selection (PDS) lasso estimator proposed by Belloni2014a, as well as DDML estimators with and without stacking. The candidate learners of the DDML estimators are linear regression, cross-validated lasso and ridge regression with interactions and second-order polynomial expansions of the controls, cross-validated lasso and ridge with no interactions but 10\textsuperscript{th}-order polynomial expansions of the controls, two versions of random forests, two versions of gradient-boosted trees, and feed-forward neural nets with three hidden layers of size five (see Table (ref) notes for details). We estimate DDML paired with conventional stacking, short-stacking and pooled stacking, and consider different methods to construct the final conditional expectation function estimator: CLS, unconstrained linear regression (OLS), selecting the single best estimator, and an unweighted average.

Table (ref) presents the mean bias, median absolute bias (MAB) and coverage rates of a 95% confidence interval associated with estimates of the effect of 401(k) eligibility on net financial assets. The left and right panels correspond to results based on data simulated from the linear (Panel A) and non-linear (Panel B) generative models, respectively. The CLS weights associated with each candidate learner are shown in Table (ref).\footnote{Further results are provided in the Appendix. Table (ref) in the Appendix gives the mean-squared prediction errors (MSPE) for each candidate learner for comparison. Table (ref) reports standard errors of the bias. Tables (ref) and (ref) show the stacking weights when using single-best and OLS as the final learner, respectively.}

sidewaystable\scriptsize\singlespacing \caption{Bias and Coverage Rates in the Linear and Non-Linear DGP} \sisetup{ table-number-alignment=left, tight-spacing=true, table-format=-4.1 } \begin{tabular}{rlS[table-format=-4.1]S[table-format=4.1]ccS[table-format=-4.1]S[table-format=4.1]ccS[table-format=-4.1]S[table-format=4.1]ccS[table-format=-4.1]S[table-format=4.1]c} \toprule \midrule && \multicolumn{7}{c}{Panel (A): Linear DGP} & &\multicolumn{7}{c}{Panel (B): Non-linear DGP}\\ \cmidrule{3-9}\cmidrule{11-17} & & \multicolumn{3}{c}{$n_b = 9\,915$} & & \multicolumn{3}{c}{$n_b = 99\,150$} & & \multicolumn{3}{c}{$n_b = 9\,915$} & & \multicolumn{3}{c}{$n_b = 99\,150$}\\ \cmidrule{3-5}\cmidrule{7-9}\cmidrule{11-13}\cmidrule{15-17} && \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{MAB} & Rate & & \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{MAB} & Rate & & \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{MAB} & Rate & & \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{MAB} &\multicolumn{1}{c}{Rate} \\ \midrule \multicolumn{4}{l}{Full sample:}\\ \IfFileExists{Simul/sim_Advantages/sim_output_folds2}{ \openin\pin@file Simul/sim_Advantages/sim_output_folds2 \setcounter{pinlineno}{1} \@whilenum\value{pinlineno}<1 \do{ \read\pin@file to\pin@line \stepcounter{pinlineno} } \addtocounter{pinlineno}{-1} \let\empty \begingroup \endlinechar\newlinechar \@whilenum\value{pinlineno}<2 \do{ \readline\pin@file to\pin@line \edef{\pin@line} \stepcounter{pinlineno} } \closein\pin@file \expandafter\endgroup \scantokens\expandafter }{ \errmessage{File `Simul/sim_Advantages/sim_output_folds2' doesn't exist!} } \multicolumn{4}{l}{DDML methods:}\\ \multicolumn{4}{l}{\it Candidate learners}\\ \IfFileExists{Simul/sim_Advantages/sim_output_folds2}{ \openin\pin@file Simul/sim_Advantages/sim_output_folds2 \setcounter{pinlineno}{1} \@whilenum\value{pinlineno}<3 \do{ \read\pin@file to\pin@line \stepcounter{pinlineno} } \addtocounter{pinlineno}{-1} \let\empty \begingroup \endlinechar\newlinechar \@whilenum\value{pinlineno}<12 \do{ \readline\pin@file to\pin@line \edef{\pin@line} \stepcounter{pinlineno} } \closein\pin@file \expandafter\endgroup \scantokens\expandafter }{ \errmessage{File `Simul/sim_Advantages/sim_output_folds2' doesn't exist!} } \multicolumn{4}{l}{\it Stacking approaches}\\ \IfFileExists{Simul/sim_Advantages/sim_output_folds2}{ \openin\pin@file Simul/sim_Advantages/sim_output_folds2 \setcounter{pinlineno}{1} \@whilenum\value{pinlineno}<13 \do{ \read\pin@file to\pin@line \stepcounter{pinlineno} } \addtocounter{pinlineno}{-1} \let\empty \begingroup \endlinechar\newlinechar \@whilenum\value{pinlineno}<24 \do{ \readline\pin@file to\pin@line \edef{\pin@line} \stepcounter{pinlineno} } \closein\pin@file \expandafter\endgroup \scantokens\expandafter }{ \errmessage{File `Simul/sim_Advantages/sim_output_folds2' doesn't exist!} } \midrule \bottomrule \end{tabular} \parbox{\linewidth}{ Notes: The table reports mean bias, median absolute bias (MAB) and coverage rate of a 95% confidence interval for the listed estimators. We consider DDML with $K=2$ cross-fit folds and the following individual learners: OLS with elementary covariates, CV lasso and CV ridge with second-order polynomials and interactions, CV lasso and CV ridge with 10th-order polynomials but no interactions, random forest with low regularization (8 predictors considered at each leaf split, no limit on the number of observations per node, bootstrap sample size of 70%), highly regularized random forest (5 predictors considered at each leaf split, at least 10 observation per node, bootstrap sample size of 70%), gradient-boosted trees with low regularization (500 trees, maximum depth of 3 and a learning rate of 0.01), gradient-boosted trees with high regularization (250 trees, maximum depth of 3 and a learning rate of 0.01), feed-forward neural nets with three hidden layers of size five. For reference, we report two estimators using the full sample: OLS and PDS lasso. Finally, we report results for DDML paired with conventional stacking, short-stacking and pooled stacking where the final estimator is either CLS, OLS, the unweighted average of candidate learners or the single-best candidate learner. Results are based on 1\,000 replications. }

Given the construction of the generative models, we would expect that linear regression performs best in the fully linear setting and that DDML with gradient boosting performs best in the nonlinear setting where the nuisance function is generated by gradient boosting. The simulation results confirm this intuition, showing that the two procedures achieve among the lowest bias and median absolute bias in the data-generating processes that are based on them. Researchers are rarely certain of the functional structure in economic applications, however, so that it is more interesting to consider their respective performance in the non-favorable setting. In the non-linear data-generating process, linear regression is among the estimators with the worst performance across all three measures. Similarly, gradient boosting-based DDML is non-optimal in the linear data-generating process. It is outperformed by linear regression and CV lasso, both of which enforce a linear functional form on the control variables, in terms of MAB.

The simulation results are consequences of the “no free lunch” theorem in machine learning wolpert1996. Informally, the theorem states that there exists no estimator that performs best across all empirical settings. Researchers must, therefore, carefully match estimators to their application. However, with limited knowledge about underlying data-generating processes and few functional form restrictions implied by economic theory, the number of plausibly suitable estimators is typically large.

The bottom section of Table (ref) reports results for DDML combined with the three stacking approaches outlined in Section (ref). For each stacking approach, we consider stacking weights estimated by (CLS) as outlined in Section (ref), set equal to $1/J$ (Average), estimated without constraint by OLS (OLS), and by selecting only the single best candidate learner (Single-best). We find that short-stacking performs similarly to, and sometimes better than, conventional and pooled stacking, while being computationally much cheaper (as shown in Table (ref)). For example, at $K=10$ and $V=5$, DDML combined with short-stacking ran around 4.3 times faster on the full sample than DDML with conventional or pooled stacking, which is roughly in line with a speed improvement by a factor of $1/V$.\footnote{The computations were performed on the high-performance cluster of the ETH Zurich. Each instance used a single core of an AMD EPYC processor with 2.25-2.6GHz (nominal)/3.3-3.5 GHz (peak) and 4GB RAM. The run time of DDML with conventional stacking was $2\,393$s on the full sample, while short-stacking ran in only 540s.}

table[table omitted — 5,159 chars of source]

The simulation set-up is favorable to using single-best as the final learner because there is one `true' candidate learner. However, single-best does not visibly outperform CLS, although single-best always selects the correct learner in the non-linear DGP if the sample size is sufficiently large (see Appendix Table (ref)). We believe that CLS is, in practice, a good default choice. In a setting where there is a single learner that does a distinctly better job approximating the target conditional expectation function, CLS should assign a very large weight to that learner and thus approximate single-best. Otherwise, when there are several learners with similar performance or learners that perform differentially well for different observations, there are potential gains from combining the different learners.

The bias of the OLS final learner is overall similar to CLS, except when employing conventional stacking under the non-linear DGP for $n_b=9\,915$ where the average bias is almost four times as large.\footnote{Appendix Table (ref) shows that the OLS stacking weights are often outside the unit interval. The weights associated with the neural net are particularly large (in absolute value), suggesting that OLS might be more sensitive to outliers than CLS.} The unweighted average appears sub-optimal for $n_b=99\,150$ under the non-linear DGP. Poor performance of unweighted averaging is to be expected in settings where, as in our setup, the candidate set includes learners that are not well-matched to the DGP.\footnote{By contrast, the unweighted average is known to often perform well in time-series settings clemen1989,timmermann2006, and in particular when the optimal weights are close to equal (see Wang2023review, Section 2.6 for a summary and wider discussion). Such a scenario is ruled out in our simulation setup, which captures the situation where the researcher does not know whether a linear or non-linear learner would be more appropriate.} We also note that the computational advantage of short-stacking with the unweighted average over short-stacking with estimated weights amounts to one (constrained) regression per conditional expectation function and is thus minimal.

The CLS weights in Table (ref) indicate that stacking approaches successfully assign the highest weights to the estimators aligning with the data-generating process (i.e., either OLS or gradient boosting) among the ten included candidate learners, illustrating the ability to adapt to different data structures. Specifically, the stacking methods applied to the linear data-generating process assign the largest weight to linear models while they assign the largest weights to the gradient-boosting estimators and the lowest weights to estimators that impose a linear functional form on the control variables in the non-linear data-generating process.\footnote{The rates at which each candidate learner is selected by the single-best final learner are shown in Table (ref) in the appendix and provide similar insights.} We conclude that DDML paired with stacking approaches reduces the burden of choice researchers face when selecting between candidate learners and specifications by allowing for the simultaneous consideration of multiple options, thus implying attractive robustness properties across a variety of data-generating processes.

DDML and Stacking in Very Small Samples

A possible concern for estimators relying on machine learning is that they might not perform well for very small samples, given that their flexibility comes at the cost of increased variance compared to parametric estimators. wuthrich2021 use two simulations to demonstrate that PDS lasso tends to underselect controls, which may result in a substantial small-sample bias. They also show that the bias heavily depends on the exact lasso penalty chosen (i.e., whether the plugin penalty of Belloni2014a, is scaled by 0.5 or 1.5), and argue in favor of OLS with appropriately chosen standard errors over PDS lasso in high-dimensional settings.

We revisit the 401(k) simulation set-up in WZ to assess if DDML with stacking suffers from similar issues in small samples and to compare the performance of DDML paired with stacking with PDS lasso and OLS. Following WZ, we run simulations on bootstrap samples of the data for $n_b = \{200, 400,$ $800, 1\,600\}$ and approximate the bias as the mean difference relative to the full-sample estimates ($n=9\,915$).\footnote{The full-sample estimates are reported in Table (ref).} WZ consider two sets of controls: two-way interactions (TWI), and quadratic splines with interactions (QSI) Belloni2017. The number of predictors is 167 and 272, respectively. Figure (ref) replicates the main results of WZ (Figure 8 in their paper). Panels (a) and (b) show the bias relative to the full sample estimate for the TWI and QSI specification based on OLS and PDS lasso with tuning parameter equal to the plugin penalty of Belloni2014a scaled by $c$ for $c \in \{0.5,1,1.5\}$. It is noteworthy that the speed at which the bootstrapped estimates converge to the full-sample estimate depends on the set of controls for the PDS lasso, but less so for OLS. While PDS lasso with $c=\{0.5,1\}$ and OLS perform similarly if QSI controls are used, PDS lasso converges much more slowly to the full-sample estimate with TWI controls.

figure[figure omitted — 924 chars of source]

The DDML-stacking framework allows us to choose between, and combine, OLS and lasso with both the TWI and QSI set of controls. Another advantage of DDML over PDS lasso is that we can leverage lasso with cross-validated penalization for a fully data-driven penalization approach. Table (ref) compares the performance of the full-sample estimators OLS and PDS lasso (shown in Panel A) to DDML-stacking estimators only relying on OLS and CV lasso with TWI and QSI controls as candidate learners (Panel B). We again consider conventional stacking, short-stacking and pooled stacking together with either CLS or single-best as the final learner. We set the number of cross-fitting folds to $K=10$ (but also consider $K=2$ below for comparison in Panel E).

table[table omitted — 5,488 chars of source]

Across all sample sizes, the DDML-stacking estimators strictly outperform both OLS specifications, as well as PDS lasso with TWI, and exhibit overall similar performance to PDS lasso utilizing QSI controls and $c=\{0.5,1\}$. The differences across DDML-stacking estimators are relatively minor. The CLS short-stacking weights reported in Table (ref), Panel A-B, reveal that CV-lasso with QSI controls receives the largest weights, while both OLS specifications contribute jointly between nearly zero (at $n_b=200$) and only up to 15% (for $n_b=1\,600$ and the estimation of $E[D|X]$). When selecting only a single candidate learner, CV-lasso with QSI is chosen in more than three-fourths of bootstrap iterations for the estimation of $E[Y\vert X]$ and $E[D\vert X]$ (Panel C-D in Table (ref)), suggesting that CV-lasso with QSI controls is strictly preferable over OLS and lasso with TWI controls in this application. This simulation exercise again highlights that relying on poorly chosen specifications that are not validated against other choices might be sub-optimal. In practice, the researcher does not know whether TWI or QSI controls perform better and whether to use OLS or lasso. Crucially, DDML paired with stacking allows for simultaneous consideration of OLS and lasso with both TWI and QSI controls and thus resolves the choice between learners and control specifications in a data-driven manner.

table[table omitted — 2,128 chars of source]

In the next step, we expand the set of candidate learners by two types of random forests, two types of gradient-boosted trees and a feed-forward neural net. In principle, widening the set of candidate learners increases robustness to a larger class of unknown confounding structures. We show the results in Panel C, Table (ref). When measuring performance based on the difference to the full-sample estimates, we find there are benefits of extending the set of candidate learners for bootstrap sample sizes of $n_b=800$ or larger. The results are generally comparable across conventional, short and pooled stacking. However, single-best exhibits a lower bias for small bootstrap sample sizes vs.\ CLS, while pooled stacking with CLS appears to perform worse. The CLS weights reported in Appendix Table (ref) illustrate how DDML-stacking estimators adapt to the sample size. For example, for smaller sample sizes, a larger weight is put on OLS in the estimation of $E[Y\vert X]$. In Panel D of Table (ref), we report results for each candidate learner individually. DDML-stacking approaches perform better than most individual candidate learners and similar to the best-performing individual learner, which is DDML with CV lasso and QSI controls. In Panel E, we also show results if we reduce the number of folds to $K=2$. The performance deteriorates drastically for smaller sample sizes, indicating that---while DDML stacking appears competitive for small sample sizes---it is important to increase the number of folds to ensure larger training samples for the CEF estimators.

A drawback of measuring the bias as the difference to the full-sample estimate is that we do not gain insights about convergence to the true parameter. We thus revisit the calibrated simulation exercise from Section (ref), which allows us to measure the bias as the difference to the true parameter. When the DGP is linear (see Figure (ref)), DDML with short-stacking or pooled stacking and using CLS performs overall similarly to OLS. DDML with conventional stacking exhibits relatively large bias with $n_b=200$. If the true DGP is non-linear, see Figure (ref), OLS and PDS-Lasso are unable to recover the true effect, while DDML with short and pooled stacking and CLS yield reasonably close approximations of the true parameter even for small sample sizes. DDML with conventional stacking is competitive only for larger samples. We provide extensive results for mean bias and coverage rates in Tables (ref)--(ref) in the Appendix.

figure[figure omitted — 1,060 chars of source]

To conclude, the results highlight the risks of relying on inappropriate functional form assumptions. DDML paired with stacking approaches---when combined with a diverse set of candidate learners---imposes weaker conditions on the underlying data-generating process compared to relying on a single pre-selected learner. Short-stacking and pooled stacking outperform conventional stacking in small samples. We conjecture the improvement is due to short and pooled stacking imposing common weights across cross-fitting folds. The use of a common set of weights imposes regularization that is consistent with learner performance being stable across subsamples, which seems like a natural benchmark. For pooled stacking, this additional regularization should reduce weight variance while coming with relatively little cost in terms of bias. Because short-stacking does not make use of the additional train and test splits from conventional stacking, it has more potential to suffer from an additional over-fitting bias. When a small number of learners is considered, this additional bias should be small relative to the variance reduction obtained from not needing to estimate fold-specific stacking weights. The simulation results provide support in favor of this conjecture.

Applications

In this section, we use two applications to illustrate how pairing DDML and stacking can increase the robustness of structural parameter estimates to the underlying structure of the data. In the first application, we estimate gaps in citations of articles in top economics journals across different gender compositions among the authors. We condition on the abstract to proxy for the content and quality of the paper and demonstrate that stacking-based DDML is a practical solution to challenging estimation problems using text data. In the second application, we revisit the UK sample of the OECD Skills Survey for Kitagawa-Oaxaca-Binder estimates of the unexplained gender wage gap where we condition on a large set of individual characteristics. Both applications pertain to the literature on gender gaps in various domains, e.g., entry to STEM programs card2021, ICT literacy siddiq2019 or wages Strittmatter2021,bonaccolto-topfer2022, and are methodologically also closely related to the broader literature on discriminatory attitudes towards minority groups hangartner2021monitoring.

Gender gap in citations

This section uses DDML with stacking to estimate a partially linear model applied to average differences in citations of articles published in top-30 economic journals from 1983 to 2020 by the gender composition of the authors. Following card2020, we distinguish between papers with (imputed) all-male, all-female, and mixed-gender authorship.\footnote{As we explain below, we impute the gender mix of authors from the authors' names.} Instead of conditioning on hand-coded characteristics such as JEL codes, we leverage the abstract text as a proxy for the topic and quality of the article. Estimating these conditional differences is a challenging statistical problem due to the non-standard nature of text data, and researchers are faced with two key decisions when operationalizing an estimator using text data: how to encode the text data into numerical features, and how to select a suitable learner given the encoded data. Both decisions are ex-ante challenging, but also practically highly relevant as text data is becoming increasingly encountered in economic applications Gentzkow2010,chen2023gender,widmer2023media. We show that these decisions can be consequential and that by simultaneously considering different encoding procedures and multiple learners, DDML with stacking provides a simple practical solution to both problems.

In documenting average differences in citations, the analysis presented also contributes to the broader literature on gender biases in academia lundberg2019,card2020,hengel2022. It is well-documented that women are under-represented in academia, especially in senior positions ceci2014,lundberg2019. A possible reason for the persistent gap in representation include is that scholarly work produced by women faces more sceptical scrutiny compared to work produced by their male counterparts hengel2022,krawczyk2016. Higher scrutiny could be, for example, reflected at the refereeing stage when a publication decision is made and, as we examine here, after publication when scholarly work is attributed by other scholars through citations card2020,roberts2020a,grossbard2021.

Throughout our analysis, we focus on a descriptive characterization of the average gaps in citations across different gender compositions of the authors as given by $\theta_0$ in the partially linear model of Equation (ref) where $Y$ denotes log-citations, $D$ is a two-dimensional vector whose first component is an indicator for all-female authorship and whose second component is an indicator for mixed-gender authorship. The vector $X$ collects the content of the abstract and a set of year-of-publication indicators. The two components of $\theta_0$ may thus be interpreted as summarizing the average relative difference in total citations between all-male and all-female authorship, and all-male and mixed-gender authorship, respectively, conditional on the article's year of publication and abstract. Throughout, we make no conditional unconfoundedness assumptions that would be necessary for causal interpretations.

We consider a sample of $27\,599$ articles that have been published between 1983--2020. The data was sourced from Scopus and is a sub-sample of the data analyzed in advani2021race, who kindly shared their data with us. For each article, we have a record of the citation count and the authors' names, which we use to infer the authors' gender.\footnote{ We use the software Namsor, which frequently ranks among the best-performing algorithms for gender classification using names sebo2021,krstovski2023inferring and is widely used in academic studies Bursztyn2021,sebo2023. In the main specification, we exclude articles of authors whose gender could not be classified with a probability of less than 70%, but we show that results are similar when we apply thresholds of 60% or 90%; see Appendix Figure (ref). Our sample includes 586 articles for which no citation is recorded. These were excluded from the analysis. We also provide results using the number of citations (instead of log-citations) in Appendix Table (ref).} In the sample, 6.3% of articles are authored by only female authors and 22.9% have authors from both genders.

Before turning to estimation, the text of the abstract needs to be transformed into a numerical vector. To admit estimation conditional on the content of the abstract, it is necessary to find a representation (referred to as an embedding) of the text that is lower-dimensional but captures its core meaning. An active literature in statistics and computer science provides solutions to this problem, suggesting a large variety of algorithms to construct text embeddings ash2023a. Thus, in addition to the choice of candidate learner, researchers intent on using text data for their analysis are faced with the additional choice of embedding algorithm. To illustrate how stacking-based DDML can help support this choice, we consider two procedures for encoding the text of the abstract into numerical features: First, we consider a bag-of-word model summarizing the text as (stemmed) word counts enke2020,Esposito2023. In our data, this results in a 211-dimensional vector of word counts for each abstract. Second, since the bag-of-word approach disregards the word order and context, we construct word embeddings generated by a pre-trained BERT model, a transformer-based large-language model devlin2018bert. In particular, for each abstract, we extract the 768-dimensional vector of weights from the last hidden layer of the BERT model that was pre-trained on a large corpus of (uncased) English text data.\footnote{The model bert-base-uncased is freely available from, among others, the Python library huggingface.} Instead of embedding individual words, BERT attempts to reconstruct both whole sentences and the context of these sentences, making it particularly suitable to characterize the content of the abstracts. Recently, bajari2023hedonic use BERT to construct embeddings of product descriptions on Amazon.com.

The numerical abstract embeddings are then used in several base learners. We consider OLS, PDS lasso and DDML with CV lasso, CV ridge, XGBoost chen2016, random forests and a feed-forward neural net (see table notes for details).\footnote{To reduce the run time, we use regression approaches both for the estimation of $E[Y\vert X]$ and $E[D\vert X]$.} The base learners are aggregated by pairing DDML with either conventional stacking or short-stacking, and with either CLS or single-best.\footnote{We omit pooled stacking from this application since the R package ddml, which was used for this application, does not currently support pooled stacking.} The final estimator thus simultaneously combines both text embedding algorithms and machine learning algorithms.

figure[figure omitted — 2,284 chars of source]

Figure (ref) shows estimates of the average relative difference in total citations between all-male and all-female authorship (top-panel) and all-male and mixed-gender authorship (bottom-panel), respectively, for different control specifications and estimators. When we only condition on the publication year, the citation penalty for all-female authorship is close to zero, while there is a large positive effect of $+22.6$% ($s.e.=1.7$) for mixed-gender authorship. We next employ PDS lasso to add the abstract text either in the form of word counts or as BERT features. Using the latter, the citation gap increases to $-7.4$% ($2.9$) for articles with all-female authorship, while the average relative difference of articles with mixed-gender authorship reduces to $+9.2$% ($1.7$). The estimates are qualitatively similar when using word counts instead of BERT features.

In the figure, we also show five cross-fitting repetitions of pairing DDML with each candidate learner (in green) and the median aggregates over these repetitions (in orange). There are considerable differences across DDML estimators, with the median estimates of the citation gap ranging between $-4.7$% ($3.8$) and $-10.6$ ($3.1$) for articles with all-female authorship and between $-1.4$ (2.2) and $+10.4$ ($1.8$) for articles with mixed-gender authorship, highlighting that different candidate learner specifications can yield vastly different effect sizes. These stark differences emphasize the need to choose and tune CEF estimators carefully. Without thoroughly validating each candidate learner, judging which results are more credible is difficult. Furthermore, it is noteworthy that some candidate learners, specifically those based on XGBoost in this example, exhibit substantial instability across cross-fitting repetitions.

We show results from pairing DDML and stacking approaches on the right-hand side of the same figure. Relative to the DDML estimates based on the individual candidate learners, the stacking approaches yield lower variability over cross-fitting repetitions. All four stacking-based approaches agree on an average relative difference in citations of around $-7.0$% ($2.9$) for articles with all-female authorship and suggest a citation advantage of between $+6.2$ ($1.7$) and $+7.0$ ($1.7$) for articles with mixed-gender authorship.

table[table omitted — 2,423 chars of source]

Table (ref) shows the stacking weights of conventional and short-stacking with constrained least squares as the final learner along with mean-squared prediction errors. The stacking estimators assign small weights to learners exhibiting a relatively large MSPE and large variability over cross-fitting repetitions. For example, in the CEF estimation of log citations, the neural nets have an MSPE that is around 50% larger than that of other learners. Stacking assigns, as desired, zero weights to the neural nets, whereas OLS leveraging BERT as one of the best-performing learners receives the largest weights. It is noteworthy that the stacking weights often vary markedly across CEFs, highlighting that there is no reason to assume that the same learner is best suited for estimating both $E[Y\vert X]$ and $E[D\vert X]$. This insight is especially important since most estimation approaches (including OLS and PDS lasso) impose the same structure for each CEF.

The results on the citation gaps in top economic journals conditional on the content of the abstract are consistent with a citation penalty for all-female authored articles, possibly due to a higher degree of skepticism towards all-female author teams compared to all-male author teams. However, similar to card2020 and maddi2021, the estimates also suggest a conditional citation advantage of articles with mixed-gender authorship.

Gender gap in wages

The gap in wages between men and women is a central measure of economic gender equality and has been the focus of an extensive empirical literature blau2017. The classic approach to estimating the unexplained gender wage gap relies on a linear version of the Kitagawa-Oaxaca-Binder decomposition (kitagawa1955,oaxaca1973,blinder1973; for an overview, see fortin2011). Several recent articles by bonaccolto-topfer2022, Strittmatter2021, boheim2021decomposition and bach2023, among others, focus instead on semi-parametric decompositions of the wage gap leveraging more flexible machine learning algorithms. Much of this literature focuses, however, on lasso-based approaches, even though there is no apparent reason to favor sparsity-based approaches over learners relying on other regularization assumptions. In contrast to the recent literature that primarily focuses on lasso-based approaches, we consider a diverse set of candidate learners and aggregate them via stacking.

The parameter of interest in this application is the unexplained gender wage gap, which is the expected difference in wages after conditioning on observed characteristics. Formally,

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

where $Y$ denotes the logarithm of wages, $D$ is an indicator equal to one for women, and $X$ is a vector of potentially many individual characteristics. The parameter is well-defined if $P(D=1\vert X ) > 0$ with probability 1.\footnote{As in the previous section, we focus our analysis on a descriptive parameter of interest and do not make conditional unconfoundedness assumptions that would be necessary for causal interpretations.}

In the absence of functional form assumptions, estimation of $\theta_0$ is a challenging statistical problem due to its dependence on unknown conditional expectation functions that need to be nonparametrically estimated. Analogous to the DDML estimator for the partially linear model outlined in Section (ref), we consider estimation of $\theta_0$ via the split-sample analogue of the efficient score function for $\theta_0$\footnote{See, e.g., equation (5.4) in Chernozhukov2018.} -- i.e.,

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

where $\hat{g}_{I^c_{k}}$ and $\hat{m}_{I^c_{k}}$ are cross-fitted estimators for $g_0(D, X)\equiv E[Y\vert D, X]$ and $m_0(X) \equiv E[D\vert X]$, and {$\hat{p}_{I^c_{k}}$ is a cross-fitted estimator of $P(D=1)$.

Following forshaw2024gender, we take the data for this application from the UK sample of the OECD Skills Survey, which was collected in 2011-12 and comes with a rich set of covariates, including age, experience, education, occupation, and industry. The final data includes $4\,836$ British respondents. We specify three sets of control variables. The simple set of controls only includes a selection of essential covariates: age (in levels and squared), years of education, a literacy and numeracy test score, years of tenure in the current job (in levels and squared), education level, hours worked per week, and number of children. The base set of controls adds, among others, management level, age of children, and parents' education level.\footnote{ The base set adds the following variables to the reduced set: area of study, part of larger organization, management position, type of contract, job satisfaction, health status, living with a partner, age of youngest child, immigration age, mother's and father's highest level of education, immigration status of parents, informal job-related education in last 12 months, informal non-job-related education in last 12 months.} Furthermore, we interact age and tenure with all categorical covariates. The extended set of controls comprises all variables and interacts each continuous covariate with each categorical covariate.

We include a diverse set of candidate learners to allow for a high level of flexibility. We employ regression approaches for the estimation of the CEF of log wages (i.e., $E[Y\vert D, X]$) and classification approaches for the CEF estimation of gender (i.e., $E[D\vert X]$). Our candidate learners are linear (or logistic) regression with the simple and base set of controls; linear (or logistic) CV-lasso and CV-ridge with the base and extended set of controls; three random forests with 500 regression (or classification) trees and minimum leaf sizes of 1, 50 and 100; two types of gradient-boosted regression (or classification) trees with and without early stopping (500 trees, maximum depth of 3); two feed-forward neural nets with hidden layer sizes of $(40,20,1,20,50)$ and $(30,30,30)$ and early stopping. Finally, we aggregate the candidate learners via conventional, short and pooled stacking, using either CLS or single-best as the final learner.

figure[figure omitted — 1,667 chars of source]

Figure (ref) reports results for individual candidate learners (on the top) and stacking approaches (on the bottom). We show results from 10 cross-fitting repetitions (in green) and the median aggregates (in orange). We again find that some candidate learners exhibit substantial variability over cross-fitting repetitions, which is also reflected in the large median-aggregate standard errors. The variability is especially large for CV-ridge with the extended set of controls and the neural nets, which are the candidate learners exhibiting the largest MSE (see Appendix Table (ref)). The stacking results are, in contrast, relatively stable over cross-fit repetitions when using CLS as the final learner. Interestingly, the results when relying on single-best as the final learner are noticeably more variable than when using CLS, indicating that a combination of candidate learners seems to better fit the data than a single learner. The stacking weights and MSE in Appendix Table (ref) confirm that there is no single candidate learner dominating the others. The instability of the single-best final learner is reflected in the stacking standard errors. Given this potential for instability of choosing a single candidate learner, we recommend favoring constrained least squares over single-best if one is not confident that one of the chosen learners will be significantly better than the rest, and thus stably selected, which seems likely to be the most common setting in practice.

Conclusion

This article assesses the performance of DDML estimators in realistic settings using applications and simulation studies calibrated to real economic data. We highlight that estimators of structural parameters based on single pre-selected (machine) learners can be highly sensitive to the underlying structure of the data and/or poor tuning, and we show that pairing DDML with stacking can help alleviate these concerns, provided that a sufficiently diverse set of candidate learners is considered.

We discuss pairing DDML with conventional stacking but also suggest two alternative stacking approaches: Short-stacking, which substantially reduces the computational burden by leveraging the cross-fitting naturally arising in the computation of DDML estimates, and pooled stacking, which decreases the variance of the stacking estimator by imposing common stacking weights over cross-fitting folds. In our simulations, both strategies are competitive with conventional stacking in settings with large and moderate sample sizes and are better in small samples. The advantages of short-stacking are particularly worth highlighting, given its substantially lower computational cost.

A key advantage of the DDML-stacking approach is that it accommodates both traditional parametric and nonparametric specifications by allowing simultaneous consideration of, for example, OLS with several sets of controls, sparsity-based learners, tree-based ensembles and neural networks. In this sense, researchers are not forcibly deviating from standard (often linear) specifications unless the data suggests there is reason to. While machine-learning-based causal methods may yield fundamentally different results from linear regression only in specific examples, the additional robustness to unexpected structures in the data thus seems to come at relatively little cost.

\printbibliography

\singlespacing

\pagestyle{empty}