A Simple Method for Estimating Interactions between a Treatment and a Large Number of Covariates

2.1 Continuous Response Model

When Y is a continuous response, a simple multivariate linear regression model for characterizing the interactions between the treatment and covariates is

where ε is the mean zero random error. In this simple model, the interaction term ${\mathit{\gamma }}_0^{\prime }\mathbf{W}(Z)·T$ models the heterogeneous treatment effect across the population and the linear combination of ${\mathit{\gamma }}_0^{\prime }\mathbf{W}(\mathbf{Z})$ can be used for identifying the subgroup of patients who may or may not benefit from the treatment. Note that since the vector W(Z) contains an intercept, the main effect for treatment is always included in the model. Specifically, under model (1), we have

i.e., ${\mathit{\gamma }}_0^{\prime }\mathbf{W}(\mathbf{z})$ measures the causal treatment effect for patients with the baseline covariate z. With observed data, γ₀ can be estimated via the ordinary least squares (OLS) method.

On the other hand, noting the relationship that

one may estimate γ₀ by directly minimizing

We call this the modified outcome method, where 2YT can be viewed as the modified outcome, which has been first proposed in the Ph.D thesis of James Sinovitch, Harvard University.

Under the simple linear model (1), estimators from both methods are consistent for γ₀, and the full least squares approach in general is more efficient. In practice, the simple multivariate linear regression model is often just a working model approximating the complicated underlying probabilistic relationship among the treatment, baseline covariates and outcome variables. It comes as a surprise that even when model (1) is mis-specified, the multivariate linear regression and modified outcome estimators still converge to the same non-random limit γ^* determined by the joint distribution of (Y⁽¹⁾, Y⁽⁻¹⁾, Z). Furthermore, the score W(z)′γ^* is a sensible estimator for the interaction effect in the sense that it seeks the “best” function of z in a functional space to approximate Δ(z) by minimizing

where the expectation is with respect to Z.

2.2 The Modified Covariate Method

The modified outcomes estimator defined above is useful for the Gaussian case, but does not generalize easily to more complicated models. Hence we propose a new estimator which is equivalent to the modified outcomes approach in the Gaussian case and extends easily to other models. This is the main proposal of this paper.

We consider the simple working model

where ε is the random error. Based on model (3), we propose the modified covariate estimator γ̂ as the minimizer of

That is, we simply multiply each component of W_i by one-half the treatment assignment indicator (= ±1) and perform a regular linear regression. Now since

the modified outcome and modified covariate estimates are identical and share the same causal interpretations. Operationally, we can perform a simple linear regression with the modified covariates without intercept. We summarize the general proposal below.

Figure 2 illustrates how the modified covariate method works for a single covariate Z. The raw data are shown on the left and the data with modified covariate are shown on the right. The regression line in the right panel estimates the personalized treatment Δ(z).

The advantage of this new approach is twofold: it avoids having to directly model the main effects and it has a causal interpretation for the resulting estimator regardless of the adequacy of the assumed working model (3). Furthermore, unlike the modified outcome method, it is straightforward to generalize the new approach to other types of outcome.

2.3 Binary Responses

When Y is a binary response, in the same spirit as the continuous outcome case, we propose to fit a logistic regression model with modified covariates W^* = W(Z) · T/2 :

If model (7) is correctly specified, then

and thus ${\mathit{\gamma }}_0^{\prime }\mathbf{W}(\mathbf{z})$ has an appropriate causal interpretation. However, even when model (7) is not correctly specified, we still can estimate γ₀ by treating (7) as a working model. In general, the maximum likelihood estimator (MLE) of the working model converges to a deterministic limit γ^* and W(z)′γ^*/2 can be viewed as the solution to the following optimization problem

where the expectation is with respect to (Y, T, Z). Therefore, if W(z) forms a “rich” set of basis functions, W(z)′γ^*/2 is an approximation to the minimizer of E{Yf(Z)T − log(1 + e^f(Z)T)}. In the Appendix, we show that the latter can be represented as

under very general assumptions. Therefore,

may serve as an estimate for the covariate-specific treatment effect and used to stratify the patient population, regardless of the validity of the working model assumptions.

As described above, the MLE from the working model (7) can always be used to construct a surrogate to the personalized treatment effect measured by the “risk difference”

On the other hand, different measures for individualized treatment effects may also be of interest. In supplementary materials we have provided an alternative approach to approximate the personalized treatment effect using the “relative risk” rather than “risk difference”.

2.4 Survival Responses

When the outcome variable is survival time, we often do not observe the exact outcome for every subject in a clinical study due to incomplete follow-up. In this case, we assume that the outcome Y is a pair of random variables (X, δ) = {X˜ ∧ C, I(X˜ < C)}, where X˜ is the survival time of primary interest, C is the censoring time and δ is the censoring indicator.

Firstly, we propose to fit a Cox regression model with modified covariates

where λ(t|·) is the hazard function for survival time X˜ and λ₀(·) is a baseline hazard function free of Z and T. When model (8) is correctly specified,

and ${\mathit{\gamma }}_0^{\prime }\mathbf{W}(\mathbf{z})$ can be used to stratify the patient population according to Δ(z), where ${\mathrm{\Lambda }}_0(t)={\int }_0^t{\lambda }_0(u)du$ is a monotone increasing function. Under the proportional hazards assumption, the maximum partial likelihood estimator γ̂ is a consistent estimator for γ₀ and semiparametric efficient. Moreover, even when model (8) is mis-specified, we still can “estimate” γ₀ by maximizing the partial likelihood function. In general, the resulting estimator, γ̂, converges to a deterministic limit γ^*, which is the root of a limiting score equation [Lin and Wei, 1989]. More generally, W(z)′γ^*/2 can be viewed as the solution of the optimization problem

where N(t) = I(X˜ ≤ t)δ, τ is a fixed time point such that P(X ≥ τ) > 0, and the expectations are with respect to (Y, T, Z). Therefore, W(z)′γ^*/2 can be viewed as an approximation to

In Appendix, we show that if the censoring time does not depend on (Z, T), then the minimizer f^* satisfies

for a monotone increasing function Λ^*(u). Thus, when censoring rates are balanced between the two arms, e.g., the hazard rate is low and most of the censoring is due to administrative reasons,

can be used to characterize the covariate-specific treatment effect and stratify the patient population even when the working model (8) is mis-specified.

2.5 Regularization for the High Dimensional Data

When the dimension of W^*, p, is high, we can easily apply appropriate variable selection procedures based on the corresponding working model. For example, L₁ penalized (Lasso) estimators proposed by Tibshirani [1996] can be directly applied to the modified data (6). In general, one may estimate γ by minimizing

where ${\left|\right|\gamma \left|\right|}_1={\sum }_{j=1}^p\mid {\gamma }_j\mid$ and

The resulting Lasso regularized estimator may share the appealing finite sample oracle properties [Van De Geer, 2008, Zhang and Huang, 2008, Van De Geer and Buhlmann, 2009, Negahban et al., 2012]. For example, for continuous outcomes, if we select the penalty parameter λ_n = O({log(p)/N}^1/2) and the covariates W(Z_i)T_i/2, i = 1, ···, N satisfy the restricted eigenvalue condition, then

where ${\left|\right|\mathit{\gamma }\left|\right|}_2^2={\sum }_{j=1}^p{\gamma }_j^2$, γ̂_λN is the corresponding Lasso regularized estimator, #γ^* is the number of non-zero components in the vector γ^*, and c_i, i = 1, 2, 3 are positive constants depending on the joint distribution of (Y, W(Z), T) [Negahban et al., 2012]. Therefore, although the dimension p may be big, γ^* can always be approximated by the Lasso estimator with an approximation error proportional to only the number of its non-zero components, which is small if γ^* is sparse. Following the Corollary 3 of Negahban et al. [2012], similar results hold when γ^* is not exactly sparse but can be approximated by a sparse vector. One consequence is that if W(·) is adequately rich such that |W(z)′γ^* − f^*(z)| is small and γ^* is either exact or approximately sparse, then ${\widehat{\mathit{\gamma }}}_{{\lambda }_N}^{\prime }\mathbf{W}(\mathbf{z})$ is a good approximation of f^*(z) with a high probability, where f^*(z) is a monotone transformation of the individualized treatment effect Δ(z). This property holds for any pair of finite N and p with a sufficiently small ratio log(p)/N and does not depend on the correct specification of the working model. A similar result for the expected utility in the patient population with the estimated optimal treatment assignment rule has been established in Qian and Murphy [2011]. This result can be extended to the binary case, with appropriate regularity condition on the covariates W(Z_i)T_i/2 to ensure that the negative log-likelihood function satisfies a form of restricted strong convexity [Negahban et al., 2012]. For the survival case, parallel results have been established for Lasso-regularized maximum partial likelihood estimator under the assumption that the multiplicative Cox model is correctly specified [Kong and Nan, 2013, Huang et al., 2013]. The finite sample properties of the Lasso estimator for mis-specified Cox model are still unknown and warrants further study.

Note that the Lasso regularization method may not be able to consistently recover all the nonzero components of γ^* [Zhao and Yu, 2006, Zou, 2006]. However, since our focus is on approximating the personalized treatment effect under a potentially misspecified working model, where the regression coefficients γ^* may not have a clinically meaningful interpretation, this limitation is not necessarily a severe drawback of the proposed method.

Regularization methods other than Lasso especially those with nonconvex penalties such as SCAD and MCP may also be used [Fan and Li, 2001, Zhang, 2010]. The composite gradient descend algorithm and the adaption thereof developed recently can be used to compute the corresponding regularized optimizer efficiently [Loh and Wainwright, 2012]. In addition, in contrast to Lasso, some nonconvex regularization methods may provide consistent estimates for the support of γ^* under suitable regularity conditions if the working model is correctly specified. [Zou, 2006, Huang and Xie, 2007, Huang et al., 2008, Zhang and Zhang, 2012].

2.6 Efficiency Augmentation

When the models (5, 7 and 8) with modified covariates are correctly specified, the MLE estimator for γ^* is the most efficient estimator asymptotically. However, when models are treated as working models subject to mis-specification, a more efficient estimator can be obtained for estimating the same γ^*. To this end, note that in general γ̂ is defined as the minimizer of an objective function motivated from a working model:

Since for any function a(z) : R^q → R^p, E{T_ia(Z_i)} = 0 due to the randomization, the minimizer of the augmented objective function

converges to the same limit as γ̂. Furthermore, by selecting an optimal augmentation term a(·), the minimizer of the augmented objective function may have smaller variance than that of the original estimator. In supplementary materials, we show that

are optimal choices for continuous and binary responses, respectively. Therefore, we propose the following two-step procedures for estimating γ^* :

1. Estimate the optimal a₀(z) :

   1. For continuous responses, fit the linear regression model E(Y|Z) = ξ′B(Z) for the appropriate function B(Z). Let

      $$\widehat{\mathbf{a}}(\mathbf{z})=-\frac{1}{2}\mathbf{W}(\mathbf{z})\times {\widehat{\xi }}^{\prime }\mathbf{B}(\mathbf{z}).$$

   2. For binary response, fit the logistic regression model logit{P(Y = 1|Z)} = ξ′B(Z) for the appropriate function B(Z). Let

      $$\widehat{\mathbf{a}}(\mathbf{z})=-\frac{1}{2}\mathbf{W}(\mathbf{z})\times \left\{\frac{e^{\widehat{{\xi }^{\prime }}\mathbf{B}(\mathbf{z})}}{1+e^{{\widehat{\xi }}^{\prime }\mathbf{B}(\mathbf{z})}}-\frac{1}{2}\right\}.$$

   Here as W(z), B(z) is a transformation of z selected by users.

2. Estimate γ^*

   1. For continuous responses, we minimize (appropriately regularized)

      $$\frac{1}{N}\sum _{i=1}^N\left\{\frac{1}{2}{(Y_i-{\mathit{\gamma }}^{\prime }{\mathbf{W}}_i^{\ast })}^2-{\mathit{\gamma }}^{\prime }\widehat{\mathbf{a}}({\mathbf{Z}}_i)T_i\right\}.$$

   2. For binary responses, we minimize (appropriately regularized)

      $$\frac{1}{N}\sum _{i=1}^N\left[-\{Y_i{\mathit{\gamma }}^{\prime }{\mathbf{W}}_i^{\ast }-log(1+e^{{\mathit{\gamma }}^{\prime }{\mathbf{W}}_i^{\ast }})\}-{\mathit{\gamma }}^{\prime }\widehat{\mathbf{a}}({\mathbf{Z}}_i)T_i\right].$$

For survival outcome, the log-partial likelihood function is not a simple sum of i.i.d terms. However, in supplementary materials, we show that the optimal choice of a(z) is

where G₁(u; z) = E{M(u)|Z = z, T = 1}, G₂(u; z) = E{M(u)|Z = z, T = −1},

Unfortunately, a₀(z) depends on the unknown parameter γ^*. On the other hand, for the high-dimensional case, the interaction effect is usually small and it is not unreasonable to assume that γ^* ≈ 0. Furthermore, if the censoring patterns are similar in both arms, we have G₁(u, z) ≈ G₂(u, z). Using these two approximations, we can simplify the optimal augmentation term as

where

Therefore, we propose to employ the following approach to implement the efficiency augmentation procedure:

1. Calculate

   $${\widehat{M}}_i(\tau )=N_i(\tau )-{\int }_0^{\tau }\frac{I(X_i\ge u)d\{{\sum }_{j=1}^N{N}_j(u)\}}{{\sum }_{j=1}^{N}I(X_j\ge u)},i=1,\cdots ,N$$

   and fit the linear regression model E(M̂(τ)|Z) = ξ′B(Z). Let

   $$\widehat{\mathbf{a}}(\mathbf{z})=-\frac{1}{2}\mathbf{W}(\mathbf{z})\times {\widehat{\xi }}^{\prime }\mathbf{B}(\mathbf{z}).$$

2. Estimate γ^* by minimizing (appropriately regularized)

   $$\frac{1}{N}\sum _{i=1}^N\left(-\left[{\mathit{\gamma }}^{\prime }{\mathbf{W}}_i^{\ast }-log\{\sum _{j=1}^N{e}^{{\mathit{\gamma }}^{\prime }{\mathbf{W}}_i^{\ast }}I(X_j\ge X_i)\}\right]{\mathrm{\Delta }}_i-{\mathit{\gamma }}^{\prime }\widehat{\mathbf{a}}({\mathbf{Z}}_i)T_i\right).$$

3.1 Continuous Responses

For continuous responses, we generated N independent Gaussian samples from the regression model

where the covariates (Z₁, ···, Z_p) follow a mean zero multivariate normal distribution with a compound symmetric variance-covariance matrix, (1 − ρ)I_p + ρ1′1, and ε ~ N(0, 1). We let (γ₀, γ₁, γ₂, γ₃, γ₄, γ₅, ···, γ_p) = (0.4, 0.8, −0.8, 0.8, −0.8, 0, ···, 0), α_ij = 0.8I(i = 1, j = 2), ${\sigma }_0=\sqrt{2}$, N = 100, and p = 50 and 1000 representing low and high dimensional cases, respectively. The treatment T was generated as ±1 with equal probability at random. We consider four sets of simulations:

1. ${\beta }_0={(\sqrt{6})}^{-1},{\beta }_j={(2\sqrt{6})}^{-1}$, j = 3, 4, ···, 10 and ρ = 0;

2. ${\beta }_0={(\sqrt{6})}^{-1},{\beta }_j={(2\sqrt{6})}^{-1}$, j = 3, 4, ···, 10 and ρ = 1/3;

3. ${\beta }_0={(\sqrt{3})}^{-1},{\beta }_j={(2\sqrt{3})}^{-1}$, j = 3, 4, ···, 10 and ρ = 0;

4. ${\beta }_0={(\sqrt{3})}^{-1},{\beta }_j={(2\sqrt{3})}^{-1}$, j = 3, 4, ···, 10 and ρ = 1/3.

Settings 1 and 2 represented cases with relatively small main effects, where the variations in responses contributable to the main effect, interaction and random error were about 37.5%, 37.5% and 25%, respectively, when the covariates were correlated. Settings 3 and 4 represented cases with relatively big main effects, where the variations in responses contributable to the main effect, interaction and random error were about 75%, 15% and 10%, respectively, when the covariates were correlated. For each of the simulated data set, we implemented three methods:

• full regression: Fit a multivariate linear regression with complete main effects and covariate/treatment interaction terms, i.e., the dimension of the covariate matrix was 2(p + 1). The Lasso was used to select the variables.

• new: Fit a multivariate linear regression with the modified covariate W^* = (1, Z)′ · T/2. The dimension of the covariate matrix was p + 1. Again, the Lasso was used.

• new/augmented: The proposed method with efficiency augmentation, where E(Y|Z) was estimated with Lasso-regularized ordinary least squared method and B(z) = z.

The selection of W(z) = B(z) = z mimicked the realistic situation where the knowledge of the true functional forms for interaction and main effects were often lacking. For all three methods, we selected the Lasso penalty parameter via 10-fold cross-validation. The outputs of all three methods are estimated scores, which are linear combinations of W(z) approximating the personalized treatment effect Δ(z) or a transformation thereof. Specifically, the score is γ̂′W(z), where γ̂ is simply the estimated coefficients for the covariate-treatment interaction terms in the full regression approach and the estimated coefficients for modified covariates in our proposal. The estimated scores can be used to rank patients according to their personalized treatment effects and select a subgroup of patients who may benefit from the treatment. Therefore, to evaluate the performance of the resulting score, we estimated the Spearman’s rank correlation coefficient between the estimated score and the “true” treatment effect

in an independently generated set with a sample size of 10000. Based on 500 sets of simulations, we plotted the boxplots of the rank correlation coefficients between the estimated scores γ̂′W(Z) and Δ(Z) under simulation settings (1), (2), (3) and (4) in the top left, top right, bottom left and bottom right panels of Figure 3, respectively. For all settings, the performance of the modified covariates method was better than that of the full regression approach with further boosting from the efficiency augmentation. The superiority of the new method was fairly obvious especially when p = 1000. For example, in the setting 3, where all the covariates were independent and the main effect was relatively big, the median correlation coefficients were 0.15, 0.46 and 0.51 for the full regression, new and new/augmented methods, respectively. Furthermore, the results on identifying covariates interacting with the treatment are reported at Tables 1 and 2 for p = 50 and 1000, respectively. The new method was also superior in terms of selecting the right covariates. For example, in the setting 3 with p = 1000, while the new method with augmentation on average selected 25.2 covariates with 2.9 true positives, i.e. covariates from Z_i, 1 ≤ i ≤ 4, the full regression method only selected 4.9 covariates with 0.9 true positive.

3.2 Binary Responses

For binary responses, we used the same simulation design as that for the continuous response. Specifically, we generated N independent binary samples from the regression model

where all the model parameters were the same as those in the case of continuous responses. Noting that the logistic regression model was mis-specified under the chosen simulation design. We also considered the same four settings with different combinations of β_j and ρ. For each of the simulated data set, we implemented three methods as those for the continuous outcome: full regression, new and new/augmented. The logistic regression working model is used for all three approaches. In the efficiency augmentation, E(Y|Z) was estimated with the Lasso-penalized logistic regression.

To evaluate the performance of the resulting score measuring the individualized treatment effect, we estimated the Spearman’s rank correlation coefficient between the estimated score and the “true” treatment effect

Although the scores measuring the interaction from the first and second/third methods were different even when the sample size goes to infinity, the rank correlation coefficients put them on the same footing in comparing performances.

In the top left, top right, bottom left and bottom right panels of Figure 4, we plotted the boxplots of the correlation coefficients between the estimated scores γ̂′W(Z) and Δ(Z) under simulation settings (1), (2), (3) and (4), respectively. The patterns were similar to those for the continuous response and confirmed our findings that the “efficiency-augmented method” performed the best or close to the best in all four settings. For example, in the setting 3, the median correlation coefficients were 0, 0.23 and 0.29 for the full regression, new and new/augmented methods, respectively, when p = 1000. In addition, while the proposed method with augmentation on average selected 13.2 covariates with 1.5 true positives, the full regression method only selected 2.0 covariates with 0.3 true positive in the same setting. Here the true positive covariates of interest are again Z_i, 1 ≤ i ≤ 4. The detailed results on covairates identification were given at Tables 1 and 2.

3.3 Survival Responses

For survival responses, we used the same simulation design as that for the continuous and binary responses. Specifically, we generated N independent survival times from the regression model

where all the model parameters were the same as those in the previous subsections. The censoring time was generated from the uniform distribution U(0, ξ₀), where ξ₀ was selected to induce a censoring rate of 25%. For each simulated data set, we implemented the following three methods as above: full regression, new and new/augmented. The Cox regression working model is used for all three methods. In the efficiency augmentation, E{M(τ)|Z} is estimated with the Lasso-regularized linear regression.

To evaluate the performance of the resulting score measuring the individualized treatment effect, we estimated the Spearman’s rank correlation coefficient between the estimated score and the “true” treatment effect based on survival probability at t₀ = 20

In the top left, top right, bottom left and bottom right panels of Figure 5, we plotted the boxplots of the correlation coefficients between the estimated scores γ̂′W(Z) and Δ(Z) under simulation settings, (1), (2), (3) and (4), respectively. The patterns were similar to those for the continuous and binary responses. For example, when p = 1000, the median correlation coefficients were 0.00, 0.40 and 0.39 in the setting 3 for the full regression, new and new/augmented methods, respectively.

Lastly, we also repeated the aforementioned numerical studies for continuous, binary and survival outcomes with the Lasso replaced by adaptive Lasso, whose empirical performance is often similar to that of nonconvex regularization such as SCAD and MCP. The simulation results (not shown here) are very similar to those presented above, demonstrating that the modified covariate approach can be coupled with regularization methods other than Lasso for the purpose of dimension reduction in practice.