Variable Selection for Causal Inference via Outcome-Adaptive Random Forest

Variable Selection for Causal Inference via Outcome-Adaptive Random Forest

\thispagestyle{empty}

\setcounter{page}{1}

Introduction

In the causal inference literature, we can classify data into two categories. The one is data from a randomized controlled trial where the researcher or practitioner has full control of the selection process. The counterexample is data from a so-called observational study. In such a setting there are confounding variables that influence both, the outcome and the probability of treatment. To construct unbiased treatment effect estimates from observational studies, propensity score (PS) methods are an increasingly popular tool to control for confounding rosenbaum1983central. One model-based approach, the inverse probability of treatment weighting (IPTW), to directly adjust for the confounding bias using propensity scores was proposed by hirano2001estimation. Estimating the propensity score can be seen as a classification problem where one seeks to have a good prediction of the assignment probability given covariates, regardless of the functional form of the distribution of the probabilities. Besides logistic regression, non-parametric methods such as random forests lee2010improving, westreich2010propensity, zhao2016propensity, neural-networks, and support vector machines westreich2010propensity have been proposed to estimate the propensity score. An interesting question is which variables should be included to estimate the propensity score. Common suggestions are, to include all pre-treatment variables that influence the treatment while excluding variables that do not affect the treatment (this can be variables that do not influence any dependent variable - they are spurious, but also variables that are only predictive on the outcome). Following this rule, we would include confounding variables since they influence the treatment and the outcome as well as variables that only predict the treatment and not the outcome but exclude variables that are only predictive of the outcome.

We give an example of the different relationships between the covariates, the outcome, and the treatment in Figure (ref). We denote $X_t$ for covariates that predict only treatment but not the outcome, and $X_o$ that predict the outcome but not the treatment probability. Of special interest are the confounding covariates $X_c$ that we need to take into account to get an unbiased estimate of the average treatment effect (ATE) while $X_s$ are spurious covariates that are uncorrelated to both, treatment and outcome. Let us illustrate the role of the variables using the vaccination against COVID 19 as an example. The treatment variable ($D$) is whether a person is vaccinated or not while the outcome variable ($Y$) is the individual probability of severe symptoms. Variable $X_t$ could be the industry sector a person is working in since it influences the probability of being vaccinated but has no influence on symptoms. $X_o$ could be whether a person smokes or not, which has an influence on the degree of symptoms but does not determine the vaccination probability. $X_c$ might be the age of a person, which is associated with vaccination and symptoms while the variable body height might be unrelated and is classified as $X_s$.

The rational behind the classification of covariates is to perform variable selection when estimating the propensity score. First, and most important to ensure unbiased treatment effects in observational studies is to achieve unconfoundedness. In theory, it is sufficient to only use covariates $X_c$ for the propensity score since the main analytic goal is to eliminate bias. This implies that we do not aim to explain the treatment assignment mechanism with high accuracy. If this would be the goal, variables $X_t$ would need to be included. However, doing so could limit the overlap assumption. These are the two main reasons why we want to exclude covariates $X_t$ in the propensity score model. Last, even if the dependency of $X_o$ and $D$ is zero in the true data generating process, the finite sample bias can be reduced if variables $X_o$ are included in the propensity score model. The finite sample bias arises due to random confounding when the sample size is small.

In a recent paper by shortreed2017outcome, the authors suggest a different approach to get unbiased treatment effects from observational studies with a focus on decreasing the variance. Their proposed outcome-adaptive lasso (OAL) approach only selects features in the propensity score estimation that have a relationship with the outcome ($X_c, X_o$) but exclude variables that are predictive on the treatment ($X_t$) as well as spurious variables ($X_s$). To do so, they first find covariates that predict the outcome by regressing the outcome on the treatment and all covariates using a linear model. In the second step, the estimation of the propensity score, they use the lasso with an additional penalty term. The penalty contains individual weights for each covariate based on the coefficients of the covariates from the first step. The result is that the lasso excludes covariates that predict the treatment but are not related to the outcome as well as spurious covariates. Consider $p$ covariates, denoted $X_j$ for $j=1,...,p$, then the OAL estimator is defined as:

where $\omega_j = |\tilde{\beta}_j|^{-\gamma}$ s.t. $\gamma>1$. The vector $\tilde{\beta}_j$ refers to the coefficient estimates from regressing the covariates on the outcome, conditioning on the treatment.

An important limitation of the outcome-adaptive lasso is that this approach is restricted to parametric models. Both, the outcome model and the propensity score model have to be correctly specified. However, to control for selection bias we would like to have as many pre-treatment characteristics as possible to condition on them. This requirement lets new datasets easily become high-dimensional. In such settings, we face two potential challenges. First, the outcome and the treatment variable might not dependent linearly on the covariates. There can be interactions between variables and complex structures. Second, even if we could include such interactions, facing many covariates leads to the problem of which covariates to include in the outcome and propensity score model? It might be the case that only a few characteristics are important - the question is which are they? Such settings call for a non-parametric method that uses regularization. Keeping the idea of outcome-adaptive regularization but accounting for non-linearity and high-dimensionality, we propose the outcome-adaptive random forest (OARF). First, we estimate a modified and standardized variable importance score using a random forest used as the penalty weight. Second, we use the modified variable importance to regularize a random forest that learns the propensity score. To do so we penalize the gain at each split and propose the use of an initial feature space using the modified variable importance. This approach allows replacing both linear models, the regression (which we will refer to as OLS), and the lasso, with a random forest to estimate the ATE when the functions are non-linear or high-dimensional. We also make use of sample splitting to avoid overfitting and apply cross-fitting to restore efficiency. Our approach of a regularized random forest further selects only those covariates in the propensity score that are predictive of the outcome and excludes spurious variables. The OARF is designed to allow for categorical variables and is robust to different amounts of levels between categorical or continuous variables. The result is a flexible non-parametric method that allows unbiased estimation of the ATE while decreasing the variance. Our approach is also fast in computation (about 30 seconds for 2000 observations and 20 covariates).

Illustrating the outcome-adaptive estimation

Let us demonstrate the outcome-adaptive approach in a simple example. Assume the true outcome model as in equation (ref). $Y$ depends on the treatment $D$ and linearly on three covariates ($X_1, X_2, X_3$). We also generate a propensity score model that includes variables $X_1, X_2, X_5$, and $X_6$ (see the Appendix for how the function is generated). The covariates are generated from a multivariate normal distribution and are independent. We set the sample size to $N = 1000$ and generate $p=20$ covariates where only the first three are dependent on the outcome $Y$. Using the notation from Figure (ref) we have the following structure: $X_1, X_2 \in X_c$, $X_3 \in X_o$ and $X_5, X_6 \in X_t$. Using a linear model (OLS), as in shortreed2017outcome, we want variables $X_1$ to $X_3$ to have the highest coefficients and penalize all others. Hence they should have a very small estimated coefficient. We run 500 Monte Carlo simulations and predict the treatment effect using the inverse probability of treatment weighted (IPTW) estimator lunceford2004stratification, hernan2006estimating. We use a logit model to estimate the propensity score, using all covariates (the full model) and only using covariates that are confounding or predict the outcome (target model). Figure (ref) shows the estimated treatment effect using both models as well as the standardized and absolute coefficients from the OLS model. The target model (which can be seen as an oracle model) shows a smaller variance around the true treatment effect of $0.5$ (indicated with the horizontal line). We also see that the OLS correctly assigns the highest coefficients to the variables that are predictors of the outcome while all other variables get smaller coefficients.

The suggested approach by shortreed2017outcome to select only variables for the propensity score model that are predictive for the outcome works as expected if the underlying outcome and propensity score functions are linear. As soon as we introduce a more complicated function, the selection process in step 1 is biased. To demonstrate this, let us now assume the following true outcome model:

Figure (ref) shows the coefficients estimated by the OLS and the random forest, respectively. Values from the linear model are absolute coefficients while values from the random forest are based on the impurity measure which we explain below. The values for both methods are standardized between zero and one to make them comparable. Say we know that the model should select variables one to three. We draw a line indicating the lowest value from the variables that should be selected. If this would be a threshold (for selection or penalization), clearly the OLS estimates higher coefficients for at least 12 additional variables (besides the two with the highest coefficients). In contrast, the random forest does find the correct variables and assigns much lower importance to all other variables. Using the same threshold method on the random forest, all spurious variables would not be selected or at least heavier penalized compared to the OLS. The problem gets more severe if we not only allow for interaction terms but non-linear structures (say, by including quadratic and trigonometric functions into the data generating process).

Method

Selecting variables using non-parametric models is not straightforward since we do not directly estimate coefficients for each feature like in the OLS. Using a random forest, we can instead create a variable importance measure to find the most predictive variables from the outcome model. Since we want to use a random forest for both steps, the outcome model and the propensity score, we have to replace the adaptive lasso by penalizing the tree-building mechanism. In this section, we show how to best estimate the variable importance measure from the outcome model and state the importance of the initial feature space. We describe in detail how the information about the covariates is then used to penalize the random forest that estimates the propensity score. This leads to a regularized version of the random forest that shrinks penalized variables to zero.

To ensure unbiased effects from a causal parameter the following assumptions from the potential outcome framework are required: Each observation has two potential outcomes, $Y^1$ if treated and $Y^0$ if not. We denote treatment by the binary indicator $D \in \{0;1\}$ and observed covariates $X \in \mathbb{R}^{^p}$. See, for example rubin1980randomization. First, ignorability: $\left(Y_{i}^{1}, Y_{i}^{0}\right) \mathrel{\text{\scalebox{1.07}{$\perp\mkern-10mu\perp$}}} D_{i}|X_{i}$. It states that the treatment assignment is independent of the two potential outcomes. Second, the stable unit treatment value assumption (SUTVA), $Y_i = Y_i^0 + D_i (Y_i^1-Y_i^0)$, ensures that there is no interference, no spillover effects, and no hidden variation. This means that the treatment status for individual $i$ does not affect the potential outcomes of individual $j$. The third assumptions describes the propensity score: $P(D_i=1|X_i=x) \overset{\mathrm{def}}{=} e(x)$, which needs to be bounded away from zero and one: $\forall x \in supp(X_i), \quad 0 < P(D_i=1|X_i=x) < 1,$. The last assumption states that the covariates are not affected by the treatment: $X_i^1 = X_i^0$.

We are interested in estimating the ATE ($\theta$), assuming a partially linear model that takes the following form:

Variable Importance in Random Forest

Building a tree based on the classification and regression tree (CART) algorithm works as follows: A subset of the feature space $\mathcal{X}$ is selected at each node $t$ in the tree. Internal nodes $t$ are labelled with a split $s_t = (X_j < c)$, which creates at least two further subsets or children $t_L$ and $t_R$ ($L$ and $R$ refer to left and right, as in a tree). This procedure is repeated until we reach terminal nodes (or leaves) that are labeled with the best prediction of the outcome variable. In the regression case, this would be the mean value $\bar{y}_t$. The predicted output for a new instance is the mean value of the leaf reached by the instance when it is propagated through the tree. Using a recursive procedure at every note $t$, a tree identifies the split $s_t = s^*$ for which the partition of the sample into $t_L$ and $t_R$ maximizes the reduction of

Let $\mathbb{D}$ be the set of observations for a specific decision or split. For the first step regression setting, we define the cost function as $\operatorname{cost}(\mathbb{D}) = \sum_{i \in \mathbb{D}}(y_i - \bar{y})^2$, where $\bar{y} = |\mathbb{D}|^{-1}\sum_{i \in \mathbb{D}}y_i$ is the mean of the outcome variable in the specified set or region. Maximising the decrease of the variance within each leaf can be seen as making the leafs pure in terms of the outcome values. Hence, the cost function is called the impurity measure.

The impurity change (or gain) through a split can be used to estimate the importance of a variable by evaluating each cost function given a specific feature $j$ on which the split $s^*$ is based. Hence, we can define the gain in terms of a feature $j$ instead of split $s$: $\Delta(j,t) \overset{\mathrm{def}}= \Delta(s,t)$. The global importance value is given by accumulating the gain over a feature, $\Delta(j)=\sum_{t \in \mathbb{S}_{j}} \Delta(j, t)$ where $\mathbb{S}_j$ represents all the splitting points used in a single tree for the $j$-th feature. This is because a feature $X_j$ can be used multiple times during the recursive growth of a tree. Since a random forest consists of $B$ such trees, the importance value ($Imp_j$) is just an average over all trees:

nembrini2018revival find that the impurity importance can be explained by two parts: The impurity reduction directly related to the true importance and a part of impurity reduction that is highly based on the structure of the feature (i.e. a dummy variable with only two levels vs. a continuous variable). The latter component introduces a bias in the impurity measure. To correct for the structure of different features, the authors propose to extend the feature space dimension $p$ to $2p$. If the selected variable is within $\{1,..., p\}$ then the variable is used as usual. If the variable is instead in the $\{p+1,..., 2p\}$ set, the variable values are permuted. This means that the levels for each observation are reordered such that each observation has a different level as before. If the feature is reordered, the importance value adds negatively to the final measure and positively if the feature is untouched. This procedure, called actual impurity reduction (AIR), allows to de-bias the total importance value by controlling for the structure of each feature.

There are other measures for variable importance like permutation importance. To calculate the permutation importance of a feature, the prediction performance is calculated for observations that are not included in the bootstrap data (the so-called out-of-bag (OOB) observations). The values of the variable are then randomly permuted for each observation. Calculating the OOB error again with the permuted feature indicates the importance of the variable. The more the prediction error changes, the more important is the feature. The permutation performance is calculated for each feature and averaged over all trees. See nembrini2018revival for an overview and comparison of different importance measures. We note that the permutation method is computationally expensive when compared with the AIR method. This is because the former method relies on OOB predictions for each feature. In terms of robustness concerning the different amounts of levels between variables, they perform similarly.

In our setting, we have at least one binary variable, the treatment assignment, which has a limited range of splitting values compared to continuous variables. Therefore, it is important to take the different structures into account by applying the AIR measure (available in the ranger package as 'impurity corrected').

Penalization parameter

The penalization parameter $\lambda_j$ should depend on the predictive power from the covariates on the outcome. Our measure of predictive power is the importance score $Imp^*_j$ which can be included as proposed by deng2013gene:

where $\lambda_0$ is a general penalization parameter and $\gamma$ is a weight parameter that determines the proportion of general and specific feature penalization. Next, we define the normalized importance score ($Imp^*_j$) as

The new normalized importance score is scaled within the interval $[0,1]$ (see the Appendix for the proof). Equation (ref) sets all negative importance scores from the AIR measure to zero. Since we only want to rely on the importance values obtained from the outcome model, we set $\lambda_0=0$ and $\gamma = 1$. This allows for the heaviest penalization based on outcome-related covariates, namely

Creating the feature space

Next, we want to use the information obtained in the first step and only use variables that have a high importance on the outcome to estimate the propensity score. Again, we want to use a random forest since the propensity score function can have a similar non-linear structure as the outcome model. To guide the feature selection, deng2012feature introduce a (guided) regularized random forest (RRF) by proposing to weight the gains of the splits during the recursive procedure. As a result of the random feature selection at each split, after $m$ splits, only a subset $\mathbb{F}$ of features are included in the tree. To limit the feature space the idea is to exclude variables not belonging to $\mathbb{F}$, unless their importance is substantially larger than the maximum of the gain for features already included. deng2013gene define the regularized gain as

where $\lambda_j \in (0,1]$ is the penalty coefficient that controls the gain for each feature $j$ if this feature was not previously used for a split. Originally, the feature space $\mathbb{F}$ is an empty set at the root note in the first tree. Only if a feature adds enough predictive information it is included. Based on equation (ref), the smaller the value for $\lambda_j$, the higher the penalty and hence the less likely it is for feature $j$ to be included in the subset. In our setting, we want to include features that may not be that predictive of the treatment but of the outcome. We also want to make sure that the important features are used in the splitting process (at least with a higher probability). Therefore, we already include features in $\mathbb{F}$ that fulfill the following criterion:

If the importance score $X_j$ is at least the mean over all importance scores, we include the variable in the initial feature space and drop values that are zero in $\mathbb{F}$. Figure (ref) illustrates the guided feature selection process. We start with a non-empty set (here variables 1 and 2 are included). The first split is based on variable $3$. If this variable is not in the feature space, the gain is multiplied by $\lambda_3$. If the penalized gain is higher than the gain from the parent the feature is included. As an illustration, we always move from top to bottom and from left to right. The next split is then on $X_4$, again if the penalized gain difference is positive, the feature is included in $\mathbb{F}$. Next, a split on $X_1$ is made, since the variable is already in the feature space the gain is not penalized. Still, it has to be higher than the gain from the parent node to keep the split. Building the first tree, we end up with a feature space containing four variables. The information of the feature space is now used to build further trees. In tree 2, we start with the initial features space as extracted from tree 1. Note that this procedure does not allow to grow trees in parallel since each tree needs the information from the former tree about the feature space to determine if the gain from variables should be penalized or not.

The covariate selection process does not depend on how often a variable can be split and hence if the variable is continuous or, for example, binary. It is sufficient if the gain from one split is above the threshold. Different from calculating the variable importance using the impurity measure we do not average over all split within a tree. This allows the covariates to be of any form, such as binary, categorical or continuous.

Sample splitting

To avoid overfitting, which can easily happen when using flexible methods such as random forest, we make use of sample splitting and cross-fitting. First, we split our sample into two equal parts, $I_A$ denotes the auxiliary sample and $I_E$ is the estimation sample. We first use the subset $I_A$ to train the propensity score function and $I_E$ to estimate the treatment effect using the predicted propensity score in the IPTW estimator. Let us denote the resulting estimator as $\hat{\theta}(I_A, I_E)$. Now we switch the roles of the auxiliary and estimation sample to obtain a second estimator, called $\hat{\theta}(I_E, I_A)$. Since both estimators are estimated on only a subset of the observations there might be an efficiency loss. Cross-fitting, which was recently introduced by chernozhukov2018double, aims to restore efficiency by simply averaging the two estimates:

This approach generalizes to $K$ folds where $I_A$ contains $K-1$ folds and $I_E$ the remaining fold. Similar to cross-validation, each fold is used to estimate the ATE by iteratively looping through the folds. The final estimator is the average over the $K$ estimators.

We use the full sample to get the variable importance from the first step. This is especially helpful if the sample size is small. When using a logit model or the lasso as a benchmark, we also use the full sample and estimate the final treatment effect in one step.

Simulation study

To evaluate the performance of our OARF method in more detail, we consider different data generating processes (DGP's) and consider the following methods for comparison and benchmarking: The OAL method by shortreed2017outcome and two generalized linear models. The first one uses all covariates in a logit model (Lo full) while the second one only uses target variables ($X_c, X_o$) to estimate the propensity score (Lo targ). We use the same benchmarks for the random forest, denoted by RF full and RF targ. We also use the regularized random forest (\textbf{RRF}) which only sets a penalty based on the first step variable importance but does not make use of an initial feature set, as proposed here by the OARF. We use the following order of variables when we look at the variable selection plots: $X_c, X_o, X_t$. If the amount of the variables are set to two, then the first two are confounders, variable three and four are regressors on the outcome and five and six are regressors on the treatment. We use the $\texttt{ranger}$ package in $\texttt{R}$ for all estimations based on the RF and the $\texttt{RRF}$ package by deng2013guided in part for the OARF. The OAL approach is based on the replication file from shortreed2017outcome. The tuning parameter $\lambda$ and $\gamma$ for the OAL method are selected using a weighted absolute mean difference (wAMD) which we describe in the Appendix.

First, we consider two linear settings and generate a binary treatment, $D$, from a Bernoulli distribution with $\textit{logit}\{P(D=1)\}=\sum_{j=1}^{p} \nu_{j} X_{j}$ and the continuous outcome variable $Y$ as $Y = \theta D + \sum_{j=1}^{p} \beta_{j} X_{j}+\varepsilon,$ where $\varepsilon \sim N(0,1)$ and $\theta = 0.5$. The two settings differ in the strength of the confounding effect. Setting 1 sets $\beta = (0.6, 0.6, 0.6, 0.6,0,0,0,...,0)$ and $\nu = (1,1,0,0,1,1,0,...,0)$, and setting 2 sets $\beta = (0.6, 0.6, 0.6, 0.6,0,0,0,...,0)$ and $\nu = (0.4,0.4,0,0,1,1,0,...,0)$. Setting 2 hence has a weaker confounding relationship than setting 1. These two settings are identical to the one used in shortreed2017outcome. Setting 3 aims to have a non-linear relationship between the covariates and the outcome but the same linear structure for the propensity score as in setting 1. Setting 3 is generated as follows: $Y = \theta D + 0.8(X_1 \otimes X_2) + 0.8(X_3 \otimes X_4) +\varepsilon$; $\textit{logit}\{P(D=1)\}=\sum_{j=1}^{p} \nu_{j} X_{j}$, with $\nu = (1,1,0,0,1,1,0,...,0)$. In all three setting we set $N=500$ and $p = 20$. We find that our OARF performs similar to the benchmark OAL method. The random forest using all covariates and the one that uses the penalization weights from the first step perform worse and are clearly biased. The reason for the bias might be that the full random forest has a higher chance to select variables $X_t$ and $X_s$ than OARF. We do see an improvement in terms of bias when using the RRF, which selects fewer of the above mentioned variables due to the penalization weights. We also find that if we weaken the confounding relationship, all methods have a smaller variance and the RF methods a smaller bias. Boxplots illustrating the IPTW estimator using different methods are shown in Figure (ref). Selected covariates from the propensity score model are shown in Figure (ref). In both settings, the OARF approach only selects the desired features. This is similar to the OAL method. For comparison, we show that the full RF selects all features and the RRF a higher proportion of all variables while always selecting the desired variables.

In setting 3, illustrated in Figure (ref), all approaches are unbiased since the propensity score function depends linearly on the covariates. Only the outcome model is non-linear which is why the variance is higher in the linear models compared to the random forest. Next to the treatment effect estimates, the selected variables are illustrated. The OAL model fails to select the correct features as was expected based on the coefficient values from Figure (ref). The RF without regularization uses all variables while the RRF always selects the correct four features but often (in about 80%) all other variables. Only the OARF selects the correct features and drops the covariates that are not of interest.

Next, we generate data processes where both functions are non-linear (settings 4 to 10). These are the settings where we would expect the OARF to perform superior. A complete list of all the DGP's is shown in Table (ref). In those settings, we set the sample size to $N = 1000$. Our results are shown for $p =20$ to allow a noticeable visualization. Varying the number of covariates to 50, 100, and 200 does not change the results of the ATE estimates nor the selection of the correct covariates. In setting 4 to 7, we keep the function on the outcome model and only change the propensity score function. Setting 1 to 7 sets the amount of variables for $X_c, X_o$, and $X_t$ to 2 while the amount of spurious variables $X_s = p-(X_c + X_o + X_t)$. Setting 8 to 10 uses $X_c = 6, X_o = X_t =2$ covariates and again the remaining set for $X_s$. The favourable covariates to select are $X_1$ to $X_4$ for settings 1 to 7 and $X_1$ to $X_8$ for setting 8 to 10. Figure (ref) shows ATE estimates for settings where all functions are non-linear and potentially complicated while Figure (ref) shows results for similar functions but with more depending covariates. The proportion of selected covariates over all simulations are illustrated in Figure (ref) and Figure (ref). Overall, we find that the OARF performs best and is close to the RF that only uses $X_c$ and $X_o$ (RF targ). Setting 4 and 5 show that all methods are biased while the OARF is closed to the true ATE. In setting 4, all methods are downward biased, which leads to a negative estimate using the OAL. Only the RF approaches estimate the correct sign of the treatment effect. In setting 6, the linear methods are slightly upward biased while the RF approaches show no bias. Setting 7 shows a similar effect of bias where only the OARF estimates unbiased treatment effects. Using more covariates allows making the functions even more complicated since the flexibility can be increased. Setting 8 and 10 again show some bias for the linear methods as well as a higher variance compared to the RF approaches. Setting 9 is comparable to setting 4 in the sense that all methods are downward biased. In this setting, even the OARF is biased and shows a higher variance. We notice that in some settings there is not so much difference between the OARF and the full RF. The advantage of the OARF remains since it achieves the same accuracy using fewer variables as the full RF. In all settings that have at least one non-linear function, only the OARF selects the correct features with high precision. The full RF selects all covariates and the RRF uses unnecessary covariates in about 80% of the cases.

Introducing a correlation between the covariates increases the bias and variance. We investigate the effect in a linear setting with moderate and heavy correlation. The results are shown in Figure (ref). If the correlation is moderate ($\rho =0.2$) the OARF is still unbiased while all other methods show a slight bias and an increase in variance. It is reassuring that the random forest can find the correct importance score even with moderate correlation among the variables. If we introduce a heavy correlation ($\rho =0.5$) also the random forest approaches are biased while the variance in the linear models increases heavily. The OARF is still closest to the target method. The variable selection is still correct as illustrated in Figure (ref). Last, we show results when the outcome function is more complicated and non-linear. Boxplots that illustrate the variance over 500 repetitions are shown in Figure (ref) while the corresponding variable coverage rates are shown in Figure (ref). Using a more complex outcome function increases the bias for the linear methods. Even the RRF shows a high bias in setting 13 and 14. The OARF is closest to the oracle RF. Setting 15 does show an equally small bias among all methods, again with the OARF and the oracle RF closest to the true treatment effect. The coverage rates regarding the selected variables show a similar picture as for all other settings. The vanilla RF does select all variables, the RRF selects the four correct variables in almost 100% of the cases but also all other variables in more than 75% of the cases. The OAL is not able to select the correct variables (e.g. variable 4 is selected in only 12.5% of the cases. For setting 13 and 14, the OARF selects the desired variables in 100% of the cases and drops all other variables with the same accuracy. Only in setting 15, one variable ($X_2$) is only selected around 37% of the time while all other variables are selected as desired.

In Table (ref) we show coverage rates of 95% confidence intervals and the width of the interval (in parentheses). Confidence intervals for IPW were constructed using a percentile‐based nonparametric bootstrap. For the OAL method, we use a smoothed non-parametric bootstrap approach that takes the model selection procedure into account. This procedure is described in efron2014estimation. The confidence intervals for all RF approaches are based on non-parametric bootstrapping. We use 500 bootstrap samples for each method. We then take the $0.025$ $(\alpha/2)$ and $0.975$ $(1-\alpha/2)$ quantile from the empirical distribution as the lower and upper bound for the confidence interval. We apply three non-parametric versions, the RF full (without regularization), the RRF (which uses regularization but no initial feature set), and our proposed OARF. We also apply the IPW method using a linear model to estimate the propensity score and the OAL method using the outcome-adaptive lasso to estimate the propensity score. We notice that the width of the confidence interval for the OARF is smaller than for the OAL method (in some settings only half as wide). The vanilla RF and the RRF do not achieve a coverage rate of 95% for any settings, while the OARF achieve the rate in 8 out of the 15 settings. The results show that some data generating processes might be too complex and hence are biased for the IPTW estimator. This is why increasing the sample size does not lead to a higher coverage rate. In all other settings, we see an increase in the coverage rate and a decrease in the width of the confidence intervals when increasing the sample size from N=500 to N=2000.

The OARF approach does not only decrease the variance when increasing the sample size but also the bias. Figure (ref) shows the mean squared error over 400 Monte Carlo replications for six different samples sizes (from 200 to 8000 observations). Here we use setting 4 and 5 as the data generating process since those settings have a low coverage rate of the confidence intervals. We compare the three random forest approaches, the full version, the regularized, and the OARF. The results show that the OARF achieves a significant decrease in MSE when increasing the sample size (e.g. for setting 5: from MSE = 0.35 for 200 observations, to MSE = 0.03 for 8000). The other two algorithms always have a higher MSE and the decrease is only slightly. The main reason for the high decrease in MSE when using the OARF is both, a decrease in variance and bias. To visualize this result we show boxplots in Figure (ref) and (ref) in the Appendix. They show the ATE estimation using the three random forest algorithms under varying sample sizes for each of the 400 Monte Carlo iterations. While the variance decreases in all algorithms when increasing the sample size, only the OARF has a noticeable decrease in bias.

In the Appendix, we discuss the possibility of tuning certain parameters in the random forest while keeping the objective to get balanced covariates rather than maximize the prediction accuracy. We also discuss a generalization of the OARF to other methods like the double machine learning, introduced by chernozhukov2018double. We find that using the OARF to estimate the propensity score decreases the bias compared to a full RF. In settings where even the OARF is biased, we additionally use the OARF to estimate the conditional mean of $Y$. This approach further decreases the bias.

Empirical Examples

In this section, we revisit two empirical examples. Both datasets are observational and contain a rich set of characteristics that may lead to selection into treatment. We use the same 5 methods as in the simulation study above. We report the ATE and a 95% confidence interval (CI). The datasets we use are freely accessible: RHC (\url{https://hbiostat.org/data/}), birth weight (\url{http://www.stata-press.com/data/r13/cattaneo2.dta}).

Re-analysis of SUPPORT data on Right Heart Catheterization

Right heart catheterization (RHC) is a diagnostic procedure used for critically ill patients. connors1996effectivenessused a propensity score matching approach to study the effectiveness of right heart catheterization in an observational setting. The authors found that after controlling for selection bias by using a rich set of covariates, RHC appeared to lead to higher mortality than not performing RHC. This conclusion is in contrast to popular belief among practitioners that RHC was beneficial.

The SUPPORT study collected data on hospitalized adult patients at five medical centers in the US. Based on information from a panel of experts, a rich set of characteristics, believed to be related to the decision on whether to perform the right heart catheterization or not, was collected. The study consists of data on 5735 individuals, RHC was performed on 2184 (the treatment group) while the remaining 3551 individuals did not get RHC (the control group). Treatment is equal to 1 if right heart catheterization was applied within 24 hours of admission, and 0 otherwise. The outcome of interest is an indicator for survival at 30 days. In total, we observe 72 covariates (including many dummy variables). After excluding variables with more than 50% missing values, we have 68 covariates left.

connors1996effectiveness matched treated and untreated patients based on the propensity score, with each unit, matched at most once. hirano2001estimation use different matching estimators (one of them is the exact matching) and their regression adjusted estimates range from $(-0.062, -0.053)$. crump2009dealing apply different samples based on the propensity score overlap and estimate an ATE from $-0.059$ to $-0.060$. ramsahai2011extending use propensity score and generic matching to investigate the effect of RHC on mortality within 180 days (note that the outcome variable is now an indicator for mortality, not survival). Their results are consistent with the previous findings, namely an estimated ATE using propensity score matching of $0.063$ and $0.046$ when using generic matching. A more recent study by li2018balancing considers different weighting strategies for covariate balancing. These strategies are based on the propensity score to estimate different target parameters. Their results also confirm previous estimates. For example, they estimate an ATE using overlap weighting of $-0.065$ and $-0.067$ when using optimal matching. All methods mentioned above use a logit model to estimate the propensity score. Given the rich set of covariates that consists of characteristics like age, race, income, and medical characteristics, there might be interaction effects and non-linear dependence. It would also be useful to know which covariates are true confounders and are selected when we do not assume any parametric form for both, the outcome model and the propensity score model.

As shown in Table (ref), all methods estimate a negative treatment effect, suggesting that performing RHC does decrease survival within the first 30 days. The results from the IPW and OAL are consistent with results from previous findings where the propensity score is estimated via a logit model. The RF results are smaller in magnitude with the lowest ATE of $-0.041$ by the OARF. The confidence intervals for the RF methods are tighter compared to the IPW or OAL method, while the latter method even includes zero. Among the 68 covariates, we find that the full RF selects the same variables as the RRF, in total 42 covariates with a proportion of $90-100\%$. The OAL method selects only 3 variables and 6 variables around 50% of the time (proportion between $50-60\%$). The OARF selects 9 variables in almost all iterations. All variables selected by OAL are also selected by OARF and additional variables are age, Duke Activity Status Index (DASI), APACHE score, Glasgow Coma Score (scoma1), white blood cell count Day 1 (wblc1), and bilirubin Day 1 (bili1). Variables that are not selected by OAL and OARF but with the RF and RRF are, among others, sex, race, education, and income. Especially characteristic variables that might be uncorrelated with a person's well being are not selected by the OARF but are still selected using the RRF using the same weights as the OARF. The complete list of covariates and their inclusion proportion is listed in Table (ref). Compared to OAL, the OARF includes age, DASI, white blood cells, and three other covariates.

Effect of smoking on birth weight

In this example, we reinvestigate the effect of maternal smoking status during pregnancy, the treatment variable, on babies’ birth weight, the outcome variable. We use a publicly available dataset that consists of 4642 singleton births in the USA. Additionally, we observe a rich set of characteristics like age, marital status, race, education, number of prenatal care visits, months since last birth, an indicator of firstborn infant, and indicator of alcohol consumption during pregnancy. All covariates are for the mother, except education, which we also observe for the father. The full dataset, containing more observations, was first used by almond2005costs who found a strong negative effect of maternal smoking during pregnancy on the weights of babies (about 200 -- 250 gram lighter for a baby with a mother smoking). In their study, the authors use a logit model to estimate the propensity score. We focus on estimating the propensity score without assuming any parametric form and again base the variable selection on features that are associated with the outcome. Table (ref) shows the ATE estimates along with 95% CI`s. We find similar negative effects as almond2005costs. The OARF estimates a decreased birth weight of -224 grams. Compared to the other RF models and the OAL method, the OARF has the tightest confidence intervals. Only the classic IPW estimator has a slightly tighter upper bound.

The OAL model and the full RF include 17 out of 19 variables while the RRF model includes 16 covariates. The OARF only includes 7 covariates and excludes, for example, an indicator for marital status, the education of mother and father, and if there were prenatal visits. The excluded variables are shown in Table (ref) and the complete list with inclusion probability is shown in Table (ref). The OARF is the only method that excludes the variable alcohol, meaning that it is not a confounding variable nor is it predictive on the outcome. The results, from a recent study by lundsberg2015low, suggest low-to-moderate alcohol exposure during early and late gestation is not associated with increased risk of low birth weight. Such findings are quite interesting since they allow a better understanding and interpretation of true confounding variables and of such that are not informative.

Discussion

We propose a non-parametric variable selection procedure for the estimation of treatment effects from observational studies. Building on outcome-adaptive penalization, we use a random forest to define variables that are predictive of the outcome. This allows the outcome model to deviate from a linear dependence on the covariates. We use a modified variable importance score to generate coefficients from the outcome model. We use the importance score to penalize variables that are spurious or that only predict the treatment but not the outcome. We show how to use penalty weights for each covariate to regularize the random forest that estimates the propensity score. Additional to the penalty, we show the importance of an initial feature space and how to include it in the RF. A Monte Carlo simulation shows that our proposed method, the OARF, has a smaller variance and produces unbiased estimates. In cases where all estimators are biased, the OARF produces the smallest bias and variance. The second goal of our proposed approach is to select only variables that have a relationship with the outcome (including all confounding variables) while excluding variables that predict the treatment and unimportant variables. Based on the simulation, we find that only the OARF selects the correct covariates and disregards all others. This holds for linear and non-linear settings and even if there is a strong correlation between the variables.

We apply the OARF and all other benchmark methods in two empirical examples. The ATE is comparable between all methods while the OARF shows tighter confidence intervals compared to the OAL and other RF methods. Regarding the variable selection, we find that the OARF selects and drops different variables compared to all other methods. This allows for a detailed evaluation of which variable might be responsible for selection bias and which variables are just spurious or only affect the propensity score.

\FloatBarrier

\phantomsection \addcontentsline{toc}{section}{References}