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.
68,680 characters · 16 sections · 5 citation commands
{ { Keywords:} Instrumental variables, local average treatment effects, Empirical Monte Carlo Study. { JEL Classification:} C21, C26. \hrule
Correspondence to: Hugo Bodory, University of St.\ Gallen, Varnb\"{u}elstrasse 14, CH-9000 St.\ Gallen, [email removed]; Martin Huber, University of Fribourg, Bd.\ de P\'{e}rolles 90, CH-1700 Fribourg, [email removed]; Michael Lechner, University of St.\ Gallen, Varnb\"{u}elstrasse 14, CH-9000 St.\ Gallen, [email removed]. **: Michael Lechner is also affiliated with CEPR and PSI, London, CESIfo, Munich, IAB, Nuremberg, and IZA, Bonn. \thispagestyle{empty}
\setcounter{page}{0}
\setcounter{page}{1}
The evaluation of the causal effect of a treatment (e.g.,\ fertility) on an outcome (e.g.,\ labor supply) is frequently complicated by endogeneity, implying that the treatment is associated with unobserved characteristics affecting the outcome (e.g.\ personality traits, preferences, and values concerning family and working life). One may nevertheless assess treatment effects in the presence of an instrumental variable (IV) which affects the treatment of (at least) some subjects in a monotonic way, does not directly affect the outcome (other than through treatment) and is as good as randomly assigned. Under these conditions, the local average treatment effect (LATE) on the compliers, the subpopulation whose treatment state reacts positively to the instrument, is identified, as discussed in \citeasnoun{Imbens+94} and \citeasnoun{Angrist+96}. In many empirical contexts, it may seem unlikely that the IV assumptions hold unconditionally, in particular when the treatment evaluation relies on observational data in which the instrument is not explicitly randomized like in an experiment. Depending on the application, it might, however, appear plausible that the IV assumptions hold conditional on covariates observed in the data. In this case, the LATE is identified and can be consistently estimated under certain conditions, see the discussions in \citeasnoun{Abadie00}, \citeasnoun{Tan2006}, and \citeasnoun{Froel2007}.
This paper assesses the finite sample performance of various parametric, semi-parametric, and non-parametric IV estimators when controlling for a fixed (i.e., pre-defined and low-dimensional) set of covariates by Monte Carlo simulations that are based on empirical labor market data from \citeasnoun{Angrist+98}. The latter study assesses the effect of fertility, defined as having at least three vs.\ two children, on mother's labor supply (for instance, a binary employment status or weeks employed per year), using twins at the second birth as instrument. The intuition for this IV strategy is that if a mother with one child get twins at the second birth, then fertility immediately increases to three rather than two children, implying a first stage effect of the twins instrument on the treatment. In the spirit of \citeasnoun{HuLeWu13}, our empirical Monte Carlo simulation makes to a certain extent use of the empirical associations in the labour market data when assessing the various IV estimators, with the aim that our analysis is more closely linked to real world applications.
We vary the simulation designs in several dimensions, including treatment effect heterogeneity, instrument selectivity across observed covariates (namely age, race, and quarter of birth), instrument strength, the outcome distribution, and sample size. We analyse the performance of a range of estimators commonly considered in treatment and policy evaluation, including two stage least squares, inverse probability weighting, matching, doubly robust estimation, and nonparametric regression. We find that overall, non-parametric estimation based on the random forest, a machine learning algorithm controlling for covariates in a data-driven way, performs best in terms of coverage rates. The latter are defined as the share of simulations in which the true LATE is included in the 95% confidence interval of a LATE estimator, where the standard error required for the construction of the confidence interval is obtained by the non-parametric bootstrap. Furthermore, the random forest-based estimator is relatively precise, implying that the confidence interval is comparably short, which (conditional on having a decent coverage) appears desirable from the perspective of statistical power. Non-parametric kernel regression as well as certain versions of semi-parametric radius matching on the propensity score, pair matching on the covariates, and inverse probability weighting also have a decent coverage, but are less precise than the random forest-based method. Concerning the average root mean squared error of LATE estimation, kernel regression performs best (and also has the smallest average standard deviations), closely followed by the random forest method, which has the lowest average absolute bias.
The remainder of this paper is organized as follows. Section (ref) discusses the identifying assumptions for IV-based LATE evaluation in the presence of covariates. Section (ref) introduces various parametric, semi-parametric, and non-parametric LATE estimators, as well as a bootstrap procedure for computing standard errors. Section (ref) presents our empirical Monte Carlo simulation approach, namely the empirical data and the simulation designs. Section (ref) presents the results on the finite sample performance of the LATE estimators. Section (ref) concludes.
In this section, we present the assumptions underlying the identification of the Local Average Treatment Effect (LATE) when controlling for covariates. To formalize the discussion, let us denote by $D_i$ a possibly endogenous treatment received by unit $i$, and by $Y_i$ the outcome variable based on which the treatment effect is to be evaluated. In their seminal paper, \citeasnoun{Imbens+94} define the LATE as the mean effect of $Y_i$ in response to a change in $D_i$ among the compliers, a subgroup whose $D_i$ reacts to an exogenous shift in the instrumental variable, which is denoted by $Z_i$. To discuss the identification of the LATE, we make use of the potential outcomes framework introduced by \citeasnoun{Rubin74}, which expresses causal effects as differences between potential outcomes under treatment and non-treatment. We adapt this concept to our instrumental variable setting with binary indicators $D_i$ and $Z_i$, and define potential outcome and treatment variables for unit $i$ in the following way:
with $d_i,z_i\in\{0,1\}$. Using this framework, \citeasnoun{Angrist+96} show that units can be divided into two subgroups, compliers and noncompliers. Compliers are those induced to take the treatment when being assigned to it. Formally, this type of units is characterized by $D_{i,1}-D_{i,0}=1$. The subgroup of noncompliers may consist of three further types, namely always-takers with $D_{i,1}=D_{i,0}=1$, never-takers with $D_{i,1}=D_{i,0}=0$, and defiers with $D_{i,1}-D_{i,0}=-1$. Note that the type of a single unit cannot be identified because the counterfactual potential treatment (that would have occurred under the alternative, rather than the factual instrument assignment) is not observed.
\citeasnoun{Abadie00}, \citeasnoun{Tan2006}, and \citeasnoun{Froel2007} consider non-parametric LATE identification and estimation when controlling for observed covariates, denoted by $X_i$. We subsequently present the identifying assumptions in this context, which consist of (i) a monotonicity restriction on the treatment, (ii) the existence of compliers, (iii) conditional independence of the instrument and the share of compliance types, (iv) conditional mean independence of the outcome and the instrument, and (v) common support.
Assumption (ref) rules out the presence of defiers, a type whose treatment never complies with the instrument. Assumption (ref) implies that the subgroup of compliers exists. Due to the conditional independence of the instrument and the shares of compliers, always-takers, and never-takers stated in Assumption (ref), the first stage effect of the instrument on the treatment is identified conditional on covariates, such that any variables affecting both the instrument and the treatment are controlled for. The conditional mean independence in Assumption (ref) rules out a direct average effect of the instrument on the outcome (exclusion restriction) and unobservables that jointly affect the instrument and the outcome when controlling for covariates. Finally, Assumption (ref) ensures that for all covariate values occurring in the population, either instrument value $Z_i\in\{0,1\}$ exists such that the instrument is not deterministic in the covariates.
Under Assumptions (ref) to (ref), the LATE, denoted as $\theta=E[Y_{i,Z_i}^1-Y_{i,Z_i}^0\vert D_{i,1}-D_{i,0}=1]$, is identified by
Based on the insights of \citeasnoun{rosenbaum1983}, \citeasnoun{Froel2007} shows that identification is also obtained by conditioning on the instrument propensity score $p(x):=P(Z_i=1 \vert X_i=x)$ rather than the covariates, because it possesses the so-called `balancing property'. That is, conditioning on the one-dimensional propensity score balances the distribution of the covariates across the states of the instrument. For this reason, the LATE is alternatively identified by
In this section, we present parametric, semi-parametric, and non-parametric methods for estimating the LATE parameter $\theta$ introduced in Section (ref). We also discuss a trimming rule that tackles limited common support in covariate values across instrument states, based on dropping observations which would obtain large weights in the estimator because their covariate values occur (almost) exclusively in only one of the instrument states. Finally, we provide an bootstrap procedure for estimating the standard errors of the LATE estimators.
One method for the estimation of $\theta$ frequently applied in empirical work is two-stage least-squares (2SLS), which is easy to implement and computationally fast. However, the linearity assumption of the 2SLS estimator implies effect homogeneity, a restriction that may not hold in empirical studies. We consider 2SLS as a benchmark method, but also include more general LATE estimators that allow for effect heterogeneity of the LATE across values of the covariates.
Equations (ref) and (ref) imply that $\theta$ can be expressed as the ratio of two treatment effect estimators that account for covariate differences in the presence and absence of the instrument. The numerator gives the reduced form effect of $Z_i$ on $Y_i$ and the denominator the first stage effect of $Z_i$ on $D_i$. Thus, a natural choice for the construction of estimators for $\theta$ is to substitute the expressions in the numerators and denominators of Equations (ref) and (ref) by estimators standardly applied in treatment or policy evaluation, see for instance the surveys by \citeasnoun{Imbens03} and \citeasnoun{ImWo08}.
Many treatment effect estimators are semi-parametric in the sense that (parametric) propensity score estimation is combined with non-parametric treatment effect estimation, using weighting, matching, or doubly robust methods. A growing number of simulation studies has investigated the finite sample behavior of such treatment effect estimators when the treatment is exogenous conditional on covariates, see for instance \citeasnoun{Froe00a}, \citeasnoun{Zh04}, \citeasnoun{LuncefordDavidian2004}, \citeasnoun{BuDNMC09}, \citeasnoun{HuLeWu13}, and \citeasnoun{FrHuWi14}. We consider such methods to estimate the LATE based on estimates of the instrument propensity score. We also vary the degree of flexibility of the estimators and implement parametric, semi-parametric, and non-parametric approaches to compute the reduced form and first stage effects in the numerators and denominators of Equations (ref) and (ref).
\citeasnoun{SmithTodd00}, among others, regard treatment effect estimators as weighted differences in outcomes. We apply this definition to the Wald formula and express the LATE as:
$n$ denotes the size of an i.i.d.\ sample of realizations of $\{Y_i,D_i,Z_i,X_i\}$ with observation $i \in {1,...,n}$. $n_1=\sum_{i=1}^n Z_i$ is the size of the subsample of those with $Z_i=1$, $n_0=n-n_1$, and $\hat{w}_i$ are weights that may depend on $X_i$ or $\hat{p}(x)$, an estimate of the propensity score $p(x)$. Next, we discuss different methods of estimating $\hat{p}(x)$ and $\hat{w}_i$.
We consider two different approaches to balance the covariates across groups for units with $Z_i=0$ and $Z_i=1$. One is to directly control for covariates $X_i$, but some LATE estimators alternatively control for estimates of $p(x)$, which is motivated by the propensity score's balancing properties discussed in \citeasnoun{rosenbaum1983}. Their results imply that $p(x)$ is capable of equalizing the covariate distributions across instrument states, such that the instrument is conditionally independent of potential outcomes and treatments given the propensity score whenever independence holds conditional on the covariates. A practical advantage of controlling for the propensity score (rather than a vector of covariates) is that it is one-dimensional and thus, avoids the curse of dimensionality.
We compute $\hat{p}(x)$ in three different ways. Firstly, we specify a probit model to estimate the conditional probability $P(Z_i=1 \vert X_i=x_i)$ by
where $\tilde{\beta}_{ML}$ denotes the estimated probit coefficients based on maximum likelihood and $\Phi(x_i^T\tilde{\beta}_{ML})$ is the cumulative distribution function of the standard normal distribution evaluated at $X_i^T\tilde{\beta}_{ML}$.
Secondly, we apply the covariate balancing propensity score (CBPS) method by \citeasnoun{ImaiRatkovic2014} to compute $\hat{p}(x)$. This methodology maximizes covariate balancing when predicting treatment assignment using the generalized method-of-moments (GMM) framework. \citeasnoun{ImaiRatkovic2014} show that the CBPS method is robust to mild misspecifications of the propensity score model, which is estimated by the following expression:
where $\tilde{\beta}_{GMM}$ are coefficients estimated by GMM and $\Lambda(x_i^T\tilde{\beta}_{GMM})$ is the cumulative distribution function of the standard logistic distribution evaluated at $x_i^T\tilde{\beta}_{GMM}$. We use the overidentified version of CBPS, with more moment conditions (based on the covariate balancing condition and the score of a logit model) than coefficients $\beta_{GMM}$, which are estimated by continuously updated GMM estimation:
$\bar{g}_{\beta}(Z,X)$ is the sample mean of the moment conditions and $\Sigma_{\beta}(Z,X)$ is a consistent variance estimator, described in more detail in Chapter 2.2 of \citeasnoun{ImaiRatkovic2014}.
Our third estimator of the instrument propensity score is fully non-parametric and based on kernel regression:
Equation (ref) corresponds to the Nadaraya-Watson (local constant) kernel estimator, where $K$ denotes the Epanechnikov kernel and bandwidth $h$ is chosen by least-squares cross-validation, i.e., by minimizing the least squares cross validation error w.r.t.\ $h$, see \citeasnoun{LiRacine06}. As an alternative to using $\hat{p}(x)^{lc}$ as weighting function, we also apply the Nadaraya-Watson estimator for estimating the outcome and treatment models in Equation (ref), see our discussion on non-parametric estimation methods in Chapter (ref).
A practically relevant issue of treatment effect methods is thin or lacking common support (or overlap) in the propensity score, which may compromise estimation due to a non-comparability across groups, see the discussions in \citeasnoun{Imbens03}, \citeasnoun{ImWo08}, and \citeasnoun{LeSt19}. If specific propensity score values among one group are either very rare (thin common support) or absent (lack of common support) among the opposite group, as it may occur close to the boundaries of the propensity score, some units may receive a very large weight $\hat{w}_i$ in LATE estimation as provided in Equation (ref). In the case of thin common support, these observations could dominate the estimator of the LATE which may potentially entail an explosion of the variance. In the case of lacking common support, this even introduces asymptotic bias by giving a large weight to observations that are not comparable to observations in the opposite group in terms of the propensity score.
\citeasnoun{HuLeWu13} and \citeasnoun{Bodory20} consider a trimming procedure to tackle common support issues in the sample also discussed in \citeasnoun{Imbens03}, which is asymptotically unbiased if common support holds asymptotically. It is based on setting the weights of those observations to zero whose relative share of all weights within either instrument state in Equation (ref) exceeds a particular threshold value in % (denoted by $t$):
We set the threshold $t$ to 5% and trim observations based on the weights of normalized IPW, see ((ref)), irrespective of the LATE estimator considered. This changes (in finite samples) the target parameter due to discarding observations with extreme weights, but ensures common support prior to estimation. Note that our bootstrap variance estimators discussed in Section (ref) account for the stochastic nature of trimming.
Inverse probability weighting (IPW) reweighs (instrument) group-specific outcomes such that the distribution of the covariates in the total population is matched, see \citeasnoun{Hirano+00} for a more detailed discussion. We consider a normalized IPW estimator in our simulations, which performed well in several simulation studies on conditionally exogenous treatments, see for instance \citeasnoun{HuLeWu13} and \citeasnoun{BuDNMC09}. The IPW-based LATE estimator corresponds to
\sloppy The normalizations $\sum_{j=1}^{n}\frac{ z_j}{ \hat{p}(x_{j})}$ and $\sum_{j=1}^{n}\frac{1-z_j}{1-\hat{p}(x_{j})}$ ensure that the weights in curly brackets add up to one. It is easy to see that ((ref)) corresponds to ((ref)) when setting $\hat{w}_i$ in the latter to $z_i n_1 \left\{\frac{\frac{1}{ \hat{p}(x_{i})}}{\sum_{j=1}^{n}\frac{ z_j}{ \hat{p}(x_{j})}} \right\}+ (1-z_i) n_0 \left\{\frac{\frac{1}{1-\hat{p}(x_{i})}}{\sum_{j=1}^{n}\frac{1-z_j}{1-\hat{p}(x_{j})}} \right\}$. IPW possesses the desirable property that it can attain the semiparametric efficiency bound (implying the smallest possible asymptotic variance) derived by \citeasnoun{Ha98}, if the propensity score is estimated non-parametrically (while this is generally not the case for parametric propensity scores). Furthermore, it is computationally inexpensive and easy to implement. However, evidence in the treatment effect literature suggests that IPW also has an important drawback: at the boundaries of the support of the propensity score, estimation may be unstable and the variance may explode in finite samples, see \citeasnoun{Froe00a} and \citeasnoun{KhTa07}.
Doubly robust (DR) estimation combines IPW with outcome regression. It reweighs outcome models for different instrument states by the inverse of the propensity scores. Denoting the conditional mean outcomes in the presence and absence of the instrument by $\mu_z^y(x):=E[Y_i \vert Z=z_i,X_i=x_i]$ and $\mu_z^d(x):=E[D_i \vert Z_i=z_i,X_i=x_i]$, the DR LATE estimator corresponds to
For non-binary outcomes, we run OLS regression to compute $\hat{\mu}_z^y(x)=x_i^T\hat{\beta}_{z,OLS}$. For binary outcome and treatment variables, we apply probit regression to compute $\hat{\mu}_z(x)=\Phi(x_i^T\hat{\beta}_{z,ML})$. The coefficients $\beta_z$ are estimated in the subgroups with $Z_i\in\{0,1\}$. Differently to IPW, which exclusively relies on reweighing by the propensity score, the DR estimator remains consistent even if either $\hat{p}(x)$ or $\hat{\mu}_z(x)$ is misspecified, as it makes use of both, the treatment and outcome models. If both are correctly specified, the DR estimator is semi-parametrically efficient, as discussed in \citeasnoun{RobinsRotnitzkyZhao1994}.
Matching is based on assigning (matching) to each observation in one instrument state one or more units in the other instrument state with comparable covariates, in order to estimate the LATE based on the ratio of average differences in the outcome and the treatment across units with and without instrument in the matched sample. We implement multiple variants of two types of matching methods, pair and radius matching, to estimate $\theta$.
Pair (or one-to-one) matching with replacement (implying that an observation may be matched several times) as discussed in \citeasnoun{Ru73a} matches to each reference observation exactly the observation with the most similar covariates in the opposite instrument state. This implies the following weights in Equation (ref):
$\varpi_{i,j}$ is the weight of the outcome (or treatment) of observation $j$ in one instrument group (e.g., $Z_j=0$) when matched to unit $i$ in the opposite group (e.g., $Z_i=1$), with $Z_k=1-Z_i$. $\mathbb{I}\{\cdot \}$ is the indicator function, which is one if its argument is true and zero otherwise. $\hat{f}(\cdot)$ is a function of the difference in covariates between observations $i$ and $j$. For example, the function could be defined as the difference in propensity score estimates of observations $i$ and $j$ in the case of propensity score matching or as a distance metric w.r.t.\ the covariate values of $i$ and $j$ like the Euclidean distance in the case of matching directly on the covariates. In pair matching, all weights are zero except for the observation $j$ with the smallest difference with reference unit $i$, which receives a weight of one. For propensity score matching, we base the weights on the distance of the one-dimensional propensity score, while for direct matching, we use a normalized Euclidean distance metric, where differences in the covariates are weighed by the inverse of the variances of $X_i$. Because only one observation is matched to each unit irrespective of the sample size and the potential availability of several suitable matches with similar covariates, pair matching is not efficient (i.e., does not attain the smallest possible variance asymptotically). On the other hand, it is likely more robust to propensity score misspecification than IPW, in particular if the misspecified propensity score model is only a monotone transformation of the true model, see for instance \citeasnoun{Zh08}, \citeasnoun{MiTc09}, \citeasnoun{Waernbaum2012}, and \citeasnoun{HuLeWu13}.
Radius matching as discussed in \citeasnoun{RosenbaumRubin1985} and \citeasnoun{DehejiaWahba99} uses all matches with propensity scores within a predefined radius around the reference unit, which trades off some bias in order to increase efficiency (or precision). This approach expectedly works relatively well if several comparable potential matches are available for a reference unit. In the simulations, we consider the radius matching algorithm of \citeasnoun{LeMiWu11}, which performed well in \citeasnoun{HuLeWu13}, who also provide details on the radius matching-related weighting function $\hat{w_i}$ in Equation (ref). The estimator combines distance-weighted radius matching, where units within the radius are weighted proportionally to the inverse of their distance to the reference unit, with a regression-based bias correction, see \citeasnoun{Ru79} and \citeasnoun{AbIm11}. For the bias correction, we apply an OLS regression adjustment for $Y$ and a probit regression adjustment for $D$ to remove small and large sample bias due to mismatches. \citeasnoun{HuLeSt2014} provide a detailed description of the estimator. As in \citeasnoun{LeMiWu11}, the radius size in our simulations is defined as a function of the distribution of distances between reference units and matches in pair matching. Namely, it is set to 3 times the maximum pair matching distance. Note that we include radius matching both with and without conditioning on the covariate `age at first birth' in addition to the propensity score to account for this influential confounder.
In our simulations, IPW, DR estimation, and matching are implemented with various degrees of flexibility in terms of parametric assumptions. We consider both semi-parametric versions based on parametric propensity score models, $\hat{p}(x)^{probit}$ and $\hat{p}(x)^{CBPS}$, as well as fully non-parametric estimators using the non-parametric propensity scores $\hat{p}(x)^{lc}$ (based on a local constant kernel regression) or when directly conditioning on $X_i$. For non-parametric DR estimation, also the conditional means of the binary treatment and binary (or non-binary) outcome $\hat{\mu}_z(x)$ are estimated by local constant (or local linear) kernel regressions.
In addition, we also consider several parametric treatment effect estimators. The first parametric approach computes the LATE by differences in the conditional mean functions $\hat{\mu}_z(x)$, which are estimated by OLS regressions for non-binary outcomes and by probit regressions for the treatment and binary outcome variables (see Section (ref)). Formally, this regression-based LATE estimator corresponds to the following expression:
Furthermore, we apply two-stage least-squares (2SLS) estimation, which was also applied by \citeasnoun{Angrist+98} for analysing the data our simulations are based on. 2SLS may be regarded as a benchmark method for instrumental variable estimation under the assumption of homogeneous treatment effects. Formally, the 2SLS estimator is given by
where $\tilde{x}_i:=(1,x_{i,1},\cdots,x_{i,K})$, $\tilde{z}_i:=(\tilde{x}_i,z_i)$, and $K$ denotes the number of covariates $X_i$. Note that in our just-identified settings with one treatment and one instrumental variable, the 2SLS estimator is numerically identical to the limited information maximum likelihood (LIML) estimator.
We analyze the performance of three further non-parametric estimation methods that do not impose any functional form assumptions on the regression functions of the outcome or the treatment.
Firstly, we apply the generalized random forest (GRF) method, a non-parametric estimator introduced by \citeasnoun{ATW18}. GRF is a variant of random forest algorithms, a machine learning approach, see for instance the discussion in \citeasnoun{LUW20} and citations therein. As described in \citeasnoun{Breiman2001}, random forests consist of averaging the predictions of many decision trees applied to different subsamples that are repeatedly drawn from the original data. In each of these samples, a decision tree partitions the space of $X_i$ into a set of rectangles and computes the fitted value of $Y_i$ as the average outcome in each of the rectangles. The partitions are chosen in a data-driven way such that the predictive performance is maximized (e.g.\ by minimizing the squared residuals based on the fitted values in each rectangle). A popular estimation algorithm for decision trees is CART (classification and regression tree), see for instance Chapter 9.2 in the textbook of \citeasnoun{hastie01statisticallearning}.
GRF shares the core features of `traditional' random forest algorithms like recursive partitioning, subsampling from the original data, and the random selection of a subset of covariates at each partitioning step. However, as a methodological twist, GRF uses a gradient-based partitioning scheme and a particular (so-called `honest') sample splitting technique (within any of the drawn sub-samples) that avoids overfitting the predictive models to the specificities of the data, see \citeasnoun{WagerAthey2018}. Using the conditional expectation function $\mu_z(X_i)$ in Section (ref) and applying the GRF to estimate the latter for the outcome and the treatment to obtain $\hat{\mu}_{z,RF}^Y(x)$ and $\hat{\mu}_{z,RF}^D(x)$ for $z\in\{0,1\}$ (where the subscript RF indicates the random forest approach), we compute the LATE as follows:
Algorithm 1 in \citeasnoun{ATW18} provides more details on the GRF method. We estimate the conditional expectations in Equation (ref) using the default options of the causal_forest function of the grf package for the statistical software R, see \citeasnoun{grf2020}.
Alternatively, we could have estimated the predictions $\hat{\mu}_{z,RF}(x)$ by standard Breiman-type random forests Breiman2001, or considered double/debiased machine learning estimators based on Neyman orthogonal scores DML18 or alternative causal forest algorithms MCF19 for estimation. Such methods would also be appropriate to evaluate the finite sample performance of LATE estimators in high-dimensional settings with many potential covariates, as they are capable of selecting control variables in a data-driven way, an interesting topic that we leave for future research.
Secondly, we use non-parametric kernel regression to estimate the conditional mean functions $\hat{\mu}_z(x)$ defined in Section (ref), see the subscript NP in the respective estimates in Equation (ref). For non-binary outcomes, $\hat{\mu}_{z,NP}^y(x)$ is estimated by local linear kernel regression, for the binary outcome and treatment variables, $\hat{\mu}_{z,NP}^y(x)$ and $\hat{\mu}_{z,NP}^d(x)$ are estimated by local constant kernel regression.
Finally, we consider a (naive) LATE estimator that is based on the mean differences of the outcome and treatment variables, respectively, across instrument states, which in contrast to the other methods does not control for the covariates. Therefore, the consistency of this approach provided in Equation (ref) generally requires that the IV assumptions hold unconditionally, i.e., without conditioning on $X$.
Table (ref) summarizes the LATE estimators analysed in our simulation study along with the corresponding conditioning sets.
Treatment effect estimation frequently relies on the non-parametric bootstrap for statistical inference Ef79,Horowitz2001. In an extensive simulation study with a conditionally exogenous treatment, \citeasnoun{Bodory20} find evidence that variance estimation of treatment effect estimators based on bootstrap procedures outperforms asymptotic variance approximations in terms of rejection and coverage probabilities in finite samples. These results even hold for matching estimators in small samples, despite the inconsistency of the non-parametric bootstrap for the (non-smooth) pair matching estimator, see the discussion in \citeasnoun{AbadieImbens06}.
For this reason, we apply the non-parametric bootstrap to estimate the standard errors of all LATE estimators. This algorithm randomly draws $B$ bootstrap samples of size $n$ (the size of a simulation sample) with replacement out of each simulation sample and estimates the LATE in every draw. Denoting the $B$ bootstrapped LATE estimators by $\hat{\theta}^b$, with $b\in\{1,2,\dots,B\}$, we estimate the standard error $\sigma$ of a LATE estimator by
In line with \citeasnoun{Bodory20}, we set $B=199$. Bootstrapping naturally accounts for heteroscedasticity as well as uncertainty due to trimming of influential observations and propensity score estimation.
Simulations often rely on randomly generated data drawn from a probability distribution that is selected by the researcher. However, the data generating processes (DGPs) of such simulations may appear somewhat arbitrary in the sense that they might be far from reflecting typical associations between variables in empirical data. To improve upon this caveat, \citeasnoun{HuLeWu13} suggest a simulation design based on empirical data, also called Empirical Monte Carlo Study (EMCS), an idea that has been subsequently applied in several papers, see for instance \citeasnoun{FrHuWi14}, \citeasnoun{HuLeMe2016}, and \citeasnoun{Bodory20}, among others. Briefly, the idea of an EMCS is to randomly draw small samples from large real data sets while relying as much as possible on the empirical associations between the variables when generating the simulation designs.
Our study follows this EMCS approach to evaluate the properties of various IV estimators of the LATE, with the aim that the simulation designs are more closely linked to real world data. However, we point out that also in an EMCS, several important choices about the simulation features are to be made by the researcher such that the DGPs are not fully determined by the data, see the caveats raised by \citeasnoun{AdSl2013}. The remainder of this section describes the implementation of our EMCS. We first present the empirical labor market data underlying our simulations and then provide the steps for generating the various simulation designs.
Our simulations are based on empirical data analysed in \citeasnoun{Angrist+98}, who aim at exploiting exogenous variation in family size to evaluate the treatment effect of fertility, defined as having at least three vs. two children, on female labor supply. This database is well suited to analyze the finite sample properties of IV estimators by means of an EMCS for several reasons. First, the data set is large, as it comprises 394840 observations and therefore easily allows one to draw many different random subsamples. Furthermore, the data contains a strong instrument that importantly affects fertility, namely twins at second birth.\footnote{There may be cases where the randomness of twin births is violated, see \citeasnoun{FGV18} for a discussion on dizygotic twinning. In our simulation study, we artificially generate random and non-random instrument assignments.} Finally, it provides demographic information on the mothers, which may be used as covariates to control for potential confounders of the instrument and the outcome.
Coming from the 1980 Census Public Use Micro Samples (PUMS), the data set contains information on young mothers aged 21 to 35, all of which gave birth to at least two children. Our analysis considers two different outcomes, the number of weeks worked within one year (with 43% zeros) and an indicator for being employed at all in that year. The binary treatment variable indicates if a mother has more than two kids (treatment is one) or two kids (treatment is zero). The binary instrumental variable is one if a mother gave birth to twins at second birth and zero otherwise. The covariates considered in our simulation include mother's age, mother's age at first birth, race, and quarter of birth.
Table (ref) reports descriptive statistics of the database, by treatment indicator (more than two kids) and the instrument (twins at second birth). The upper part presents descriptives for the two labor market outcomes `weeks worked' (in weeks) and `worked for pay' (binary). There are large differences between the outcomes of the treated and non-treated in terms of the standardized difference statistic as suggested by \citeasnoun{RosenbaumRubin1985} (the literature considers values around 20 and above as severely unbalanced). The line underneath the outcomes in Table (ref) gives details on the treatment variable. Not surprisingly, the treatment fully complies with the instrument if the latter equals one, because all mothers with twins at second birth ($Z_i=1$) necessarily have more than two children ($D_i=1$). The subsequent row of Table (ref) provides information on the instrument. It reveals that 2% of women with at least three children have twins at their second birth. Considering the covariates, the standardized differences show that mothers' characteristics are partly unbalanced across treatment states, whereas they are well balanced across instrument states, in line with a randomly assigned instrument. The randomness of the instrument is also supported by the pseudo-R2 statistic with a value of 0.2% when regressing the instrument on the covariates.
Data generating processes (DGPs) may differ in (infinitely) many dimensions. We select ten practically relevant dimensions for varying the specifications of our simulation models. These dimensions include: effect homogeneity vs.\ heterogeneity, randomness vs.\ non-randomness of the instrument, varying levels of instrument strength, binary vs.\ non-binary outcome distributions, and different sample sizes. Summary statistics of all DPGs are presented in Table (ref).
We start by assuming homogeneous treatment effects with a randomly assigned instrument and the empirically observed instrument strength. To evaluate the performance of the estimators under these conditions, we define a new population for which the true LATE is equal to zero. To this end, we drop all 3380 observations from the database who receive the instrument ($Z_i=1$). Among the remaining 391460 observations with instrument state $Z_i=0$ (no twins at second birth), there is no reduced form effect of the instrument on the outcome or first stage effect of the instrument on the treatment, such that there exists no LATE. After that, we create a pseudo-instrument and artificially assign $Z_i=1$ to those who are similar to the 3380 discarded observations in terms of observed characteristics. This similarity is determined by $1:M$ matching on the covariates without replacement. By setting $M=58$, we assign $Z_i=1$ to approximately half of the observations, see column 4 of Table (ref). In addition, we set the treatment state of everyone with $Z_i=1$ to $D_i=1$ (as in the original database) to maintain the empirically observed instrument strength. Finally, we draw small samples from our new population to compare the finite sample properties of alternative LATE estimators.
To simulate specifications with a weaker instrument, we reduce the first stage effect by lowering the impact of $Z_i$ on $D_i$. Instead of setting all observations with $Z_i=1$ to $D_i=1$, we change the treatment status from zero to one only for those with $Z_i=1$ for which the condition $D_i=\mathbbm{1}(u_i>1.25 )$ holds. $\mathbbm{1}(\cdot)$ denotes the indicator function which is one if its argument is true, otherwise it is zero, and $u_i$ is a standard normally distributed random variable. Column 7 of Table (ref) displays the first stage coefficients for the different DGPs.
The randomness of the instrument implies that the covariates are balanced across groups. To mimic a non-random assignment of the instrument, we increase the magnitude of instrument selectivity in the following way. We first estimate the propensity score $\hat{p}(1.5X_i)^{probit}$ (see Equation (ref)) using the original database. Then, we change the instrument status $Z_i$ from zero to one for observations with characteristics similar to the 3380 observations dropped from the original database (with $Z_i=1$). We obtain such similar matches by $1:M$ matching on the estimated propensity score $\hat{p}(1.5X_i)^{probit}$, with $M=22$. Next, we assign $D_i=1$ to all observations with $Z_i=1$. Based on this modified data set, we increase the selection into the instrument by discarding the best matches for the newly created observations with $Z_i=1$ among observation with $Z_i=0$. To find the best matches to be discarded, we apply $1:M$ matching on a newly estimated propensity score $\hat{p}(X_i)^{probit}$ (with the modified instrument assignments), where $M=3$. The selectivity of the instrument is provided in columns 5 and 6 in Table (ref).
To model a scenario with non-constant treatment effects, we introduce effect heterogeneity with respect to age and race as follows. We add to the existing control variables squared and cubic terms of both age variables (`age' and `age at first birth') and interact the unmodified age variables with the indicator variable for African Americans. This new set of control variables for settings with effect heterogeneity is denoted by $X_i^{het}$ for each unit $i$. We generate $Y_i$ and $D_i$ in each simulation sample according to the rules $Y_i=Y_{i,1}^d Z_i+Y_{i,0}^d (1-Z_i)$ and $D_i=D_{i,1}Z_i+D_{i,1}(1-Z_i)$. To this end, we compute the non-binary potential outcomes based on the equation $Y_{i,z}^d=X_i^{het}\hat{\beta}_{OLS}+\hat{\sigma}v_i$, where $v_i$ is a standard normally distributed random variable. $\hat{\beta}_{OLS}$ and $\hat{\sigma}$ are the coefficients and residual standard deviation of OLS regressions in subsamples by instrument state $Z_i\in\{0,1\}$ of our new population. The binary potential outcomes are computed based on $Y_{i,z}^d=\mathbbm{1}(X_i^{het}\hat{\beta}_{probit}+v_i>0)$, where $\mathbbm{1}(\cdot)$ is the indicator function and $\hat{\beta}_{probit}$ are the coefficients estimated from probit models in subsamples by instrument state of our new population. The potential treatments $D_{i,1}$ are set to one, whereas $D_{i,0}$ is computed analogously to $Y_{i,z}^d$ in the binary outcome case.
We combine these variations in the DGPs with respect to effect heterogeneity, instrument strength, and instrument selectivity with smaller and larger sample sizes of 1000 and 2000, respectively, and with binary and non-binary outcome distributions. We run 2000 simulations for the smaller and 1000 simulations for the larger samples. Table (ref) presents summary statistics of the DGPs considered in our simulation study.
This section presents results about the finite sample performance of various LATE estimators across different DGPs. We rank the estimators by their coverage rates, which are defined as the share of simulations in which the true LATE is included in the 95% confidence interval of the respective LATE estimator. We recall that the standard errors for computing those confidence intervals come from the non-parametric bootstrap, as discussed in Section (ref). For the sake of brevity, we subsequently only discuss a selection of our results, which conveys the main message of our findings. In the Appendix, we include more detailed results.
Table (ref) provides the average coverage rates and lengths of confidence intervals across all DGPs of any parametric, semi-parametric, or non-parametric LATE estimator which performs best (in terms of coverage) in at least one of the ten DGPs discussed in Section (ref). We find that only the non-parametric random forest-based LATE estimator described in Equation (ref) attains exactly the nominal coverage size of 95% on average. Furthermore, its average length of confidence intervals is the second shortest among the estimators analyzed in Table (ref), 8% larger than the average interval of the parametric regression estimator, the nominal size of which is 96%. Conditional on obtaining a decent coverage, a short confidence interval is desirable in terms of precision, as it implies a lower estimation uncertainty. Three out of the four LATE estimators whose average coverage rates come closest to 95% are non-parametric, with those of non-parametric kernel regression (94.8%) and pair matching on the covariates (94.6%) having a minor under-coverage. Also semi-parametric radius matching on the propensity score performs decent in terms of coverage rates, with the probit-based version attaining an average rate of 95.1%, and two further versions achieving 94.4% and 94.1%, respectively. Furthermore, also IPW using the CBPS method for propensity score estimation reaches a satisfactory average coverage rate of 95.5%. However, the average length of the confidence intervals of radius matching, kernel regression, pair matching, and IPW is substantially larger than that of parametric regression or of the random forest.
The coverage accuracy of the different estimation methods is related to the bias and variance of the LATE estimators, as well as the bias of the bootstrap-based standard error. Table (ref) provides details on these statistics. We find that two non-parametric methods perform best when considering the bias of LATE estimation, the standard deviation, and the mean squared error (i.e., the sum of the squared bias and the variance), as well as the bias of the standard error (relative to a LATE estimator's true standard deviation). The random forest-based LATE estimator has on average the smallest deviation from the true LATE, with its absolute bias amounting to 0.6. The non-parametric kernel regression estimator has the smallest average standard deviation among the estimators in Table (ref), amounting to 7.0. It also performs best in terms of root mean squared errors across DGPs with an average value of 7.2. When considering the median bias of the bootstrap standard errors relative to the true standard deviations of the respective LATE estimators, the inference method of the random forest-based LATE estimator performs best, with an average median bias of 1.7. The averages of the median biases of its competitors are on average at least 100% larger.
Our findings suggest that the coverage accuracy is mainly driven by a LATE estimator's bias. This is for instance the reason why the mean differences estimator (which ignores covariates), the bias of which exceeds the bias of the random forest-based LATE estimator by 440.1%, shows a poor coverage in Table (ref). Also the OLS estimator performs poorly in terms of coverage, due to its high bias, while its variance is small (results not presented but available on request). The performance of all LATE estimators by DGP is presented in Section (ref) of the Appendix. Tables (ref)-(ref) provide details on the coverage rates, biases, standard deviations, and root mean squared errors.
Table (ref) lists the best performing LATE estimators in terms of average coverage, separately for each of ten DGP features (see the rows) as well as for parametric, semi-parametric, and non-parametric methods (see the columns). The results suggest that semi-parametric and non-parametric estimators come closest to the nominal coverage rate of 95%. The radius matching algorithm of \citeasnoun{LeMiWu11} (with or without controlling for a covariate in addition to the propensity score) most frequently performs best both among the semi-parametric LATE estimators (in 70% of cases) and overall (in 50% of cases). Radius matching achieves the best average coverage in settings with effect homogeneity, a strong instrument, non-binary and binary outcomes, and under a larger sample size. For specifications with standard and strong selection into the instrument, the non-parametric random forest-based estimator is closest to the nominal size. Considering scenarios with effect heterogeneity, a weaker instrument, and a small sample size, the best performers are LATE estimators based on 2SLS, DR, and non-parametric regression, respectively.
The performance of the best performing LATE estimators across the ten DGP features is presented in Section (ref) of the Appendix. Tables (ref)-(ref) provide information on the average coverage rates, biases, standard deviations, and root mean squared errors.
This paper presented a simulation study based on empirical labor market data to investigate the finite sample properties of a range of point estimators of the local average treatment effect (LATE) when controlling for a fixed (and low-dimensional) set of covariates. The structure of these estimators is inspired by the Wald estimator, consisting of the ratio of the estimated reduced form effect of the instrument on the outcome and the estimated first stage effect of the instrument on the treatment. Furthermore, we applied the non-parametric bootstrap to estimate the standard errors and the 95% confidence intervals of the LATE estimators. We find that among the LATE estimators considered, non-parametric kernel regression has the smallest average root mean squared error across the different simulations, closely followed by the random forest-based approach, which has the lowest average absolute bias. The random forest method also performs very competitive in terms of average coverage rates, while at the same time having relatively narrow confidence intervals, which is attractive in terms of precision. Specific versions of semi-parametric radius matching on the propensity score, nonparametric kernel regression, inverse probability weighting, and pair matching on the covariates perform decently in terms of coverage, too, but have substantially wider confidence intervals.
{ \setcounter{equation}{0}