Regularized Outcome Weighted Subgroup Identification for Differential Treatment Effects

1. Introduction

Subgroup identification for differential outcomes is becoming an increasingly important methodological research topic in this era of personalized medicine and high medical costs (Ruberg, Chen, and Wang, 2010). Identification of subgroups via valid predictive markers can facilitate tailored intervention to maximize effectiveness. The controversy and challenges of subgroup identification have been well documented (Brookes et al., 2004; Lagakos, 2006; Ruberg et al., 2010). Yet there are still great needs for descriptive and exploratory subgroup identification using existing data to develop new hypothesis, to facilitate study population selection for further investigation, and to generate new evidence for medical decision making (Varadhan et al., 2013).

The classical approach to identify subgroups involves fitting of a parametric outcome model which includes covariate main effects as well as interactions between covariates and treatment. Subgroups are identified based on significance of specific interaction terms (Kehl and Ulm, 2006). Besides the multiplicity issue, misspecification of the main effects in the presence of many covariates interferes with the covariate–treatment interaction estimation and thus impedes valid subgroup identification. Nonparametric approaches based on the classification and regression tree methodology (Breiman et al., 1984) were attempted to separate the main effects from the covariate–treatment interaction effects, either through multiple testing strategies (Su et al., 2008, 2009; Lipkovich et al., 2011; He, 2012) or prediction (Foster, Taylor, and Ruberg, 2011). For example, Su et al. (2009) developed an interaction tree method by building splitting rules based on covariate– treatment interaction tests. Foster et al. (2011) developed the virtual twins method by first estimating treatment differences based on random forests and then applying tree to those estimates to form subgroups. Alternatively, Cai et al. (2011) built in calibration steps to alleviate the model misspecification problem when using parametric or semiparametric models.

In this article, we propose an approach that works directly with differential treatment effect estimation. Instead of modeling the outcome, we use a causal inference formulation (Zhao et al., 2012) that leads to a target function directly reflecting the covariate–treatment interaction effect. This spares the need to model the covariate main effects. The function uses patient outcomes as weights rather than modeling targets. Consequently, our method can deal with binary, continuous, time-to-event, and possibly contaminated outcomes in the same fashion, and can be implemented using existing packages. Foster et al. (2011) pointed out that instead of focusing on finding the optimal or largest subgroup with enhanced treatment effects, it may be more desirable to determine a sub-optimal subgroup from simple and interpretable predictive variables. Similar to Imai and Ratkovic (in press), we focus on linear rules and build in smoothness constraints that allow the resulting subgroup to be more interpretable. Recent applications often involve a large number of covariates. Considering covariate–treatment interaction terms also increases the dimension of the covariate matrix. Thus, we add a regularization step in the subgroup identification process. The framework and details of our methodology are given in Section 2. We note that such regularization has been used in dynamic treatment regimens (Song et al., in press).

Our motivating example is a mammography screening randomized trial that tested the efficacy of interventions to promote mammography screening among eligible women who were 51 years or older (Champion et al., 2007). Data were collected from 09/01/96 to 11/30/02. Different from the usual drug trial, this study intended to change health behavior of the study subjects, which was more complicated and challenging. Besides establishing an overall efficacy, it was also highly desirable to determine whether the interventions were more effective for certain subgroups determined by their social and behavioral variables. For this reason, in addition to demographics and social economic variables, the study collected many health behavioral variables such as self-efficacy, susceptibility, knowledge, fear, fatalism, benefits, and barriers based on the well-established Health Belief Model (Janz and Becker, 1984). These variables intended to capture corresponding behavioral constructs. If interventions benefit subjects differently according to their characteristics, then they can be applied more cost-effectively. In addition, newer interventions can be designed to target relevant behavioral constructs and other modifiable characteristics. It is therefore highly desirable that the resulting treatment assignment rules are interpretable and simple. For this purpose, we focus on finding one-way and two-way interactions among the covariates and interventions. Nevertheless, a large number of variables still need to be investigated due to the many levels of the behavioral variables, prompting us to focus on parsimonious and interpretable models.

In Section 2, we argue that for many popular models, it suffices to find the optimal linear rule. The linearity allows regularization that can lead to sparse solutions, in contrast with Zhao et al. (2012). In Section 3, we develop an inference procedure on the enhanced treatment effects to further evaluate the benefit of the identified subgroup. Such a step is important to aid practical decision making. Our method is compared with some other recent subgroup identification methods, including a sparse version of Zhao et al. (2012), through simulation studies in Section 4, and is applied to two real life examples in Section 5. We conclude with a discussion of our findings and future work in Section 6.

2. Subgroup Identification

Assume data are collected from a randomized trial with a binary treatment T ∈ {−1, 1 being according to P(T = 1) = π a clinical outcome Y, and a p-dimensional covariates vector X = (X₁,…, X_p)′. We allow X to contain derived covariates such as interaction terms as well as dummy variables for categorical covariates. Let $\mathcal{D}(\mathit{X})$ be an assignment rule based on X, P the distribution of (Y, T, X) in the data and ${\mathrm{P}}^{\mathcal{D}}$ the distribution of (Y, T, X) given that $T=\mathcal{D}(\mathit{X})$. Then, the expected outcome ${\mathrm{E}}^{\mathcal{D}}(Y)$ under the rule $\mathcal{D}$ is given by (Qian and Murphy, 2011; Zhao et al., 2012; Zhang et al., 2012)

where I(·) is the indicator function. Our main goal is to assign each individual to an appropriate treatment based on X to optimize the clinical outcome. Assuming larger Y is preferable, we want to find the optimal rule ${\mathcal{D}}^{\ast }$ that maximizes ${\mathrm{E}}^{\mathcal{D}}(Y)$, or equivalently

In general, the optimal rule ${\mathcal{D}}^{\ast }(\mathit{X})$ could be a very complicated function of X. However, we show in the following proposition that, under mild conditions, the optimal rule ${\mathcal{D}}^{\ast }$ only depends on the sign of a linear function of X.

3. Inference on Enhanced Comparative Treatment Effects

To evaluate the effect of subgroup identification, we consider the following quantity measuring the improvement of response when assigning subgroups of subjects to appropriate treatments:

The actual comparative treatment effect for the subgroup identified for T = 1 via $\widehat{\mathcal{D}}(\mathit{X})=1$ is reflected by $d(\mathit{X},1,\widehat{\mathit{\beta }})\equiv d(\mathit{X},T=1,\widehat{\mathit{\beta }})$, and the mean effect is $E_{\mathit{X}}\{d(\mathit{X},1,\widehat{\mathit{\beta }})\}$ for the identified subgroup. The expectation E_X is with respect to X and we treat βˆ as “fixed” for future studies, similar to Foster et al. (2011). Correspondingly, the actual comparative treatment effect for the subjects identified for T = −1 is reflected by $d(\mathit{X},-1,\widehat{\mathit{\beta }})$, and the mean effect is. $E_{\mathit{X}}\{d(\mathit{X},-1,\widehat{\mathit{\beta }})\}$ Ideally, both $E_{\mathit{X}}\{d(\mathit{X},1,\widehat{\mathit{\beta }})\}$ and $E_{\mathit{X}}\{d(\mathit{X},-1,\widehat{\mathit{\beta }})\}$ should be as large as possible. However it is also likely that one of them is small, for example when T = –1 is an inactive control. Inference about both $E_{\mathit{X}}\{d(\mathit{X},1,\widehat{\mathit{\beta }})\}$ and $E_{\mathit{X}}\{d(\mathit{X},-1,\widehat{\mathit{\beta }})\}$ is important and can aid practical decision making regarding the identified subgroups.

In the following, we propose a bootstrap method to construct confidence intervals (CIs) for $E_{\mathit{X}}\{d(\mathit{X},1,\widehat{\mathit{\beta }})\}$and $E_{\mathit{X}}\{d(\mathit{X},-1,\widehat{\mathit{\beta }})\}$. If such intervals only contain positive numbers, we conclude that subgroup identification is indeed beneficial for treatment enhancement. We consider cases when P(X′β^† = 0) = 0 and when P(X′β^† = 0) > 0. The first case is regular and we can use the usual bootstrap. The second case is irregular (Shao, 1994; Laber and Murphy, 2011) due to a non-negligible probability concentrated at the discontinuous points of $d(\mathit{X}\mathit{,}T\mathit{,}\widehat{\mathit{\beta }})$ as a function of $\widehat{\mathit{\beta }}$. In this case, we adopt the m-out-of-n bootstrap (Shao, 1994; Bickel, Götze, and van Zwet, 1997; Chakraborty, Laber, and Zhao, 2013). The method uses the bootstrap sample size m smaller than n to asymptotically avoid the non-negligible probability concentrated at the discontinuous points. Usually m is chosen as n^α for some α < 1. Chakraborty et al. (2013) used α = 0.8. We also adopt this choice in our simulation and data analysis. Some preliminary sensitivity analysis on different choices of α is presented in the Web Appendix D. In practice, if many covariates are continuous, we can use the ordinary bootstrap (Moodie and Richardson, 2010). If many covariates are categorical, we propose to use m-out-of-n bootstrap.

Let {(Y_i, T_i, X_i), i = 1,…, n} be the original dataset. The m-out-of-n bootstrap algorithm comprises the following two steps.

• Step 1: Generate a bootstrap sample with size $m,\{(Y_k^{(b)},T_k^{(b)},{\mathit{X}}_k^{(b)}),k=1,\dots ,m\}$, from the original dataset. Obtain the corresponding estimates ${\widehat{\mathit{\beta }}}^{(b)}$. Let $A_1=\{k|{\mathit{X}}_k^{\prime }{\widehat{\mathit{\beta }}}^{(b)}\ge 0,T_k=1\}$, $A_2=\{k|{\mathit{X}}_k^{\prime }{\widehat{\mathit{\beta }}}^{(b)}\ge 0,T_k=-1\}$, $A_3=\{k|{\mathit{X}}_k^{\prime }{\widehat{\mathit{\beta }}}^{(b)}<0,T_k=-1\}$, $A_4=\{k|{\mathit{X}}_k^{\prime }{\widehat{\mathit{\beta }}}^{(b)}<0,T_k=1\}$and |A_j| be the size of the set A_j for j = 1, 2, 3, and 4. Calculate $\stackrel{\sim }{d}(1,{\widehat{\mathit{\beta }}}^{(b)})={\displaystyle {\sum }_{k\in A_1}{Y}_k/|A_1|-}{\displaystyle {\sum }_{k\in A_2}{Y}_k/|A_2|}$ and $\stackrel{\sim }{d}(-1,{\widehat{\mathit{\beta }}}^{(b)})={\displaystyle {\sum }_{k\in A_3}{Y}_k/|A_3|-{\displaystyle {\sum }_{k\in A_4}{Y}_k|A_4|}}$.

• Step 2: Repeat Step 1 B times. Let ${\widehat{l}}_{+}$ and ${\widehat{u}}_{+}$ be the α/2 and 1−α/2 quantiles of $\{\stackrel{\sim }{d}(1,{\widehat{\mathit{\beta }}}^{(b)}),b=1,\dots ,B\}$, and ${\widehat{l}}_{-}$ and ${\widehat{u}}_{-}$corresponding quantiles of $\{\stackrel{\sim }{d}(-1,{\widehat{\mathit{\beta }}}^{(b)}),b=1,\dots ,B\}$ respectively. Then, the (1 − α) CI of $E_{\mathit{X}}\{d(\mathit{X},1,\widehat{\mathit{\beta }})\}$ is given by $[{\widehat{l}}_{+},{\widehat{u}}_{+}]$ and the (1 − α) CI of $E_{\mathit{X}}\{d(\mathit{X},-1,\widehat{\mathit{\beta }})\}$is given by $[{\widehat{l}}_{-},{\widehat{u}}_{-}]$.

4. Numerical Results

In this section, we consider binary outcomes. Some of our simulation studies for time-to-event outcomes are presented in the Web Appendix B. We numerically compare our proposed regularized outcome weighted subgroup identification (ROWSi) method with some other recent subgrouping methods: (1) the FindIt by Imai and Ratkovic (in press), (2) the virtual twins (VT) by Foster et al. (2011), (3) the interaction tree (IT) by Su et al. (2009), (4) the logistic regression with LASSO penalty (logLASSO) by Qian and Murphy (2011), 5) the L₁-penalized support vector machine (L1-SVM) method adapted from Zhao et al. (2012).

To implement the FindIt, we used the R package FindIt developed by Imai and Ratkovic with optimal tuning parameters chosen by their algorithm. We only considered two-way or three-way interactions among the treatment and covariates. To implement the VT, we first used the R function<monospace>randomForest</monospace>to estimate individual treatment effect $\widehat{P}(y_i=1|T_i=1,{\mathit{X}}_i)-\widehat{P}(Y_i=1|T_i=-1,{\mathit{X}}_i)$ for subject i, and the average treatment effect $\widehat{P}(y=1|T=1,)-\widehat{P}(Y=1|T=-1)$ Then, we assigned subjects to treatment T = 1 if their individual treatment effects were greater than the average effect plus a constant, 0.05, as in Foster et al.(2011). We also followed Foster et al.(2011) for the choices of the terminal node size and complexity parameter. For the IT method, we first grew a large regression tree, where each split was induced by a threshold on some selected covariates. We then pruned the tree and selected the best-sized subtree based on an interaction-complexity measurement and an additional penalty on the total number of internal nodes (Su et al., 2009).

In presence of a large number of interaction terms especially when many categorical variables exist, we adopted two options in our method. For the first option, we used a screening step similar to the sure independence screening method (Fan and Lv, 2008) to ease the computational burden of solving (4). In particular, we used a group LASSO procedure (Yuan and Lin, 2005) to pre-screen the interaction terms at the variable level in (4). For example, for the interaction between one variable with C_j levels and one variable with C_k levels, we treated the total C_j × C_k levels of interactions as a group and imposed a group LASSO penalty on it. For the second option, we abandoned this interaction screening step and directly solved (4). We chose tuning parameters λ_1n and λ_2n by the Akaike Information Criterion (AIC). We labeled our method using the latter option as “ROWSi_–onestep” which allows us to investigate the advantage from the “screening step.”

All of our simulated data sets contained five binary covariates X_A, X_B, X_C, X_D, X_E, five ordinal covariates X_a, X_b, X_c, X_d, X_e and one continuous covariate X_Ca. Here we use capital letters in the subscript for binary variables, small letters for ordinal variables, and Ca for the continuous variable. For binary covariates, we use, for example X_A1, to denote the case level of X_A. For ordinal covariates, we use, for example X_a3 to denote the third level of X_a. All binary covariates were simulated from the Bernoulli distribution with success probability of 1/2. All ordinal covariates had four uniformly distributed levels. The continuous covariate was simulated from the standard normal distribution. The treatment T was independent of the aforementioned covariates and took values of −1 and 1 with equal probabilities. We examined three different sample sizes: n = 500, 1000, and 2000 for each simulated data set and ran 200 simulations in each setting.

We considered a variety of models for logit P(Y = 1) as follows.

1. 0.5X_C1 + 2(X_B1 + X_a3∗X_A1)∗T;

2. 0.5X_C1 + 2{X_B1 + X_a3∗(X_b2 + X_b3)}∗T

3. 0.05(−X_A1+X_B1)+{(X_a2 +X_a3)+(X_b2 +X_b3)∗X_Ca}∗T;

4. log log{(X_b3 + X_c3) + 5(X_a2 + X_a3+X_A1X_B1)∗T + 20}²;

5. (X_A1 + X_B1) + 2∗;

6. 0.5X_A1 + 0.5X_B1 + 2I{X_Ca < 5, X_a < 2}∗T

7. 0.5X_A1 + 0.5X_B1 + 2I{X_Ca < 5, X_Cb < 2}∗T

8. 0.5X_Ca + 0.5X_Cb + 2I{X_Ca < −2, X_Cb > 2}∗T

Models (A)–(C) are logistic models with various interaction structures. In particular, Model (A) has one main effect (X_C1), one first-order interactions (X_B1 ∗ T), and one second-order interaction (X_{a3 *}X_A1*T). Model (B) is similar to (A) but has different second-order interactions. Model (C) has more complicated structure of main effects and interaction effects. Model (D) is not a logistic model and can be used to compare the robustness of different methods in case of link function misspecification. Model (E) is a logistic model without covariate–treatment interactions. Models (F)–(H) take treelike forms which are not linear in terms of covariates.

We evaluated the performance of each method by four criteria: sensitivity, specificity, prediction accuracy, and the expected outcome $E^{\widehat{\mathcal{D}}}(Y)$ under the resulting rules $\widehat{\mathcal{D}}$ (a.k.a the value function, Zhao et al. (2012)). Sensitivity is the proportion of true non-zero interactions (including both first-order and second-order interactions) being estimated as non-zero. Specificity is the proportion of true zero interactions being estimated as zero. Prediction accuracy is the proportion of empirical treatment assignments, $\mathrm{sign}({\mathit{X}}^{\prime }\widehat{\mathit{\beta }})$, that are consistent to the “oracle” assignments sign(X′β^†). Finally, interaction term selection accuracy is the proportion of a particular interaction term being selected.

To evaluate the sensitivity and specificity, Figure 1 plots the product of sensitivity and specificity for each method. Our method appears to outperform the competitors notably. Figure 2 plots prediction accuracy. Again, our method seemingly surpasses all other methods in almost all settings except under Model (H), where the true assignment rule is not linear. We present comparison of the value function in the Web Appendix C.

We used Models (A) and (C) to illustrate the construction of CI for enhanced comparative treatment effects. Model (A) does not satisfy the property that P(X′β^† = 0) = 0 whereas Model (C) does. Hence, we used the m-out-of- n bootstrap for Model (A) and the regular bootstrap for Model (C). The bootstrap resampling was repeated for 1000 runs, from which a CI was obtained. We ran the above procedures for 1000 repetitions and report the mean length of CIs in Table 1. The length of CIs decreases as n increases. We also see that all CIs only include positive numbers, which indicates that treatment effects are enhanced in the subgroups identified by our method. In addition, we also report the coverage probabilities of the CIs for $E_{\mathit{X}}\{d(\mathit{X},1,\widehat{\mathit{\beta }})\}$and $E_{\mathit{X}}\{d(\mathit{X},-1,\widehat{\mathit{\beta }})\}$. Both of the two expectations were evaluated from a large independent testing data set with a sample size of 100,000. The performance seems satisfactory when the same size is 2000. There is some over-coverage when n = 1000. We noted that the performance depended on the choice of α in the m-out-of-n bootstrap. For this reason, we also investigated the sensitivity of this choice in Web Appendix D. Instead of a subjective choice, a double bootstrap procedure (Davison and Hinkley, 1997; Chakraborty et al., 2013) can be used for choosing α in a data-driven manner. Due to its computational intensity, we do not pursue here.

5. Data Analysis

In this section, we apply our proposed method to two real studies and compare it with the FindIt, VT and IT methods. The first is the National Supported Work (NSW) study (LaLonde, 1986) that appeared in Imai and Ratkovic (in press) and the FindIt package. The second is the mammography screening study discussed in the introduction section.

6. Conclusion and Discussion

In this article, we proposed a regularized outcome weighted method for subgroup identification. The method differs from many existing methods for subgroup identification in that it is based on a target function that directly relates to covariate-treatment interaction. Such separation from the covariate main effects is important as these main effects usually affect the estimation of interactions. Because the clinical outcomes are used as weights instead of as modeling targets, our method can be extended to deal with continuous, time-to-event, and possibly contaminated outcomes. Another advantage of our method is that it could be easily implemented by using standard packages. In fact, under the logistic surrogate loss ϕ (x) = log(1 + e^−x), if the regularization part of (4) is ignored, we can essentially solve

based on a logistic regression model for T weighted by outcomes. This can be conveniently fitted using standard packages such as the SAS LOGISTIC procedure by using subject-specific weights. For more general settings, we have presented some codes in the Web Appendix F. To keep it easier for the readers to implement, we have only presented codes that use LASSO as the penalty. We are developing an R package that will also incorporate the fused LASSO penalty.

We demonstrated the impressive performance of our method compared with other methods in various settings. Moreover, we not only identified the subgroups, but also carried out an additional inference procedure to assess how much benefit such a subgroup identification procedure could bring. Our approach for subgroup identification is mainly used to determine the rules for treatment assignment. This should be distinguished from other approaches for “subgroup identification” that focus on identifying some natural regions in the covariate space with enhanced treatment effects with no implication that the complement of that region constitutes another subgroup with a negative treatment effect (Kehl and Ulm, 2006; Lipkovich et al., 2011; Dusseldorp and Van Mechelen, 2014). This distinction is important as it may lead to different strategies in personalized medicine: one is to find the right drugs for patients and the other to find right patients for drugs.

For simplicity of presentation, we have assumed a complete randomization scheme from clinical trial data, that is P(T = 1) = π. Our method also works when the treatment assignment depends on X, that is, P(T = 1|X) = π(X) by replacing Tπ + (1 − T)/2 in (1)–(4) with Tπ(X) + (1–T)/2. The approach can also be extended to multiple interval treatments. The framework will be similar to dynamic treatment regimes, which are customized sequential decision rules for individual patients that can adapt over time to an evolving illness (Murphy, 2003; Zhao et al., in press). We can generalize the proposed method based on dynamic programming, which identifies the best rule at each interval backwards and recursively.

Zhang et al. (2012) showed that the optimal treatment regime could be regarded as the solution of the classification problem that ${\mathcal{D}}^{\ast }=\arg {min}_{\mathcal{D}}\mathrm{E}(|C(\mathit{X})|{[I\{C(\mathit{X})>0\}-\mathcal{D}(\mathit{X})]}^2)$, in which C (X) = E(Y|T = 1, X) − E(Y|T = −1, X) is the treatment contrast. Our formulation in (1) can be generalized by using C(X) based subject-specific weight. However, C(X) is unknown and needs to be estimated. Zhang et al. (2012) proposed a doubly robust estimator of C(X), which protects the estimation against model misspecification, especially when treatment assignment probabilities are known as in a randomized clinical trial. It will be of great interest to further investigate such a procedure as well as its theoretical properties.

Supplementary Material