Variable Selection in Causal Inference using a Simultaneous Penalization Method

1 Introduction

In the analysis of observational data, when attempting to establish the magnitude of the causal effect of treatment (or exposure) in the presence of confounding, the practitioner is faced with certain modeling decisions that facilitate estimation. Should one take the parametric approach, at least one of two statistical models must be proposed; (i) the conditional mean model that models the expected outcome as a function of predictors, and (ii) the treatment allocation model that describes the mechanism via which treatment is allocated to (or, at least, received by) individuals in the study, again as a function of the predictors [1, 2].

Predictors that appear in both mechanisms (i) and (ii) are termed confounders, and their omission from models (i) and (ii) is typically regarded as a serious error, as it leads to inconsistent estimators of the treatment effect. Thus practitioners usually adopt a conservative approach, and attempt to ensure that they do not omit confounders by fitting a richly parameterized treatment allocation model. The conservative approach, however, can lead to predictors of treatment allocation only – and not outcome – being included in the treatment allocation model. The inclusion of such “spurious" variables in model (ii) is usually regarded as harmless. However, the typical forfeit for this conservatism is inflation of variance of the effect estimator [3, 4]. Moreover, [5] and [6] showed that inclusion of variables that are only predictors of treatment may cause bias [7, 8]. This problem also applies to the conditional mean model, but is in practice less problematic, as practitioners seem to be more concerned with bias removal, and therefore more likely to introduce the spurious variables in model (ii). Little formal guidance as to how the practitioner should act in this setting has been provided. In this paper, we refer to variables that only predict the treatment and are not associated with the outcome as treatment predictors.

As has been shown by [9], it is plausible that judicious variable selection may lead to appreciable efficiency gains, and several approaches with this aim have been proposed [10, 11]. However, confounder selection methods based on either just the treatment assignment model or just the outcome model may fail to account for non-ignorable confounders which barely predict the treatment or the outcome, respectively [12, 13]: in this manuscript, we use the term weak confounder for these variables. [14] shows that confounder selection procedures based on AIC and BIC can be sub-optimal and introduce a method based on the focused information criterion (FIC) which targets the treatment effect by minimizing a prediction mean square error (see also the cross-validation method of [15]). [16] introduces a Super Learner estimator which is computed by selecting a candidate from a set of estimators obtained from different models using a cross-validation risk [17, 18].

Bayesian adjustment for confounding (BAC) is a parametric variable selection approach introduced by [19]; BAC specifies a prior distribution for a set of possible models which includes a dependence parameter, w∈[1,∞], representing the odds of including a variable in the outcome model given that the same variable is in the propensity score model [20]. If we know a priori that a predictor of treatment is in fact a confounder, then w can be set to ∞ [13, 21]. [22] proposes a decision-theoretic approach to confounder selection that can handle high-dimensional cases; a Bayesian regression model is fit and using the posterior credible region of the regression parameters, a set of candidate models is formed. A sparse model is then found by penalizing models that do not include confounders. This method is conservative in the sense that it may include treatment predictors that may inflate the variance of the treatment effect [23]. Also, tuning the penalty function can be challenging. Doubly robust procedures have also been proposed where variables are selected using both the outcome and the treatment assignment models [12, 24].

Asymptotically, it is known that penalizing the conditional outcome model, given treatment and covariates, results in a valid variable selection strategy for causal effect estimation; however, for small to moderate sample sizes, it may result in the omission of weak confounders [25]. The objective of this manuscript is to improve the small sample performance of the outcome penalization strategy while maintaining its asymptotic performance (Table 2). We present a covariate selection procedure which facilitates the estimation of the treatment effect in the high-dimensional cases. Specifically, we propose a penalized objective function which considers both covariate-treatment and covariate-outcome associations and has the ability to select even weak confounders. This objective function is used to identify the set of non-ignorable confounders and predictors of outcome; the resulting parameter estimates do not have any causal interpretation. We derive the asymptotic properties of procedure and show that under some mild conditions the estimators have oracle properties (specifically, are consistent and asymptotically normally distributed). We utilize the selected covariates to estimate the causal effect of interest using a doubly robust estimator.

2 Preliminaries & notation

In standard notation, let Y(d) denote the (potential) outcome arising from treatment d, and let D denote the treatment received. We consider for illustration the case of binary treatment. The observed outcome, Y, is defined as DY(1)+(1−D)Y(0). We restrict attention here to the situation where each predictor can be classified into one of three types, and to single time-point studies. We consider

1. treatment predictors (X1), which are related to treatment and not to outcome.

2. confounders (X2), which are related to both outcome and treatment.

3. outcome predictors (X3), which are related to outcome and not to treatment;

see the directed acyclic graph (DAG) in Figure 1.

Figure 1

Covariate types: Type-I: X1, Type-II: X2 and Type-III: X3.

In addition, as is usual, we will make the assumption of no unmeasured confounders, that is, that treatment received D and potential outcome to treatment d, Y(d), are independent, given the measured confounders, i.e., X2. In any practical situation, to facilitate causal inference, the analyst must make an assessment as to the structural nature of the relationships between the variables encoded by the DAG in Figure 1.

3 Penalization and treatment effect estimation

In this section, we present our proposed method for estimating the treatment effect in high-dimensional cases. We separate the covariate selection and treatment effect estimation procedure. First, we form a penalized objective function which is used to identify the important covariates, and establish the theoretical properties of the resulting estimators (i.e., minimizers of the penalized objective function). Note that because of the special characteristics of this function the estimators do not have any causal interpretation and are used just to prioritize variables. Second, treatment effect estimation is performed using a doubly robust estimator with the selected covariates. We use a simple model structure to illustrate the methodology, and assume that a random sample of size n of observations of outcome, exposure and covariates is available.

In order to present the method, we initially assume that columns of X are orthogonal; this simplifying assumption is relaxed in Appendix C of the Supplementary Materials.

3.1 Penalized objective function

Consider the following linear outcome model under a binary exposure d

Y=θd+xαy+ϵ,

where x is a 1×r standardized covariate vector, and ϵ is standard normal residual error. Assuming a logit model for the propensity score, we have

π(x,αd)=p(D=1|x,αd)=exp{xαd}1+exp{xαd}.

Let X denote the corresponding n×r design matrix. In the formulation, αy and αd denote r×1 vectors of parameters and θ is the treatment effect.

First, we form an objective function, M(α), in which the coefficients αy and αd are replaced by a single common vector of coefficients α that is proportional to a weighted sum of |αy| and |αd|, where |.| denotes the componentwise absolute value. This will guarantee that the objective function satisfies the following condition:

Argmin Condition: An element, αˆj, of the minimizer αˆ of M(α) converges to zero as n→∞ if the corresponding covariate is not associated with both treatment and outcome.

In standard linear regression where we regress Y on a vector of covariates X (with no treatment variable) that are presumed orthogonal, an example of an objective function that estimates |αy| is

M(α)=12n∑i=1nxi⊤yi−nα⊤∑i=1nxi⊤yi−nα

yielding the estimating equation and estimate

∂M(α)∂α=∑i=1nxi⊤yi−nα=0⟹αˆ=1n∑i=1nxi⊤yi.

We have used the fact n−1x⊤x is an identity matrix due to standardization. It is clear that as n⟶∞, αˆ⟶p|E[xi⊤Yi]| Similarly, in a binary exposure setting, |αd| can be estimated using the objective function −∑i=1nxidiα+∑i=1nlog(1+exp{xiα}) which is based on minus the log likelihood for a logistic regression. Now, we combine these two objective functions and show that the resulting objective function have some interesting features. First, we obtain the least squares (or ridge) θ˜ of θ by regressing Y on D and the vector of covariates \textbf{X}. Because of the inclusion of spurious and treatment predictors in the model, θ˜ is not efficient but is but is consistent for θ under the no unmeasured confounders assumption and a correctly specified outcome model. Define Y˜=Y−θ˜D. Let

(2)M(α)=12[|∑i=1n\textbf{x}iy˜i|−α⊤∑i=1n(\textbf{x}i⊤\textbf{x}i)]∑i=1n\textbf{x}i⊤\textbf{x}i−1[|∑i=1n\textbf{x}iy˜i|−α⊤∑i=1n(\textbf{x}i⊤\textbf{x}i)]⊤+1τ[−|∑i=1nxidi|α+∑i=1nlog(1+exp{xiα})],

where τ is a positive constant. Under orthogonality of X and standardization, we have

M(α)=12n[|∑i=1nxiy˜i|−nα⊤][|∑i=1nxiy˜i|−nα⊤]⊤+1τ[−|∑i=1nxidi|α+∑i=1nlog(1+exp{xiα})].

Note that for each j, the parameter αj corresponding to xj is the same in both models. Now, we show that in this function the absolute values play a critical rule in satisfying the Argmin condition. Let αˆ=argminαM(α) and α˜y be the least squares estimate of the parameters in the outcome model. By convexity of M(α), αˆ must be the unique solution to ∂M(α)/∂α=0, which implies that

(3)nαˆ+1τ∑i=1n\textbf{x}i⊤exp{\textbf{x}iαˆ}1+exp{\textbf{x}iαˆ}=∑i=1n\textbf{x}i⊤y˜i+1τ∑i=1n\textbf{x}i⊤di=nα˜y+1τ∑i=1n\textbf{x}i⊤di.

This equation shows that αˆ=0 (i.e. the Argmin condition is satisfied) if cov(X,Y)=cov(X,D)=0. Note that cov(X,D)=0 implies cov(X,Y˜)=cov(X,Y). Let α˜d be the maximum likelihood estimate of the parameters in the treatment model. Then, using the first two terms of a Taylor expansion, we have ∑i=1nxi⊤di≈nα˜d/4, exp{xiαˆ}/(1+exp{xiαˆ})≈1/2+αˆxi/4, and thus

αˆ≈4τ4τ+1α˜y+14τ+1α˜d.

Hence, αˆ is a weighted sum of |α˜y| and |α˜d|. The constant τ controls the contribution of components α˜y and α˜d to the estimator. For example, for τ=2, αˆ≈811α˜y+111α˜d; while for τ=0.1, αˆ≈0.41.4α˜y+11.4α˜d. Thus, τ gives a flexibility to our objective function such that as it decreases to zero the proposed estimate αˆ converges to α˜d. See Section 3.2 for more detailed discussion.

Figure 2 visually presents how αˆj behaves, for a fixed τ=0.5, when αˆjy converges to zero as sample size increases, i.e., the jth variable becomes insignificant. Specifically, this figure presents a case where there is just one covariate (i.e., j=1) and the coefficient of this covariate in outcome and treatment models are αjy=1/n and αjd=0.3, respectively, where n is the sample size. As expected, αˆ corresponding to this covariate does not converge to zero as sample size increases. The same behavior would be observed if αjd→0 as n→∞ and αjy is a non-zero constant.

Figure 2

Performance of the modified objective function based estimator for different sample sizes n. The dashed and solid lines are 1/n and the estimated coefficient αˆ using the modified objective function.

Therefore, penalizing objective function (2) results in selecting covariates that are either related to outcome (Type-III) or treatment (Type-I). However, this is against our goal of keeping variables that are either predictors of the outcome or non-ignorable confounders and excluding treatment predictors (Type-I). To deal with this problem, we present a weighted lasso penalty function that is tailored specifically for causal inference variable selection:

λ∑j=1r|αj|(α˜jy)2(1+|α˜jd|)2,

where α˜y and α˜d are the least squares (or ridge) and maximum likelihood estimates of the parameters in the outcome and treatment models, respectively, and λ is a tuning parameter. Thus, the proposed modified penalized objective function is given by

(4)Mp(α)=M(α)+λ∑j=1r|αj|(α˜jy)2(1+|α˜jd|)2.

We refer to αˆ=argminαMp(α) as penalized modified objective function estimators (PMOE). The magnitude of the penalty on each parameter is proportional to its contribution to the outcome and treatment model. Note that as αˆjy→0, our penalty function puts more penalty on the jth parameter while considering the covariate-treatment association. For example, when a covariate barely predicts the outcome and treatment, our proposed penalty function imposes a stronger penalty on the parameter compared to a case where a covariate barely predicts the outcome and is strongly related to treatment. This is an important feature of the proposed penalty function which allows selecting such weak confounders for small sample sizes. The proposed penalized estimator asymptotically selects the same covariates as the adaptive lasso on the outcome model [25, 32]. Thus, the proposed method improves the small sample performance of the outcome model penalization strategy while maintaining the asymptotic properties of this strategy.

4 Simulation studies

In this section, we study the performance of our proposed variable selection method using simulated data. We compare our results with BAC method introduced by [19], the Bayesian credible region (Cred. Reg.) introduced by [22], the collaborative targeted maximum likelihood estimators (CTMLE) introduced by [24], a doubly robust approach (Belloni) proposed by [12] and outcome penalized estimator (Y-fit) where we regress the outcome against the treatment and covariates, but only penalize the covariates. Our simulation includes a scenario in which there is a weak confounder that is strongly related to the treatment but weakly to the outcome. We consider linear working models for g() throughout this section.

We generate 500 data sets of sizes 300 and 500 from the following two models:

1. [1.] Y∼Normal(D+2X1+0.5X2+5X3+5X4,4)

2. [2.] Y∼Normal(D+2X1+0.2X2+5X3+5X4,4)

where in both models

and Xk has a N(1,4) for k=1,...,100. Note that in the second scenario, X2 is considered as a weak confounder. Results are summarized in Table 1; the Y-fit row refers to the estimator obtained by penalizing the outcome model using LASSO penalty, and the Oracle row refers to estimator obtained by including (X1,X2,X3,X4) in the propensity score and the outcome model (5). We studied the performance of PMOE for different values of constant τ=0.1,0.5,1,2,4,100 and the optimal one that is chosen by the GCV criterion, i.e., PMOEτ=opt. In the first scenario there is no weak confounder and the Y-fit is omitted since the results are similar to the PMOEτ=20 row. The Bayesian adjustment for confounding (BAC) has been implemented using the R package BEAU with ω=∞ and the Bayesian credible region (Cred. Reg.) has been implemented using the R package BayesPen with flat prior. The doubly robust variable selection method, i.e., Belloni, is also implemented using the R package hdm. Finally, we implemented CTMLE.

Figure 3

Proportion of times that the algorithm selects certain values of τ. The left panel presents the first simulation scenario where there is no weak confounders while the right panel presents the second simulation scenario where there is a weak non-ignorable confounder.

Figure 3 shows the proportion of times that our algorithm selects a certain value of tuning parameter τ. In scenario 1, where there is no weak confounder, the algorithm tends to select higher values for τ which means the procedure adaptively selects the important covariates mostly based on the outcome model; while in scenario 2, where X2 is a non-ignorable weak confounder, the algorithm puts more weights on the treatment model by selecting smaller values for τ.

The variance of the estimator in the BAC is too large due to the inclusion of spurious variables that are not related to the outcome. The PMOE and Cred. Reg. estimators, however, are unbiased and have smaller variance. In fact, PMOE has the lowest variance compared to BAC and Cred. Reg. methods regardless of the value of τ. In the second scenario, the Y-fit estimator is biased because of under selecting the confounder X2. Also, the bias of the proposed PMOE estimator increases by increasing τ that is expected because as τ increases the proposed method should perform similarly to the outcome based variable selection methods such as Y-fit. In general, doubly robust procedures, e.g., CTMLE and Belloni, are sensitive to the positivity assumption. Specifically, when the estimated propensity score using the selected variables is close to boundaries, such methods lead to unstable estimates and possibly under covered confidence intervals. Both CTMLE and Belloni produce confidence intervals that are under covered, and CTMLE seems to be affected more seriously by the issue discussed above. The main drawback of Belloni’s method is that it fails to eliminate variables that only predict the treatment from the set of selected variables. As our results in Table 1 shows this deficiency leads to inflating the standard error of the corresponding estimates. The estimator PMOEopt outperforms all the other estimators and has similar performance as the oracle estimator which confirms that our method successfully tunes τ.

Table 2 presents the average number of coefficients set to zero correctly and incorrectly under the second scenario. There are four non-zero coefficients in our generative model so the number in the correct column should be 96 and in the incorrect column should be 0. This table shows that Cred. Reg., and BAC are somewhat conservative and include some of the variables that should not be included which can be the source of the observed variance inflation in Table 1. This is mostly due to inclusion of variables that are predictors of treatment model but have no association with the outcome. Also, increasing τ in PMOE, increases the chance of setting the coefficient of X2 to zero. The Y-fit row shows that this method is setting a nonzero coefficient to zero (i.e., coefficient of X2) which explains the bias in Table 1. The tuned estimator PMOEopt has satisfactory performance in ignoring covariates that are only predictors of the treatment mechanism and selecting all the non-ignorable confounders. This, in fact, highlights the importance of our proposed method. We couldn’t obtain the information needed in Table 2 using the available CTMLE’s and Belloni’s R codes, and thus, they are omitted.

Simulation studies presented in Appendix D of the Supplementary Materials study the performance of our covariate selection and estimation procedure when the covariates are non-orthogonal. Web Tables 2 and 3 presented in Appendix D show that the proposed method is still outperforming the other methods. Moreover, in Appendix B of the Supplementary Materials, we study cases that the number of covariates is larger than the sample size (r>n). We also investigate cases where either of the working models of the propensity score or the outcome model is misspecified. Our results show that the proposed method performs well, and outperforms Y-fit (Web Table 1).

5 Application: the effect of life expectancy on economic growth

In this section, we use the proposed method to analyse the cross-country economic growth data reported by [47] and compare our findings with similar studies on the same data set. For illustration purposes, we focus on a subset of the data which includes 88 countries and 35 variables. Additional details are provided in [48]. We are interested in selecting non-ignorable variables which confound the effect of life expectancy (exposure variable) as a measure of population health on the average growth rate of gross domestic product per capita in 1960-1996 (outcome).

The causal (or, at least, unconfounded) effect of life expectancy on economic growth is controversial. [49] find no evidence of increasing life expectancy on economic growth while [50] shows that it might have a positive effect. We dichotomize the life expectancy based on the observed median, which is 50 years. Hence, the exposure variable D=1 if life expectancy is below 50 years in that country and 0 otherwise.

We applied our method described in Section 3 to selection important confounders. In our analysis, Y-fit refers to the case where just the outcome model is penalized using LASSO to select the significant covariates. We also implement the BAC with ω=∞ and Credible region with flat prior methods. Belloni’s and CTMLE methods are also implemented. Our procedure selects τ=0.1 as the optimal value for this tuning parameter. The estimates obtained by PMOEτ=0.1, i.e., 0.617 with CI: (−1.27,1.361), and CTMLE, i.e., 0.719 with CI: (0.512,0.927) are fairly close. However the standard error of the latter one is 70% smaller that the former, which raises concern about the validity of the corresponding confidence intervals (see Table 1). Y-fit and Belloni lead to the smallest and largest effects, respectively. Based on our simulation studies, we conjecture that this is because Y-fit may ignore some of non-ignorable confounders and Belloni’s may include some of the treatment predictors in the model. The BAC estimate, i.e., 0.524 with CI:(−0.228,1.286), is also close to the PMOEτ=0.1. The Credible region method leads to the second highest effect estimate of 0.820 with CI: (−0.203,1.341).

Table 3 presents the list of variables that are selected at least by one of the methods. For brevity, we focus on Y-fit, BAC, Cred. Reg., and PMOE. The proposed method selects 12 and 9 variables depending on the value of τ while Y-fit, BAC and Cred. Reg. select 7, 6 and 11 variables, respectively. Y-fit and PMOEτ=5 are susceptible to under selecting some of non-ignorable confounders that barely predict the outcome. Specifically, our results suggest that Population Density 1960, and Initial Income are such non-ignorable confounders which are known to be important confounders in the economics literature. Table 4 shows that ignoring these variables leads to a substantially different treatment effect estimates. Moreover, because there are evidence that variables such as Oil producing Country, Land Area Near Navigable Water, and Nominal Government Share may confound the effect of the life expectancy on the economic growth [48, 49, 51, 52, 53, 54, 55, 56], PMOEτ=0.1 seems to select more reasonable covariates than PMOEτ=0.5. Despite the small discrepancies in the selected variables between the BAC, Creg. Reg. and PMOEτ=opt, these methods lead to the same conclusion and suggest that while the life expectancy has positive effect on the average growth rate of gross domestic product per capita in 1960-1996, the effect is not statistically significant.

To gain insight into the effect of parameter τ on the selection results, we plot the estimated coefficients α for different values of τ for given tuning parameters. In Figure 4, the blue solid (dark dashed) lines correspond to coefficients that their estimated value is (not) greater than 0.05 for the entire range of τ. Also, displays (a)–(d) correspond to tuning parameter values λ=0.0, 0.0001, 0.001, and 0.01, respectively. Figure 4 shows that, for different values of τ, the selected set of covariates may vary slightly. For example, in Figure 4(c) where a moderate penalty function is imposed, the proposed method suggests to include Colony (the dashed line) to the selected set for larger values of τ. Also, as the penalty becomes stronger, the estimated coefficients become more stable across values of τ (Figure 4(d)). This is because, for larger values of λ, the penalty function has a more dominant rule on the variable selection than τ.

Figure 4

Economic growth data: The plot of the estimated parameters using the proposed method given different values of parameter τ. The tuning parameter λ is fixed at 0.0, 0.0001, 0.001, and 0.01 in plots (a), (b), (c) and (d), respectively. The blue solid (dark dashed) lines correspond to coefficients that their estimated value is (not) greater than 0.05 for the entire range of τ.

Table 4 also reports the standard errors of the estimated effect of life expectancy using different penalization methods. The standard error of the PMOE estimator is approximated using an idea similar to [46]. Specifically, we bootstrap the sample and in each bootstrap force the components of the penalized estimator αˆ to zero whenever they are close to zero and estimate the treatment effect using the selected covariates, i.e., we define αˆ†=αˆ1(|αˆ|>1/n). We utilize this thresholded bootstrap method to approximate the standard error of the treatment effect. Although more investigation is required to validate the asymptotic properties of this method, we only use this standard errors to shed light on the behavior of the penalized estimators. For example, the estimators correspond to PMOEτ=5 and Y-fit have the lowest standard errors that may support the possibility of under-selecting important covariates. Our results suggests that although the effect of life expectancy is positive, it is unlikely to be significant that is consistent with [49].

6 Discussion

We have established a two-step procedure for estimating an unconfounded treatment effect in high-dimensional settings. First, we deal with the sparsity by penalizing a modified objective function which considers both covariate-outcome and covariate-treatment associations. Then, the selected variables are used to form a doubly robust regression estimator of the treatment effect by incorporating the propensity score in the conditional expectation of the outcome. The selected covariates may be used in other causal techniques as well as the proposed regression method. The proposed method may also be used to identify candidate instrumental variables [57, 58,59, 60]. Specifically, one can penalize the treatment model first and record the selected variables and then apply the proposed method and identify the candidate instrumental variables as those that are selected by the treatment model but not the proposed method [23].

As described in section 3.5, any covariate selection procedure which involves the outcome variable affects the subsequent inference of the selected coefficients [61, 62, 63, 64, 65]. This is because the selected model itself is stochastic and it needs to be accounted for. [66] proposes a method to produce a valid confidence interval for the coefficients of the selected model in the post-selection context. In our setting, although we do not penalize the treatment effect, the randomness of the selected model affects the inference about the causal effect parameter through confounding. Moreover, note that the oracle property of the penalized regression estimators is a pointwise asymptotic feature and does not necessarily hold for all the points in the parameter space [67, 68]. In this manuscript, we assume that the parameter dimension (r) is fixed while the number of observation tends to infinity. One important extension to our work is to generalize the framework to cases where the tuple (n,r) tends to infinity [69]. Analyzing the convergence of the estimated vector of parameters in the more general setting requires an adaptation of restricted eigenvalue condition [70] or restricted isometry property [71].