Valid causal inference with unobserved confounding in high-dimensional settings

1 Introduction

In contrast to randomized experiments, observational studies are prone to the presence of confounding variables that are not balanced among treated and control individuals [1]. In such studies, all confounders are often assumed to be observed for identifiability of the causal parameter of interest [2,3]. Efforts to make this assumption plausible and the use of flexible nuisance models often yield a high-dimensional setting, where the number of nuisance parameters to be fitted may be at least of the order of the sample size. Nonetheless, some confounders may still be unobserved. Therefore, well-conducted observational studies should investigate the sensitivity of the inference to the assumption of no unmeasured confounding [4].

In this article, we study sensitivity analysis for semiparametric inference on the average causal effect of a binary treatment in high-dimensional observational studies. In this context, under certain conditions, augmented inverse probability weighting (AIPW) [5] and targeted learning [6] estimators yield uniformly valid inferences even when nuisance models are fitted with flexible machine learning algorithms (see, e.g., Farrell [7], Van der Laan and Gruber [8], Chernozhukov et al. [9] and Belloni et al. [10], Moosavi et al. [11] for a recent review). Roughly, uniformly valid inference means that the finite sample behavior of the estimator can be well approximated by its asymptotic distribution, even if preliminary analysis, such as variable selection, has been performed on the nuisance models prior to final estimation. Farrell [7] studied uniformly valid inference associated with the AIPW estimator when flexible estimators of nuisance models are weakly consistent and fulfill multiplicative rate conditions (Assumption 2). Here, we build on these results to allow for unobserved confounding and study the uniform validity of the resulting inference. For this, we first specify the confounding bias of the AIPW estimator of the causal effect as a function of unobserved confounding. We then propose uncertainty intervals for the causal effect that account for such bias [12,13]. Our suggested inference on the causal effect using the uncertainty intervals ignores the finite sample bias and randomness in the estimation of the confounding bias. It provides uniformly valid inference given assumptions on the amount of unobserved confounding relative to the sample size; the latter is formalized in terms of convergence rate. In particular, unobserved confounding cannot be arbitrarily large. While this may seem like a restriction, our simulation experiments show that valid inference is obtained with relevant amounts of unobserved confounding.

This work contributes to the sensitivity analysis literature that originated in Cornfield et al. [4] and has been studied further by many others (to cite only a few: [14–20]). More specifically, we contribute to research that considers the sensitivity parameter to be the correlation between the potential outcomes and the treatment assignment given the observed covariates induced by unmeasured confounders [13,21–23]. Our contribution is to consider post-model-selection inference, including high-dimensional situations, and inference following machine learning fits of the nuisance models. In this respect, our work is related to Nabi et al. [24], which also considers inference on a causal parameter based on flexible estimation of nuisance models and allows for unobserved confounders via sensitivity analysis. Their approach differs in how the confounding bias is parameterized. We contrast the methods using a case study in Section 4.

The rest of this article is organized as follows. Section 2 describes the context and states our result on the uniform validity of the inference for a post-model-selection estimator of a causal parameter, for a given amount of unobserved confounding described by a correlation parameter. In practice, this correlation is unknown and a plausible range of correlation values may be considered as a sensitivity analysis. In Section 3, simulations are provided to illustrate the relevance of the theory for finite samples, in both low- and high-dimensional settings. Section 4 presents a case study of the effect of maternal smoking on birth weight. This illustrates the use of the herein proposed methods implemented in the R-package hdim.ui (https://github.com/stat4reg/hdim.ui). Section 5 concludes this article. All proofs are given in the Appendix.

2 Theory and method

We aim to study the causal effect of a binary treatment T on a continuous outcome of interest Y . Let Y ( t ) , t ∈ { 0 , 1 } , be the potential outcomes that would have been observed under a corresponding treatment level, and assume the observed outcome is Y = T Y ( 1 ) + ( 1 − T ) Y ( 0 ) , for all units in the study. The average causal effect of the treatment is defined as E ( Y ( 1 ) − Y ( 0 ) ) . Without loss of generality, we focus on the parameter τ = E ( Y ( 1 ) ) ; (see the Appendix for other parameters of interest). We denote the vector of pre-treatment covariates and their transformations by X = ( 1 , X ( 1 ) , … , X ( p − 1 ) ) . For simplicity, we call this the set of covariates. The dimension p of this vector is considered to increase with the sample size n .

Consider O i , n = { ( X i , n , Y i , n , T i , n ) } i = 1 n to be an i.i.d. sample that follows a distribution P n , which is allowed to vary with the sample size n . We drop the subscript n , which implies dependence on n , when it is clear from the context. Moreover, let n t denote the number of individuals with observed treatment T = 1 . We use the notation E n [ W i ] = n − 1 Σ i = 1 n W i , E n [ W i ] q = ( n − 1 Σ i = 1 n W i ) q , and E n t [ W i ] = n t − 1 Σ i = 1 n T i W i .

We consider the augmented inverse propensity weighting estimator (AIPW; see Robins et al. [5] and Scharfstein et al. [25]):

where m ˆ ( X ) and e ˆ ( X ) are the arbitrary estimators of m ( X ) = E ( Y ( 1 ) ∣ X , T = 1 ) and e ( X ) = E ( T ∣ X ) , respectively. Let us consider the following assumptions.

The following theorem gives the asymptotic behaviour of the AIPW estimator and its asymptotic target parameter. The theorem generalizes Theorem 3 in Farrell [26] by avoiding the assumption E ( Y ( 1 ) ∣ X = x , T = 1 ) = E ( Y ( 1 ) ∣ X = x ) , i.e., allowing for unobserved confounders.

Assumption 1 can be checked empirically. Assumption 2 is central. It allows data-adaptive fits of the propensity score and outcome model with rates lower than n 1 ⁄ 2 . For instance, Assumption 2b is fulfilled if both nuisance functions are fitted at the n 1 ⁄ 4 rate. Such rates are attainable by machine learning algorithms [27], e.g., Lasso under sparsity conditions [e.g., 7,11], Deep neural networks under smoothness conditions [26,28], and other nonparametric smoothers and combinations thereof via super-learners [8]. Thus, Theorem 1 ensures a parametric rate of convergence for AIPW, e.g., in sparse high-dimensional situations.

Following Genbäck and de Luna [13], we study the consequences of the presence of unobserved confounders of the treatment–outcome relationship on inference about τ by postulating a sensitivity model:

It is generally agreed upon that a sensitivity model should make as few restrictions as possible on the observed data distribution [29,30]. Assumption 3 makes one such slight restriction, i.e., assuming the normality of η (probit model for the propensity score). This restriction is, however, rather weak since it is typically the case that changing the link function, e.g., from probit to logit (corresponding to the logistic distribution for η ) does not substantially change the fitted probabilities. We study the effect of misspecifying the distribution of η in our simulation experiments.

The case ρ = 0 corresponds to the inclusion of all confounders in X . Assumption 3, therefore, considers a departure from the commonly used identification assumption E ( Y ( 1 ) ∣ X = x , T = 1 ) = E ( Y ( 1 ) ∣ X = x ) (no unobserved confounders) in order to perform a sensitivity analysis. Assumption 3 allows us to consider a potential amount of unobserved confounding through an interpretable parameter ρ , as well as to compute the size of the confounding bias of the AIPW estimator.

Note that under the conditions of Theorem 1, the AIPW estimator (1) is a semiparametric efficient estimator of τ − [see, 31, Example 5]. However, if confounders of the treatment–outcome relationship are not observed (not included in X ), ρ and consequently, b is nonzero, and hence, the AIPW estimate is not a consistent estimate of τ . It can be shown that for given values of ρ and σ , the efficient influence function of the bias term has the form ρ σ E ( λ ( g ( X ) ) ) − b . This motivates us to estimate the confounding bias b empirically by plugging in the estimated g ( X i ) values. For a more realistic scenario, the parameter σ must also be estimated. An intuitive choice is to use σ ˆ = E n t [ ( Y ( 1 ) i − m ˆ ( X i ) ) 2 ] 1 ⁄ 2 . This estimator of the variance, σ ˆ 2 , is biased due to both high dimensionality/variable selection and unobserved confounders. Theorem A1 in the Appendix provides the asymptotic properties of τ ˆ AIPW − b ˆ as an estimator of τ , assuming approximate sparsity and known ρ .

Theorem 3 introduces a corrected estimator of σ . Given the true value of the sensitivity parameter, this corrected estimator of the variance gives us a consistent estimate of the causal parameter τ .

Theorem 4 provides further asymptotic properties of τ ˆ AIPW − b ˆ c as an estimator of the average causal effect τ .

As a corollary, the following ( 1 − α ) % confidence interval for τ given ρ is uniformly valid.

In practice, ρ is not known and a sensitivity analysis is obtained by considering the 95% uncertainty interval constructed as the union of all 95% confidence intervals obtained by varying ρ ∈ ( ρ min , ρ max ) , a plausible interval for ρ . If the latter interval contains the true value of ρ , then this uncertainty interval covers the parameter τ with at least 95% probability [32, Corollary 4.1.1].

When no parametric models are assumed and/or the number of covariates is greater than the number of observations, lasso regression [33] can be used in the first step to select low-dimensional sets of variables for fitting linear and probit nuisance models. In other words, a linear (in the parameters) regression can be fitted using variables selected by a preliminary lasso regression. A probit regression for the treatment can be fitted using variables selected by a preliminary probit-lasso regression. When it comes to the rate conditions in Assumption 2, one can use a post-selection linear model fit for the outcome under common sparsity assumptions on the true data-generating process (for post-lasso regression, see [7], Corollary 5). We are unaware of any similar result for an estimator in a sparse probit model. However, we investigate the performance of a post-lasso probit regression in the simulation section.

In the following sections, we use the R package hdim.ui (https://github.com/stat4reg/hdim.ui), where the estimators mentioned previously have been implemented and can be used to obtain uncertainty intervals and perform sensitivity analysis in high-dimensional situations. The package builds on code from the ui package (https://cran.r-project.org/web/packages/ui/index.html) [13].

3 Simulation study

The simulation study in this section aims to illustrate the finite sample properties of the asymptotic approximations obtained in the previous section. Here, we consider a high-dimensional setting, while the results for low-dimensional settings are reported in the Supplementary Material. Data were generated according to the model in Assumption 3. All the covariates are generated to be independently and normally distributed with mean 0 and variance 1. The error terms are generated using ( η , ξ ) ∼ N ( 0 , Σ ) with Σ = 1 ρ ρ 1 . We let p = n and simulate outcome and treatment as

where both models include covariates that are weakly associated with the dependent variable. Such covariates are likely to be missed in a variable selection step. To study the method behavior in a case where the probit link for the treatment model is misspecified, we also simulate the error η from a logistic distribution with mean 0 and variance 1, while keeping the correlation with ξ equal to ρ .

In order to investigate the performance of the inference under variable selection, we use τ ˆ AIPW refit − b ˆ refit (using σ ˆ , Theorem A1) and τ ˆ AIPW refit − b ˆ c refit (using σ ˆ c , Theorem 4). In these estimators, m ˆ ( X ) is found by refitting the linear regression using variables selected by preliminary linear-lasso regression and g ˆ ( X ) is found by refitting the probit model using variables selected by preliminary probit-lasso regression. The linear-lasso regression is implemented using the hdm package in R [34,35]. The probit-lasso is implemented using the glmnet package [36], where the regularization parameter is found by cross-validation. Then, 95% confidence intervals are constructed for a range of values of the sensitivity parameter ρ according to Corollary 1 using the influence curve-based standard error estimator. If empirical coverages are close to nominal, then the corresponding uncertainty intervals will be conservative by construction, as mentioned previously [see 13,32]. We are particularly interested in investigating coverages of confidence intervals when varying ρ as sample size increases ( n = 500 , 1,000, 1,500). This is because Theorem A1 assumes n ρ 3 = o ( 1 ) when using the biased estimator of outcome error variance, while the weaker condition ρ = o ( 1 ) is used in Theorem 4 when the corrected variance estimator is used. The results of the simulation study are based on 500 Monte-Carlo replications.

Table 1 reports, respectively, biases and empirical coverages of the 95% confidence intervals for τ using τ ˆ AIPW refit − b * , τ ˆ AIPW refit − b ˆ refit and τ ˆ AIPW refit − b ˆ c refit . Here, b * is an approximation of the confounding bias b using the true values of σ and ρ and the Monte Carlo estimate of E ( λ ( g ( X ) ) ) .

Biases are as expected negligible when correcting with the true bias b * . Also as expected, when the bias b is estimated, biases in the estimation of τ are largest for the larger values of ρ . The confidence intervals constructed using τ ˆ AIPW refit − b * have minimum empirical coverage of 0.91. When a plug-in estimator of the confounding bias is used low empirical coverage may arise when ρ is too large (relative to the sample size). However, because of the weaker condition on ρ in Theorem 4, the coverage for the estimator τ ˆ AIPW refit − b ˆ c refit is, as expected, closer to nominal level compared to the estimator τ ˆ AIPW refit − b ˆ refit , and this improvement is more pronounced for larger sample sizes.

4 Case study: Effect of maternal smoking on birth weight

We re-examine a study that aims to assess the effect of smoking during pregnancy on birth weight [24,37,38]. The data come from an open-access sample of 4996 individuals, a sub-sample of approximately 500000 singleton births in Pennsylvania between 1989 and 1991 (see Nabi et al. [24] and resources therein for more details on the data). Weight at birth in grams is the outcome variable. The treated group consists of pregnant women who smoked during pregnancy, and the control group consists of women who did not smoke.

We use the same set of covariates as in Nabi et al. [24]. As they argue, sensitivity analysis is necessary to account for potential unobserved confounders, such as genetic factors not observed. The observed covariates that we consider are maternal data including numerical (age and the number of prenatal visits) and categorical variables (education -less than high school, high school, more than high school- and birth order -one, two, larger than two) and binary variables including white, hispanic, married, foreign and alcohol use. Higher-order and interaction terms of the numerical variables of order up to three are also considered. Finally, all the second-order interactions with categorical variables are added, which gives a total of 80 terms.

The analysis in Almond et al. [37] estimates that smoking has an average effect of − 203.2  g on birth weight, determined by a regression model without taking into account the effect of unobserved confounding. In Nabi et al. [24], a semiparametric approach yields an estimate of − 223  g (95% confidence interval [ − 274 , − 172 ]). However, after accounting for unobserved confounding in a sensitivity analysis, the latter study suggests that the average effect is not more extreme than − 200  g (see (2) for details on the clinical assumptions used).

In our analysis, we use the AIPW estimator of average causal effect, E ( Y ( 1 ) ) − E ( Y ( 0 ) ) . The results presented in this article are directly applicable to E ( Y ( 0 ) ) as well, and, therefore, to E ( Y ( 1 ) ) − E ( Y ( 0 ) ) under corresponding assumptions. In particular, there are now two sensitivity parameters through the sensitivity model of Assumption 3, for both Y ( 1 ) and Y ( 0 ) , denoted, respectively, ρ 1 and ρ 0 . We use subset = ‘refit’ in function ui.causal to exploit the built-in variable selection function (ui.causal included in the R-package hdim.ui). The resulting AIPW estimation gives an estimate of the average causal effect of − 221  g (95% confidence interval: [ − 273 , − 168 ]).

In their sensitivity analysis, Nabi et al. [24] argued that the following inequalities make clinical sense:

For instance from the second line, under smoking, non-smokers are expected to have higher average birth weight than smokers, E ( Y ( 1 ) ∣ T = 0 ) > E ( Y ( 1 ) ∣ T = 1 ) . We choose to also use these prior clinical assumptions, which yield bounds for our sensitivity parameters ρ 1 and ρ 0 . That is because the terms on the right- and left-hand sides of the inequalities in (2) can be unbiasedly estimated using sample averages since the data are fully observed, whereas the terms in the middle are functions of the sensitivity parameter. Appendix provides a description of how they are estimated. According to the range of sensitivity parameters obtained from the above restrictions (Figure 1), we have derived estimates and the uniformly valid confidence intervals along with a resulting uncertainty interval for E ( Y ( 0 ) ) and E ( Y ( 1 ) ) (Figure 2). The resulting uncertainty interval for the average causal effect of maternal smoking on birth weight is [ − 333 , 49]. Thus, while our AIPW estimation and 95% confidence interval under the assumption that all confounders are observed are similar to those of Nabi et al. [24], the sensitivity analysis is quite different, since they suggest that the average effect is not more extreme than − 200  g while our analysis does not discard a potentially much larger adverse effect (up to 333 g weight loss) of maternal smoking depending on the strength of confounding.

Figure 1

Left plot: estimated expection of birth weight for smokers under not smoking, E ( Y ( 0 ) ∣ T = 1 ) . This estimate is constrained by the sample average estimations of E ( Y ( 1 ) ∣ T = 1 ) in red and E ( Y ( 0 ) ∣ T = 0 ) in blue, as (2) suggests. Right plot: the estimation of E ( Y ( 1 ) ∣ T = 0 ) . See Appendix for a description of the estimators. From these results, plausible values for the sensitivity parameters ρ 0 and ρ 1 are ( − 0.25 , 0.05) and ( − 0.2 , 0.06) respectively; i.e., where crossing with the red and blue bounds occur in respective plots.

Figure 2

Estimates of average birth weight under smoking and non-smoking during pregnancy are given in the right and left plots, respectively; using output from ui.causal function (hdim.ui package). Solid black lines for bias-corrected AIPW estimates, τ ˆ AIPW refit − b ˆ refit , and 95% confidence intervals (gray area). Red intervals are confidence intervals assuming no unobserved confounding ( ρ t = 0 , t = 0,1 ), and blue intervals are 95% uncertainty intervals (using − 0.25 ≤ ρ 0 ≤ 0.05 and − 0.20 ≤ ρ 1 ≤ 0.06 ).

The difference in conclusions may be due to the fact that Nabi et al. [24] analysis excludes the unconfoundedness situation as one of the possible scenarios while ours do not. This discrepancy might be due to different modeling assumptions and thereby different influence functions and corresponding estimation procedures. For instance, Nabi et al. [24] sensitivity model establishes a link between the observed and unobserved potential outcome densities, which necessitates the requirement of common supports; a condition stating that the support of each of the missing potential outcomes must be a subset of the support of the corresponding observed potential outcome. This assumption may affect the validity of the results (see Section 7 in [16]). Furthermore, our method only uses clinical assumptions for bounding sensitivity parameters, whereas Nabi et al. [24] uses those assumptions for both selecting a valid tilting function (defining the sensitivity model) and bounding sensitivity parameters. The effect on their analysis of using alternative tilting functions is unclear to us.

5 Discussion

Unobserved confounding cannot be discarded nor empirically investigated in observational studies, and therefore, sensitivity analysis to the unconfoundedness assumption should be common practice. Moreover, high-dimensional settings are typical in observational studies using large data sets and machine learning for nuisance models. We have presented here a novel method to conduct sensitivity analysis in such situations using uniformly valid estimators, which is essential in high-dimensional setting. In particular, the sensitivity analysis proposed is based on the construction of an uncertainty interval for the causal effect of interest which we show has uniformly valid coverage. Finite sample experiments confirm the asymptotic results.

We use a sensitivity model with a sensitivity parameter, a correlation, which is easy to interpret and discuss with subject-matter scientists [23]. As all sensitivity models, ours describes potential departures from the unconfoundedness assumption. If sensitivity is detected as is the case in the presented application on the effect of smoking on birth weight, then this is important information. If no sensitivity is detected, then one might argue that this does not preclude the analysis to be sensitive to other departures from the unconfoundedness assumption. Note that our results show that we need to let ρ tend to zero with increasing sample size unless we can estimate bias due to unobserved confounding, and hence the propensity score, with a root-n convergence rate. An alternative to bias estimation is to fix the bias itself as the sensitivity parameter, at the expense of interpretability. Finally, future directions for research include generalizing our results from binary to continuous treatments, from continuous to binary outcome, as well as considering ultra-high dimensional situations ( p ≫ n ), which have their own challenges [39].

6 Supplementary material

Supplementary material contains additional simulation results and is available online.