Medicine Mathematics
Computer Science-2021.gz:143337
Computer Science
0fdc65062b08c181e1dc142135ad475ad6128e67
Robust Estimation of Heterogeneous Treatment Effects using Electronic Health Record Data

Estimation of heterogeneous treatment effects is an essential component of precision medicine. Model and algorithm-based methods have been developed within the causal inference framework to achieve valid estimation and inference. Existing methods such as the A-learner, R-learner, modified covariates method (with and without efficiency augmentation), inverse propensity score weighting, and augmented inverse propensity score weighting have been proposed mostly under the square error loss function. The performance of these methods in the presence of data irregularity and high dimensionality, such as that encountered in electronic health record (EHR) data analysis, has been less studied. In this research, we describe a general formulation that unifies many of the existing learners through a common score function. The new formulation allows the incorporation of least absolute deviation (LAD) regression and dimension reduction techniques to counter the challenges in EHR data analysis. We show that under a set of mild regularity conditions, the resultant estimator has an asymptotic normal distribution. Within this framework, we proposed two specific estimators for EHR analysis based on weighted LAD with penalties for sparsity and smoothness simultaneously. Our simulation studies show that the proposed methods are more robust to outliers under various circumstances. We use these methods to assess the blood pressure-lowering effects of two commonly used antihypertensive therapies.


Introduction
The ultimate goal of precision medicine is to optimize therapeutic outcomes by tailoring medical treatment and care provision according to individual patient characteristics. In practice, such tailoring must be guided by causal treatment effects expressed as functions of the observed patient characteristics x 1 , which account for patient heterogeneity in a given clinical population. But in reality, the true treatment effect function τ 0 (x) is almost never known and cannot be easily ascertained from clinical trials.
There is a sizable literature on the estimation of treatment effects in the form of τ 0 (x). With covariates averaged out, τ 0 (x) is reduced to the average treatment effect (ATE) τ 0 = τ 0 (x)f (x)dx, which can be estimated from clinical trials as well as observational studies 2 . While randomized experiments provide by far the most straightforward estimation of τ 0 , valid estimates can also be ascertained from observational data, by using the Neyman-Rubin causal model under appropriate assumptions 3 . Estimating treatment effect in the presence of heterogeneity, however, is a much involved task. Popular approaches include the advantage or A-learning methods that directly model the contrasts among treatments 4,5 , and the quality or Q-learners that regress the outcomes on patient characteristics 6,7 . Under the general umbrella of A-learners, Tian, Chen, and colleagues described a covariate-modification method 8,9 . More recently, Nie and Wager proposed a two-step learning algorithm that possesses a quasi-oracle property for estimating τ 0 (x) 10 . Xiao and colleagues further improved the algorithm for enhanced robustness 11 .
The performance of the above causal estimators is often influenced by the features of the observed data.
An attractive and readily available data source for causal inference is electronic health records (EHR), digitalized medical records collected and maintained by health care organizations 12 . While statisticians have long recognized the values of EHR data in causal analysis 13 , they are also keenly aware of the challenges presented by such data, including data outliers and high dimensionality. The former could result in biased estimation and questionable inference, whereas the latter leads to a "curse of dimensionality" 14 .
In this research, we address the above issues in a broader context of heterogeneous treatment effect estimation.
Specifically, we put forward a general estimation framework based on weighted score equations. The new formulation unifies many of the existing learners, while retaining the flexibility to accommodate different loss functions, permitting for example robust least absolute deviation (LAD) regression. The estimating formula enhances modified-covariate method's capacity against outliers 9 and extends the robust R-learner's ability to handle higher dimensionality 11 , giving each an improvement.
The approach's direct targeting of τ 0 (x) relates it nicely to the concept of the A-learning methods. We performed extensive simulation studies to investigate the new methods' operational performance, in comparison with the existing ones. We also described a real data application to illustrate the use of the proposed methods.
2 Proposed Methods

Models and Assumptions
We consider the estimation of τ 0 (x), the conditional average treatment effect (CATE), within the Neyman-Rubin potential outcome framework 15 . The binary treatment indicator T takes values 1 or −1, i.e., T ∈ {±1}. We let Y (1) and Y (−1) be the potential outcomes under T = 1 and T = −1, respectively.
We assume that data {(Y i , T i , X i )} n i=1 are independent and identically distributed (i.i.d.), where the pre-treatment covariates X i could be high dimensional as in EHR analyses. We require the stable unit treatment value assumption (SUTVA) 16,17 and write the observed outcome as Y = I(T = 1)Y (1) + I(T = −1)Y (−1) , where I(·) is an indicator function.
Within this framework, we focus on where the last part comes from the ignorability assumption defined below. This makes CATE estimation possible when X contains all confounders. When T ∈ {±1}, we can always express the conditional mean outcome as . This leads to a general interaction model where ε i is subject to Assumption 3 below, along with the other assumptions stipulated by Rubin and Rosenbaum 15,18 .
In the existing literature, τ 0 (x) is often depicted by a simple parametric model 8,11 . With µ(x) = Robinson decomposition used by the R-learner 10 .

Assumption 3 (Conditional Independence Error)
The error is independent of the treatment assignment conditional on covariates, i.e. {ε i ⊥ ⊥ T i |X i }. We further assume that the conditional expectation of error exists.

A Unified Formulation for Heterogeneous Treatment Effect Estimation
There are two general strategies for estimating τ 0 (x) in (2). The first is to depict the conditional mean function µ t (x) = E[Y |X = x, T = t] with a regression model and then obtain the treatment effect estimator τ (x) =μ 1 (x) −μ −1 (x). For example, from the objective function one can estimate β and γ simultaneously, and then achieves a CATE estimateτ (x) = τ (x;β) 19 . Such an approach is often referred to as the Q-learning, because its objective function plays a role similar to that of the Q or reward function in reinforcement learning 19 . The frequently used Two-or Single-learners (T or S-learners for short) are variants of this approach 20 .
An alternative strategy, one that we follow in the current research, is to directly target τ 0 (x) in a predefined objective function. This approach is often referred to as the A-learning 21 . A-learning first emerged in the context of dynamic treatment regime 4,5 , and was later generalized to one-stage case for treatment effect estimation 8,9 . In this paper, we show that there exists a unified formulation for the objective function, written in the form of score equations, that covers many of the existing learners.
Before introducing the general formulation, we first review the existing methods to highlight their connections.
1. The modified outcome methods. Certain transformations of Y could be used to facilitate the estimation of τ 0 (x). Estimation methods relying on such transformations are collectively known as the modified outcome methods. This class of methods includes the inverse propensity score weighting (IPW) 22,23 and the augmented IPW (AIPW) methods 24 . A common feature of this class of methods is to express the true treatment effect τ 0 (x) as a conditional expectation of the transformed outcome variables.
For IPW and AIPW, the transformations are Writing the modified outcome as Y * one has E(Y * |X) = τ 0 (X). An estimate can therefore be obtained by minimizing the square error loss, i.e., min τ (·) 2. The modified covariates methods. An alternative set of methods, collectively known as the modified covariates methods, have been derived from the model (2). The central idea of this approach is to estimate τ 0 (X) by re-weighting the loss function instead of the response variable 8,9 where the two weight functions w(x, t) and c(x, t) are subject to the following constraints for all x and t: One can show that the existing estimation methods, including MCM, MCM-EA, RL, IPW, and AIPW, are all covered by this general formulation. In Appendix A.1, we show that for each of the above methods, the corresponding functions c and w meet the three conditions.
A few additional remarks are in order for this general expression: Conditions C1-C3 are put in place to assure E(S|X) = 0. It can be shown that under the square error loss function, the estimates derived from (4) are indeed minimizers of the target function, i.e., in the same section of the appendix.
Remark 2. For given w(x, t) and c(x, t), one might be able to choose an appropriate g(x) to achieve robustness to model mis-specification. For example, the g(x) = (1 − p(x))µ 1 (x) + p(x)µ −1 (x) in the augmented inverse probability weighting (AIPW) method with equation leads to double robustness. Specifically, AIPW is robust against mis-specification of either propensity score model or both µ −1 (x) and µ 1 (x).

Remark 3.
When an additional condition c(x, 1) = 1 − p(x) holds and g(x) = µ(x), the score equation in (4) leads to an estimator with the minimized variance. For an R learner, we have c(x, 1) = 1 − p(x), and the choice of g(x) = µ(x) leads to the most efficient estimator. For MCM, this additional condition also holds when p(x) = 1 2 , as in the case of randomized clinical trials.

Remark 4.
With the unified formulation for the score functions, new estimators can be derived, for is an arbitrary augmented function of X.
With the score function expressed as in (4), we propose an estimation procedure for CATE τ (·), where ρ(·) is a user-specified loss function, and Λ n (·) is a structural penalty function for τ (·). This general procedure covers most of the existing methods for heterogeneous treatment effect estimation through a unified formulation.

Estimation Methods under the L 1 Loss
The estimation procedure described in (5) is general and flexible in the sense that it allows the analyst: (1) to choose different estimators through the specification of w(·) and c(·); (2) to select g(·) for efficiency enhancement; and (3) to specify a loss function ρ(·) that is most appropriate for the application. This general formulation provides a natural remedy to two practical issues in EHR data analysis: (1) lack of robustness of the L 2 -based methods against outliers, (2) lack of accommodation of the high dimensionality of X, and nonlinearity of τ (X).
Specifically, we put forward a class of robust estimators within the confines of the general estimating function (5). The method accommodates nonlinearity in τ (·), further enhancing the modeling flexibility.
Estimation is implemented under the usual causal inference Assumptions 1 − 3.
Under the L 1 -loss function, we show in Appendix A that with Conditions C1-C3, we have To increase efficiency, we opt to use g(X i ) = µ(X i ) in proposed methods.
Herein, we consider the following penalized least absolute deviation estimator where Λ n is added to ensure sparsity at the function level. For simultaneous variable selection and smooth estimation, we adopt a similar penalty term in (6) as described by Meier 28 .
We further assume an additive structure for the treatment effect function τ (·): where α is the intercept, and m j (·) is the jth additive component corresponding to x j . We write are the B-spline basis functions, K n and q are number of knots and degree.
Rewriting the spline bases and coefficients as vectors, we have τ ( (β 11 , β 12 , · · · , β 1(K+q) , · · · , β p(K+q) ) T . For simplicity, we choose a common K n + q for all spline components. Following a suggestion of K n √ n + 4 by Meier 28 , we use K n = √ n/2, which is of the same order and not too large for implementation.
With this, we define the penalty term in (6) as where ||m j || 2 n = 1 n n i=1 m 2 j (X i ) = 1 n β T j D j β j is for variable selection in a group-wise manner, and and Ω j respectively. And P λ1,γ (·) is the smoothly clipped absolute deviation (SCAD) penalty defined by its first derivative with γ > 2 and P λ1,γ (0) = 0. We use γ = 3.7 as suggested by Yuan and Lin 29 .
Hence, optimization of (6) can be expressed as a general group SCAD problem With these transformations, the optimization of (6) becomes an ordinary least absolute deviation (LAD) regression with a group SCAD penalty where ||β j || is the Euclidean norm, The estimation of CATE is thereforeτ (x) =α +β TB (x).

A Computational Algorithm
To optimize (9), one has to estimate µ(·) and p(·), as they are involved in the weight functions w(x, t) and c(x, t). Herein, we use pre-estimatedμ(·) andp(·) as plug-in estimates for solving (9). Estimation accuracy of these quantities, however, can be impeded by the dimension of x i and the uncertainty of the functional forms of the x i 's associations with T i and Y i . To remedy, we use a gradient boosting machine (GBM) 30 to estimate these two functions , with packages gbm 31 and caret 32 . In cases of ultra-high dimensional x i , one could first use non-parametric independence screening (NIS) method 33 to reduce the dimensionality to a moderate one (n − 1 or log(n)) as suggested by Fan and Lv 34 ), before applying our proposed method.
With the plug-in estimates of µ(·) and p(·), we solve the L 1 optimization problem in (9), by using R package rqPen 35 , which is designed for penalized quantile regression in general. The nonconvex group penalized optimization with quantile loss is solved by the extension of quantile iterative coordinate descent (QICD) algorithm proposed by Peng and Wang 36 . For comparison purposes, we also use R package oem 37 to ascertain the L 2 estimators.
The main steps of the procedure are described in Algorithm 1.

Algorithm 1:
Input: Outcome Y , treatment assignment T , and pre-treatment covariates X 1 Data screening: Screen covariates with NIS when in situations of ultra-high dimension.
2 Nuisance quantity estimation: Estimate p(x) by using GBM with cross-validation (CV) and estimate µ(x) by using L 1 -based GBM with CV. 3 Data transformation: Construct the B-spline design matrix B(X), calculate w(X i , T i ), c(X i , T i ), and g(X i ) following Conditions C1-C3, and transform B(X) toB(X) using (8).
4 Optimization: Solve penalized LAD regression (9) with a group SCAD penalty to achieve estimates of α andβ with regularization parameters selected by CV.

Asymptotic Properties ofτ (x)
For theoretical examination, we consider the simple case of a univariate covariate X i ∈ R: where ρ(·) is a loss function that is convex and has unique minimizer at origin. This simplification will not diminish the contribution of the asymptotic analysis, which is complicated by the B-spline approximation and the various loss functions including the L 1 , L 2 , Huber, and Bisquare loss functions.
With a B-spline approximation, we write τ (x) := where q is the degree of the B-splines and K n is the number of knots, which we assume depending on sample size n. Zhou et.
al. 38 provided the L ∞ approximation error for B-splines. In particular, with τ * (x) := B(x) T β * as the best L ∞ approximation to the true function τ 0 (x), it satisfies with {κ k } Kn k=0 are the knots in the B-spline approximation, τ (q+1) 0 (x) is the (q + 1)th order derivative of τ 0 (x), and Br q (x) is the q-th Bernoulli polynomial.
We focus on the asymptotic theory of the L 1 spline estimatorτ (x) = B(x) Tβ , wherê The error forτ (x) can be decomposed as a summation of the estimation error and approximation error We only need to study the estimation errorτ (x)−τ * (x) = B(x) T (β −β * ) thanks to the L ∞ approximation result by Zhou et al. 38 .
To show a pointwise asymptotic normality of √ a n (τ (x) − τ * (x)) with a convergence rate a n to be specified later in Appendix A, we only need to prove the convergence of Then the minimizerδ n of U n (δ) is simply our target, i.e., δ n = √ a n (β − β * ).
If one regards {U n (δ)} as a sequence of random functions and the finite-dimensional distributions of U n (δ) converge in distribution to those of some random function U (δ) which has a unique minimum, then it will follow thatδ n = √ a n (β − β * ) → d argmin(U (δ)), as n → ∞ per Hjort and Pollard 39 , and Geyer 40 .
With a given loss function ρ(·), we define Φ(s|X = x, Let Φ (s|X = x, T = t) and Φ (s|X = x, T = t) be the first and second derivative of Φ(s|X = x, T = t) with forδ n ; respect to s. Several additional conditions are required for the proof of asymptotic normality: C4. X is distributed as Q(x) on a compact set in R. Without loss of generality, we assume X ∈ [0, 1].
C7. The function ρ(u) is convex, it has a unique minimizer at zero, and its first and second derivatives exist. To describe the asymptotic normality of the spline estimatorτ (x), we introduce two matrices: We and another square matrix D of the same dimension with its (i, j)-th element being .
Remark 6. With the order of K n larger than O(n 1 2q+3 ), the B-spline approximation error b a (x) can be ignored relative to the order of its variance.
For the rest of the paper, we focus on the LAD loss where Conditions C7-C10 are naturally satisfied, and C11 can be simplified as the following: C12. There exists a constant γ ≥ 0 such that E |w(X, T )c(X, T )| 2+γ |X = x < ∞.
To describe the asymptotic normality of the spline estimatorτ (x) under the L 1 loss, we write matrix D with the (i, j)-th element being where f 1 (y|x) and f −1 (y|x) are the conditional density functions of Y (1) and Y (−1) given X = x, respectively. We give the following theorem for the spline-based LAD regression: Remark 7. For inference concerning τ 0 (x), the variance of the estimator can be obtained by using resampling methods, as the asymptotic variance is difficult to work with. In a simulation experiment in Appendix B.2.2, we show that the bootstrap C.I. consistent with theoretical C.I.. 13 We conducted an extensive simulation study to evaluate the finite sample performance of the proposed methods. We considered a large number of parameter settings, including four different learners under two different loss functions: (1) a robust version of the modified covariate method with efficiency augmentation We used the A-learner described by Lu and colleagues 26 . The objective functions of the A-learning methods 3 and 6 shared the same structure, except for the loss function ρ where γ and p(x) are estimated in advance. We estimated γ by regressing Y on X using a linear regression, and p(x) by regressing T +1 2 on X using GBM. The objective functions of the Q-learning methods 7 and 8 shared the same structure where we used L 1 or L 2 loss function for ρ. Note that a difference between the A-learner and R-learner is the choice of the augmentation. For A-learner we used a linear function as suggested by Lu 26 to estimate We designed the simulation study to assess the robustness of the L 1 and L 2 -based methods, and to contrast the performance of the A and Q-learners. We also examined the performance of the methods under different sample sizes, dimensionality, and proportions of outliers.
We assessed the performance of the methods using the standard metrics, including bias, variance, mean square error as well as mean absolute error. In addition, we compared the value function Q(η) = E(Y (η)), i.e., the expected average outcome under treatmentη, whereη(x) = 2I(τ (x) > 0) − 1, as recommended by each method 41 . To estimate the Q(η) for a given regimen, we conducted a Monte Carlo simulation using model Y (η) = b 0 (X) +η 2 τ 0 (X) + ε, replacing T in (2) byη, and we set the number of replicates is 10 6 . The value function calculated based on the true treatment effects was E[Y (η opt )] = 1.25, where We also assessed the sensitivity and specificity for variable selection under our penalty. With the number of simulation replication R, we defined

Data Generation
We generated data as follows, the dimension of the covariates was indexed by p:  was not a linear function. And there were little practical differences between the robust R-learner and robust MCM-EA. The Q-learner performed the best under MSE and MAE because it is a one-step estimation procedure, and thus avoiding the errors associated with the nuisance quantity estimation. This is consistent with the observations made by Schulte 21 that the Q-learner tended to perform better than the standard A-learner when all models were correctly specified. We conducted a separate simulation for a setting where the Q-function was mis-specified. The results reported in Appendix B.2.3 showed that in the presence of outliers, bias in the mis-specified L 1 -QL was larger than that of the L 1 -MCMEA, L 1 -RL, and L 1 -AL. The same was also true for MSE and MAE. In terms of the value function Q(η), L 1 -QL had smaller Q(η) values than methods under the A-learning umbrella; findings were consistent with MSE.  Table B.5.

Simulation Results
To illustrate the methods we propose, we estimated the treatment effects of two different antihypertensive therapies by analyzing the observed clinical data set from the Indiana Network of Patient Care, a local EHR system. The data were a subset of a previous study assessing the blood pressure (BP)-lowering effects of various antihypertensive agents 42 . This analysis compared the BP effects of angiotensin-convertingenzyme inhibitors (ACEI) alone and a combination of ACEI and hydrochlorothiazide (HCTZ). We considered those on ACEI alone as in treatment group A, and those on ACEI+HCTZ as in group B. The primary outcome of interest is clinically recorded systolic BP in response to these therapies. Independent variables included the demographic and clinical characteristics, as well as medication-use behaviors of the study participants. Data from 882 participants were used in the current analysis. Among these, 350 were on the monotherapy of ACEI, and 532 were on the combination therapy of ACEI+HCTZ. Characteristics of the study participants are presented in Table 1. There were four continuous variables (pulse, BMI, age, and medication adherence) and 12 binary variables (gender, race, and ten comorbidities). The continuous variables were standardized before the analysis and expressed as linear combinations of splines.
We expressed the treatment effect of treatment B, in comparison against treatment A, as a function of the patient characteristics x where Y (A) , Y (B) represented the potential systolic BP of ACEI alone group and ACEI+HCTA group.
Since the antihypertensive effect of a therapy is measured by its ability to lower BP, a negative τ (x) indicates a superior effect of the combination therapy over the monotherapy, for a given x. An important covariate of interest was the level of medication adherence, which we measured with the proportions of days covered (PDC) by the medication.
Preliminary data examination showed that the observed systolic BP was right-skewed in both groups.
The Shapiro-Wilk's test further confirmed that the systolic BP was not normally distributed, and there were outliers in the observed outcome (ACEI alone: W = 0.9912, p = 0.035; ACEI+HCTZ: W = 0.9617, p = 1.498e − 10). We, therefore, used the L 1 -based methods with additive B-splines to analyze the data. Here the B-splines were used to accommodate the possible nonlinear influences of the independent variables on the treatment effect.
Naive comparison of the systolic BP-effects between the two treatment strategies suggested that the combination therapy (ACEI+HCTZ) was significantly worse than the monotherapy (ACEI alone) in its ability to lower systolic BP ( BP is regulated by hormones in the renin-angiotensin-aldosterone system (RAAS) 44 . ACE inhibitors block the conversion of angiotensin I to angiotensin II, diminishing the latter's effects on aldosterone production and sodium retention and causing BP reduction. Thiazide diuretics lower BP by suppressing the extracellular fluid volume, which in turn reduces aldosterone secretion. Together, the two drugs are expected to have additive effects in lowering BP. In clinical practice, the two are often used concurrently.
( Table 1 goes here) A closer examination of the characteristics of the patients on these therapies showed that patients on the combination therapy were older, more likely to be female, and overweight. Using GBM described in Section 2.3, we examined the mean function of systolic BPμ(x) and the propensity of a patient receiving the combination therapyp(x). The estimated propensity score distributions were clearly different for the two treatment groups, whereas the mean functions were similar. See Appendix C. More specifically, the histogram of mean functions overlapped, indicating no apparent differences between the mean systolic BP between the two treatment groups. The different propensity score distributions of the two groups clearly showed that non-random treatment assignment. The importance levels of the covariates from GBM and additional modeling details were summarized in Appendix C. The systematic differences in patient characteristics between the two treatment groups suggested that a naive comparison was not appropriate and should not be trusted.
We then analyzed the data with the proposed methods. Importantly, both the L 1 -MCM-EA and  Lowerτ indicated a stronger efficacy of the combination therapy than the monotherapy. Although decreasing trends were observed in both L 1 and L 2 -based methods, the L 2 methods failed to detect any differences between the two therapies, as the 95% confidence intervals forτ (P DC) consistently covered zero. The L 1 -based estimators, however, showed a superior blood pressure-lowering effect of the combination therapy, but only when PDC> 90%. The fact that treatment effects varied with medication adherence should not be surprising. As the former US Surgeon General, Dr. C. Everett Koop, wisely observed, "Drugs don't work in patients who don't take them." 45 In this analysis, we do not expect significant differences between the treatments when patients are not adherent to the prescribed regimen.
Findings such as this are not unexpected in comparative effectiveness analysis of EHR data. Because unlike well-controlled clinical trials, few measures are in place to ensure patients faithfully take their medications in the real-world of clinical care. In the current application, the fact that the L 1 -based estimators detected significant differences highlights the proposed methods' advantage. Using L 1 -based estimators, we also examined the influences of BMI on τ , which did not reach the level of statistical significance (data not shown).
( Figure 3 goes here) To check the conditional independence error assumption, we performed the invariant residual distri-  46 , and invariant residual prediction test (IRP-test) 47 . The conditional independence error assumption held for both proposed methods at the significant level of 0.05 (see Table 2).
The estimates based on L 2 -MCMEA and L 2 -RL were slightly higher than the corresponding L 1 -based methods. This results in Table 3 showed that the SBP could be reduced if treatment were to be assigned in accordance with the therapy recommended by the estimated treatment regime. Table 4 showed among the patients included in the analysis, 100 (11.3%) had PDC above 90%. We further examined the numbers of patients assigned to the two different treatment groups based on the estimated treatment effects.
More patients would be assigned to the combination therapy group because it had a significantly greater blood pressure efficacy when patients take their medications. On the other hand, had we used the L 2 based methods, almost all of the patients would have been assigned to the monotherapy group, which contradicts the recommendations from the current clinical guidelines.
( In summary, the naive and L 2 -based methods showed that the combination therapy of ACEI and HCTZ had a worse BP-lowering effect than the monotherapy of ACEI, a finding that contradicts the recommendations of the current clinical guidelines of hypertension treatment. The L 1 -based methods have produced results that are better explained by the existing clinical and biological evidence. The analysis showed that treatment effects tended to improve when patients adhere to their prescribed medications.

Discussion
We started this work searching for a robust estimator for heterogeneous treatment effects that could be used in EHR analysis, where outliers often undermine the validity of estimation. In the process, we       Note: Patients in the application data set whose P DC > 0.9 are reassigned treatments by estimated treatment effect, i.e.η(x) = 2I(τ (x) < 0) − 1. Under the L1-based methods, most of the patients will be assigned to the combination therapy group, consistent with the results in Figure 3. However, under the L2-based methods, most of the patients will be assigned to the monotherapy group, which is counter-intuitive, because when patient adhere to prescription, the combination therapy is known to be more efficacious.    (1) MCM-EA. The objective function of L 2 -MCM-EA method is We write Then which shows the c and w functions satisfy Conditions C1 and C2. Condition C3 (w > 0 and c = 0) is clearly met. The same set of parameters can be used in L 1 loss. The verification is the same.
(2) R-Learning. The objective function of L 2 -based R-learning method is We write Then Therefore, Conditions C1-C3 are met. The same specification works for L 1 loss. The verification of A-learning remains the same.
We write Therefore, Conditions C1-C3 are met. The same specification works for the L 1 loss function.
(4) AIPW. The verification of AIPW method is the same.

A.2. Basic properties of the general formulation
Proof of Property 1.
Proof of Property 2. We provide the optimal g(·) in this section, by optimality we mean the g(·) that minimizes the variance of estimator. Let S(Y i , X i , T i ; τ (X i )) be the derivative of the objective function 2 , with respect to τ . Then the estimating equation is , the solution of the augmented score equation always converges to τ 0 (·) in probability. Following 8,9 , selecting the optimal g(·) is equivalent to minimizing the conditional variance of , which satisfies the equation for any function η(·). By interaction model (1) and Condition C2, the expression of g 0 (x) can be further simplified to g 0 (x) = µ(x) + [1 − p(x) − c(x, 1)]τ 0 (x). As τ 0 (·) is the unknown target, when c(x, 1) = 1 − p(x), the optimal augmentation function is mean outcome function, i.e., g 0 (x) = µ(x).
Proof of Property 3.

A.3. Asymptotic Properties
To prove Theorem 1, we first introduce two lemmas.

Lemma 1.
Under the same assumptions as Theorem 1, W n is asymptotically equivalent to the (K + p)dimensional normal with mean 0 and variance G.
with respect to X i is as follows. First, we calculate the conditional expectation with respect of X i and T i , From (1) to (2) is based on Condition C1. Then, based on the interaction model and the distance between whereõ (1) uniformly holds for all x by the distance between τ 0 (x) and B(x) T β * . From (3) to (4) is where α (1) , α (−1) ∈ (0, 1) are from Taylor expansion. The first two terms in (5) =o(1).
As the order of b a (x) is o(K −(q+1) n the conditional expectation is of o(1). This means the conditional expectation of the score function with loss functions satisfy Conditions C7-C11 goes to zero.
where the last two steps are derived by Condition C11. The conditional variance of Z n respect to X (n) can be calculated as following.
where G is the variance of W n . Here the derivation from (7) to (8) So it follows that from Lyapunov's Theorem. The conditional expectation of Z n respect to X (n) can be calculated as The last step is from the proof of Lemma 6.10 of 48 and equation (6), for j = −p + 1, ..., K n , we have Proof of Lemma 2. The (i, j)-component of D is From the fundamental property of B-spline function (Lemma 6.1 in 38 ), we have From the property of B-spline function 49 , D is positive definite matrix. Therefore, Now, we are ready to prove Theorem 1. For simplicity we write a n as ∼ b n , where random sequence {a n } and {b n }, if a n /b n = O P (1).
Proof of Theorem 1. The objective function of proposed method is As this minimization problem doesn't have explicit solution, for the convergence of √ a n (β − β * ), we modify the objective function L n (β) as follows: Then the minimizerδ n of U n (δ) can be obtained aŝ where X (n) represents all the observed X. We have E[R n (δ)|X (n) , T n ] = 0 from the straight calculation. Let Then by Condition C9-10 with s = α n c( Therefore, we have from K n = o(n 1/2 ), E[R n (δ) 2 ] = 1 n V [r 1 ] = o(1) and R n (δ) = o P (1). By the definition of Φ(t|X, T ), the Taylor expansion of Therefore, the conditional expectation of U n (δ) given X n can be written as Thus, we have U n (δ) as The minimizer of U n (δ) isδ which is the solution of ∂Q n (δ)/∂δ = 0. Hence, becauseδ = 1 αn (β − β * ), we have The asymptotic variance ofτ (x) is similar to that ofτ (x) − τ * (x) because W n is the only random vector in the asymptotic form ofτ (x), it is easy to show that due to the Riemann integral, fundamental asymptotic property of B-spline basis, and Lemma 2. Under the condition K n = O(n 1/(2q+3) ), we . This completes the proof.
To prove Theorem 2, we introduce two lemmas as well.
. Under the same assumptions as respect to X i can be calculated as following.
The derivation from the first equation to the second equation is by Condition C1, that from the third equation to the fourth equation is by Model (1), that to the fifth equation is by proposed Condition C2, the last two steps are by Taylor expansion and the order of b a (x i ).
where the last two steps are derived by Condition C12. The conditional variance of Z n respect to X (n) can be calculated as following.
So it follows that since γ ≥ 0. This leads to from Lyapunov's Theorem. The conditional expectation of Z n respect to X (n) can be calculated as (1)).
Therefore we obtain For i = 1, ..., n, we have Therefore the variance of R n can be evaluated as and hence, Lemma 4 holds.
Now we are ready to prove Theorem 2.
Proof of Theorem 2. The objective function of proposed method is Then the minimizerδ n of U n (δ) can be obtained asδ n = … n K n (β − β * ).
Following the Knight's identity, we can write U n (δ) as Therefore, for both methods we obtain Because U 0n (δ) is convex with respect to δ and has unique minimizer, the minimizer of U n (δ) converges to δ 0 = argmin δ {U 0n (δ)}. This fact is detailed in 50 . Hence, we have Sinceτ (x) − τ * (x) = B(x) T (β − β * ), we obtain for x ∈ (0, 1), as n → ∞, Under the condition K n = O(n 1/(2q+3) ), we have . This completes the proof. We designed and conducted simulation on an additional parameter setting, to assess the performance of the proposed methods in high-dimension situations. As described in Section 2.3, we used NIS to screen the covariates in the first step.

B.1. Simulation Results of Settings 1-3
We generated data as follows. the dimension of the covariates was indexed by p: where p o is the proportion of outliers, n = 1000, p o ∈ {0, 0.05}, and p ∈ {1000, 3000, 5000}

B.2.2. Effects of Smoothness Penalty
Finally, we investigated the effects of an added smoothness penalty. We considered a situation involving a univariate covariate x. We visualized the performance differences of the L 1 and L 2 methods, with and without the smoothness penalty.

B.2.3. Comparison of Q-learning and A-learning when model is mis-specified
We investigated the performance of Q-learning and proposed methods when model is mis-specified.
We summarize the MSE and MAE of the Q-learning and proposed methods with combination of L 1 and L 2 loss when there is a small amount of outliers.

C.1. Existence of outliers
The following figure shows that the outcome observations in both treatment groups are beyond normally distributed.  In the application of proposed method, the importance levels from GBM are consistent with the result from regression. The importance levels from GBM and the linear and logistic regression results are summarized in the following two tables.