Doubly robust adaptive LASSO for effect modifier discovery

2.1 The framework

The observed data, { ( W i , A i , Y i ) } i = 1 n , are comprised of independent and identically distributed samples of O = ( W , A, Y) ∼ P ₀, where W is the baseline covariates of a patient, A is the binary treatment which equals 1 if the patient received treatment and 0 otherwise, and Y is the observed outcome (binary or continuous). Let V represent the subset of the variables in W that represents the potential EMs of interest. We use O _i = ( W _i , A _i , Y _i ) to represent the ith observation of the data. In order to define the target parameter, we use the counterfactual framework of Rubin [25]. Let Y ^a denote the potential (or counterfactual) outcome that would have occurred under the treatment value A = a. In this paper, we focus on marginal models for the CATE. If we assume that we observe Y = Y ^a when A = a (consistency [26], no interference, positivity and no unmeasured confounders [27]), the CATE can be defined and identified nonparametrically as:

where E ₀ is the expectation with respect to the outcome and E _W|V is the expectation conditional on the baseline covariates. In this work, we choose to model the CATE using a linear regression model defined as ψ ̃ 0 ( V ) = β 0 + V T β V where the relevant subset of V will be selected using adaptive LASSO [24]. Our goal here is to identify the true EMs among the set V , and estimate their associated coefficients. One could use non-linear models or machine learning methods to estimate ψ ̃ 0 ( V ) , which is important when the goal is prediction [37] (e.g. for personalized medicine). However, if interpretation of the coefficient associated with each V ^(s) is important, it may be beneficial to use a linear model rather than a black box approach [28].

2.2 Adaptive LASSO

The adaptive LASSO [24] is an extension of the traditional LASSO of Tibshirani [29] that uses coefficient specific weights. Zou [24] showed that the adaptive LASSO estimator has the oracle property which roughly means that the algorithm identifies the right subset of variables (consistency of variable selection) and that the coefficient estimators of the selected variables are asymptotically normal. In a prediction (non-causal) setting, let Y be an observed outcome and V a set of covariates. Under the linear model, we can select predictors of Y by solving the equation below:

where β ′ = β 1 ′ , … , β p ′ , w ̂ j = 1 / | β ̃ j ′ | γ , for some γ > 0 and β ̃ j ′ is a n -consistent estimator of β j ′ . The selected variables are the positions of the non-zero entries of the solution of (2). When the sample size grows, the weights associated with the zero-coefficient predictors tend to infinity, while the weights corresponding to true predictors converge to a constant. Thus, true-zero coefficients are less likely to be selected by the adaptive LASSO than by the standard LASSO, which does not have the oracle property [24].

2.3 Highly adaptive LASSO (HAL)

Assume E(Y|V) a regression function where Y is the observed outcome and V is the set of covariates. Consider a map of V onto a set of binary indicator basis functions. For example, if V is scalar, we generate for an observation v, ϕ * ( v ) = ( ϕ 1 * ( v ) , … , ϕ n * ( v ) ) T , where ϕ i * ( v ) = I ( v ≥ V i ) , for i = 1, …, n. With two dimensions, V = ( V ( 1 ) , V ( 2 ) ) T , we need to include the second order basis functions ϕ i * ( v ) = I ( v 1 ≥ V i ( 1 ) , v 2 ≥ V i ( 2 ) ) , for i = 1, …, n. The HAL estimator [30] is obtained by fitting a L ₁-penalized regression of the outcome Y on these basis functions, with the optimal L ₁-norm chosen via cross-validation. The HAL estimator of the regression function E(Y| V ) converges to the true regression function in L ₂-norm no slower than n ^−1/4 regardless of the dimension of V , under the assumption that the regression function has bounded variation norm.

2.4 Selective inference

Let β ̂ ′ be the solution of (2) and β ̂ M ̂ ′ the non-zero subvector of β ̂ ′ where M ̂ ⊆ { 1 , … , p } corresponds to the positions of the non-zero entries. Suppose that we are interested in making inference for β ̂ M ̂ ′ in the prediction model of Section 2.2. A naive way to obtain inference after selecting the covariates in the model is the standard hypothesis tests for linear regression that treat M, representing the non-zero entries of β ′ and thus the true model, as known. It is easy to see that β ̂ ′ depends on the selected model M ̂ . Therefore, Lee et al. [31] studied the conditional distribution β ̂ M ′ | { M ̂ = M } and showed that this conditional distribution is a truncated normal Gaussian. They constructed a pivotal statistic for β ̂ M ̂ ′ which can be used for hypothesis testing and therefore by test inversion, to construct a confidence interval. Let F(y; μ, σ ², l, u) be the CDF of a normal N(μ, σ ²) truncated to the interval [l, u], e _j the unit vector for the jth coordinate so that β ̂ M ′ j = η M T Y , η M = V M T V M − 1 V M T T e j and σ * 2 = σ 2 η M T η M . In the linear regression setting where Y ∼ N μ , σ 2 I n , Lee et al. [31] showed that F ( β ̂ M ′ j ; β M ′ j , σ * 2 , ν − , ν + ) | { M ̂ = M } ∼ Unif ( 0,1 ) , where [ν ⁻, ν ⁺] is defined in [31] as a function of Y and the model M. By inverting the hypothesis testing, we can find a (1 − α) confidence interval for β ̂ M ′ j , conditional on M ̂ = M , by finding [L*, U*] such that

and

In this next section, we will explain how this result is applied in our setting.

2.5 The model

3.1 Data generation and parameter estimation

To evaluate the performance of the proposed method in finite samples, we conducted a simulation study under four scenarios. We simulated data O = ( W , A, Y) representing baseline covariates W , a binary exposure A, and a continuous outcome Y. The baseline covariates W include three confounders (X, V ⁽¹⁾, V ⁽²⁾), one instrument Z (pure cause of treatment), and two pure causes of the outcome (V ⁽³⁾, V ⁽⁴⁾). All covariates were generated independently with the Bernoulli distribution with success probability p: X ∼ B(p = 0.4), V ⁽¹⁾ ∼ B(p = 0.5), V ⁽²⁾ ∼ B(p = 0.6), V ⁽³⁾ ∼ B(p = 0.5), V ⁽⁴⁾ ∼ B(p = 0.7) and Z ∼ B(p = 0.45).

We varied the strength of the relationship between covariates, outcome and treatment across three low-dimensional scenarios. In the first, we used an outcome model where the covariates were strongly predictive, and a treatment model where the covariates were weakly predictive. The treatment mechanism g ₀ was set as a Bernoulli with the probability generated linearly in the three confounder variables and single instrument,

where expit(x) = 1/{1 + exp(−x)}. The observed continuous outcome Y was linearly generated as:

The effect modification arises due to interaction between treatment and covariates.

The second scenario has the same data generation except that the coefficient of the interaction term V ⁽¹⁾ V ⁽²⁾ V ⁽³⁾ is 0 instead of 4. In the third scenario, we use an outcome model where the covariates are weakly predictive, and a treatment model where the covariates are strongly predictive. We focus here on the first scenario and describe all other simulations settings and results in the Appendix.

We thus have two EMs (V ⁽¹⁾, V ⁽³⁾), where the first is a confounder and the second is a pure cause of the outcome. In practice, we are not aware of the true data generating mechanism. So we have a potential set of EMs: V = (V ⁽¹⁾, V ⁽²⁾, V ⁽³⁾, V ⁽⁴⁾). Let ψ 0 ( V ) = E P 0 ( Y 1 − Y 0 | V ) be the true (nonparametric) CATE, which we model as an MSM: ψ ̃ 0 ( V ) = β 0 + β 1 V ( 1 ) + β 2 V ( 2 ) + β 3 V ( 3 ) + β 4 V ( 4 ) . Our goal here is to identify among the set V , the true EMs and estimate their associated coefficients. Given the data generated, the true values of the coefficients are β _v = (0.5, 0, 1, 0). We set n = 1000 and then 10 000. We also add a smaller sample size n = 100 with results in the appendix.

To evaluate the performance of our method in high-dimensional settings, we also extend the first scenario by adding 50 pure binary noise covariates (unrelated to treatment or outcome) to our set of covariates, which are included as potential confounders and EMs. The true values of the coefficients in the MSM are thus β _v = (0.5, 0, 1, 0, …, 0).

Under each low-dimensional scenario, we tested our proposed method under four different implementations:

1. Qcgc: Both of the models for Q ̄ and g are correctly specified using generalized linear models (GLMs).

2. Qc: Only the GLM for Q ̄ is correctly specified. g is misspecified using a logistic regression of treatment A on variable X.

3. gc: Only the GLM for g is correctly specified. Q ̄ is misspecified using a GLM of treatment Y on variables A and V ⁽³⁾.

4. HAL: Both Q ̄ and g are estimated using the Highly Adaptive LASSO (HAL) [30, 37]. We use the package default setting.

For comparison, we also tested two implementations of a linear regression model for the outcome to directly assess effect modification:

1. NLin: Linear regression with main terms (treatment and all covariates) and interactions between treatment and covariates. Only first-order interactions were included.

2. CLin: Linear regression with a correctly specified outcome model.

Standard confidence intervals are presented for the linear model case and, in our summary, a p-value of less than 0.05 is used as a criterion for a variable to be selected.

In the higher dimensional scenario, only HAL was used to estimate Q ̄ and g.

3.2 Simulation results

For each scenario, we produced boxplots of the MSM coefficient estimates. We also present the percent selection, the coverage proportion of the confidence intervals and the false coverage rate in order to summarize the average performance of each estimator and implementation. The percent selection for our LASSO method was obtained as the percentage of estimated coefficients that are non-zero throughout the 1000 generated datasets, and for the linear regression, the percentage of p-values < 0.05 . The coverage for each true effect modifier was obtained as the number of times the true model was selected and the corresponding confidence intervals contained the true coefficients, divided by the number of times the true model was selected. For the linear regression, the percent coverage was instead calculated for each coefficient and defined as the proportion of the confidence intervals that contained the true coefficient throughout the 1000 generated datasets. The false coverage rate (FCR) for our LASSO model was obtained as the number of non-covering confidence intervals among the selected coefficients, divided by the number of the selected coefficients throughout the 1000 generated datasets [31].

For the first low-dimensional scenario, Figures 1 and 2 contain the boxplots of the MSM coefficient estimates for the true EMs V ( 1 ) , V 3 and non-EMs (V ⁽²⁾, V ⁽⁴⁾), respectively. Table 1 (in the Appendix) contains the numerical results. As shown in the first two boxplots in Figures 1 and 2, the implementations (1) Qcgc and (2) Qc performed very well. We obtained unbiased estimates and confidence interval coverage that tended to be around 95% as sample size increased as shown in Figure 5. The FCR was close to the optimal 0.05. In the third boxplot, corresponding to implementation (3) gc, where only the propensity score was correctly specified, the estimator was more biased for both sample sizes but had higher coverage rates and lower FCR. In the fourth boxplot where the estimator was implemented with HAL, the estimator performed well across all measures. Overall, as shown in Figure 5, the percent coverage of the true EMs V ( 1 ) , V 3 was around the nominal 95% as sample size increased or when at least the outcome model was correclty specified or machine learning methods were used to estimate both nuisance parameters. In all implementations the true effect-modifiers (V ⁽¹⁾, V ⁽³⁾) were selected around 100 percent of the time except when only the propensity score was correctly specified for the smaller sample size (gc). The percent selection of variables that are not effect-modifiers (V ⁽²⁾, V ⁽⁴⁾) was around 20% for n = 1000. In implementations (1), (2), and (4), the percentage was almost halved for n = 10,000. The FCR was controlled around the nominal 0.05 level in all situations even when only one nuisance model was correctly specified. This supports the double robustness of the proposed estimator and the appropriateness of the post-selection confidence intervals. In implementation (5) NLin, the naive linear model with a misspecified term performed poorly, even when increasing the sample size. On the other hand, when the linear model was correctly specified in implementation (6) CLin, the coefficient estimates were unbiased on average and the coverage was near-optimal. For the two other data generating scenarios described at more length in the Appendix, the results (Tables 2 and 3) look similar to those in the first scenario.

Figure 1:

Simulation results illustrations (data generating scenario 1). Box plots of 1000 MSM coefficients estimates for the true EMs V ( 1 ) , V 3 . The true values of the coefficients are (0.5, 1). Notation: Qcgc: Models for Q ̄ and g are correctly specified, Qc: Q ̄ is correctly specified, gc: g is correctly specified, HAL: Q ̄ and g are estimated with HAL, NLin: Naive linear model, CLin: Correct linear model. %sel: percent selection of a covariate × 100, %cov: coverage rate of the confidence interval of a coefficient estimate × 100, FCR: False coverage rate of the model × 100.

Figure 2:

Simulation results illustrations (data generating scenario 1). Box plots of 1000 MSM coefficients estimates for the non-EMs V ( 2 ) , V 4 . The true values of the coefficients are (0, 0). Notation: Qcgc: Models for Q ̄ and g are correctly specified, Qc: Q ̄ is correctly specified, gc: g is correctly specified, HAL: Q ̄ and g are estimated with HAL, NLin: Naive linear model, CLin: Correct linear model. %sel: percent selection of a covariate × 100, %cov: coverage rate of the confidence interval of a coefficient estimate × 100, FCR: False coverage rate of the model × 100.

Table 4 in the Appendix contains the results with the small sample size n = 100. The performance of the proposed methods decreased across all measures except for V ⁽³⁾ where there was a higher coverage rate when Q ̄ and g were correctly specified or estimated with HAL. The results of the high-dimensional setting are presented in Figures 3 and 4. Q ̄ and g were estimated with HAL. The estimates were taken over 100 generated datasets and look similar to Figures 1 and 2 for the covariates V = (V ⁽¹⁾, V ⁽²⁾, V ⁽³⁾, V ⁽⁴⁾) in common. For the noise covariate coefficients, the estimates, given in the density plot of Figure 4, were unbiased for 0. The noise covariates had a low percent selection (see Table 5). Using median statistics, the noise covariates were selected around 14% of the time and that proportion decreased to 13% as we increased the sample size. The FCR exceeded the nominal 5% level and was around 15%.

Figure 3:

Simulation results for high-dimensional setting (data generating scenario 1). Box plots of MSM coefficients estimates over 100 simulations for the potentials EMs V = (V ⁽¹⁾, V ⁽²⁾, V ⁽³⁾, V ⁽⁴⁾). The true values of the coefficients are (0.5, 0, 1, 0). Notations: HAL: Q ̄ and g are estimated with HAL, %sel: percent selection × 100, %cov: coverage rate × 100, FCR: False coverage rate × 100.

Figure 4:

Illustrations for high-dimensional setting. Box plots of the MSM coefficients estimates over 100 simulations for the 50 noise covariates (for both n = 1000 and n = 10,000). The true values of the coefficients are (0, …, 0).

In summary, Table 1 demonstrates that in low-dimensional settings, the proposed algorithm is able to produce unbiased estimates and control the FCR around the nominal level. In contrast, Table 5 demonstrates that in the context of high-dimensional covariates with many candidate EMs, the FCR is generally much larger than the nominal level. Similar results were obtained by Zhao et al. ([10], Figure 2). In addition, at least some non-EMs were always selected by the algorithm at the sample sizes investigated.

4.1 Data

Our data were obtained from a cohort (Firoozi et al. [38]) of deliveries of pregnant women with asthma in order to study the effect of using inhaled corticosteroids (ICS) during pregnancy on birth weight. The population of interest is pregnant women with mild asthma and a singleton delivery in Québec, Canada between 1998 and 2008, aged ≤ 45 years. For simplicity, We considered only the first delivery for each woman in this period. Asthma severity was defined according to an index that is based on the Canadian Asthma Consensus Guidelines (Cossette et al. [39]). A total of 4707 pregnancies in our database fell into this category. ICS exposure was classified in two categories: “use”(a woman who filled at least one prescription of ICS during pregnancy) and “no use”(a woman who did not fill any prescription of ICS during pregnancy). The outcome of interest is birth weight (continuous in kilograms). We identified a variety of maternal baseline variables. These potential confounders measured in the year before pregnancy include demographic characteristics (e.g. income security provider and place of residence), chronic diseases (e.g. hypertension and diabetes) and variables related to asthma (e.g. at least one hospitalization for asthma, at least one emergency department visit for asthma, and oral corticosteroids). We also included the cumulative daily dose of ICS in the year before pregnancy and sex of the newborn as potential confounders. A full list of measured potential confounders can be found in Table 6 in the Appendix. As we do not know which variables are effect modifiers, we included a wide range of variables in the set V , 22 variables in all. Specifically, these variables were: In the year before pregnancy: at least one dose of inhaled short-acting β ₂-agonists (SABA) taken per week, medication for epilepsy, use of warfarin, use of beta blockers, asthma exacerbation, oral SABA use, oral corticosteroids, leukoteriene-receptor antagonists, intranasal corticosteroids, at least one hospitalization for asthma, at least one emergency department visit for asthma, and welfare recipient; At the start of the pregnancy: chronic obstructive disease, cyanotic heart disease, obesity, uterine disorder, antiphospholipid syndrome, sex of the newborn, rural/non-rural residence indicator, hypertension, diabetes, and chromosomal anomalies.

For our pregnancy cohort, the average treatment effect is the expected difference in the mean counterfactual birth weight if all women were exposed to ICS during pregnancy versus the counterfactual birth weight if all women were not [40]. The target parameters are the coefficients β _j , j = 1, …, 22 of the MSM defined as: ψ ̃ 0 ( V ) = β 0 + ∑ j = 1 22 V ( j ) β j , with V = (V ⁽¹⁾, …, V ⁽²²⁾) the set of potential EMs. Taking the sex of the newborn as an EM for example (V ^(j) = sex), β _j is the difference in the CATE for women having male versus female children.

4.2 Results

Baseline characteristics of the pregnancy cohort are presented in Table 6. We first implemented a standard linear regression with main terms for all potential confounders and interaction terms between the treatment and the set V . The estimates of the coefficients of the interaction terms are given in Table 7. A variable was considered to be selected as an EM in the standard linear regression if the coefficient of the interaction term between that variable and the treatment had a p-value < 0.05 . This model concluded that leukoteriene-receptor antagonists and chromosomal anomalies are EMs. In addition, we implemented our LASSO methods using HAL for the estimation of the outcome expectation and propensity score. All of the covariates were included in the propensity score model as well as in the outcome model. Due to larger weights, a 5% truncation for the values of g _n was used. The selected coefficients of the MSM and their estimated values are presented in Table 8. Three covariates (leukoteriene-receptor antagonists, warfarin one year before pregnancy, and chromosomal anomalies) were selected using the adaptive LASSO and two of them were significant (leukoteriene-receptor antagonists and chromosomal anomalies) using post-selection inference. Leukoteriene-receptor antagonists and chromosomal anomalies were thus selected as EMs in the association of taking ICS during pregnancy on birth weight. Although the naive linear model and our algorithm generate very similar sets of EMs, the coefficients of the selected EMs are different (compare Table 7 with Table 8). For example, the estimated coefficient of leukoteriene-receptor antagonist is around −0.17 in the adaptive LASSO while it is −0.365 using the linear model.