---
author:
- |
    Wei Xiao$^\dagger$, Hao Helen Zhang$^\ddagger$, and Wenbin Lu$^\dagger$\
    \
    \
bibliography:
- 'OTR.bib'
title: Robust regression for optimal individualized treatment rules
---

> [*Abstract:*]{} Because different patients may response quite differently to the same drug or treatment, there is increasing interest in discovering individualized treatment rule. In particular, people are eager to find the optimal individualized treatment rules, which if followed by the whole patient population would lead to the “best” outcome. In this paper, we propose new estimators based on robust regression with general loss functions to estimate the optimal individualized treatment rules. The new estimators possess the following nice properties: first, they are robust against skewed, heterogeneous, heavy-tailed errors or outliers; second, they are robust against misspecification of the baseline function; third, under certain situations, the new estimator coupled with pinball loss approximately maximizes the outcome’s conditional quantile instead of conditional mean, which leads to a different optimal individualized treatment rule comparing with traditional Q- and A-learning. Consistency and asymptotic normality of the proposed estimators are established. Their empirical performance is demonstrated via extensive simulation studies and an analysis of an AIDS data.
>
> [*Key words and phrases:*]{} Optimal individualized treatment rules; Personalized medicine; Quantile regression; Robust regression.

Introduction
============

Given the same drug or treatment, different patients may respond quite differently. Factors causing individual variability in drug response are multi-fold and complex. This has raised increasing interests of individualized medicine, where customized medicine or treatment is recommended to each individual according to his/her characteristics, including genetic, physiological, demographic, environmental, and other clinical information. The rule that applied in personalized medicine to match each patient with a target treatment is called individualized treatment rule (ITR), and our goal is to find the “optimal” one, which if followed by the whole patient population would lead to the “best” outcome. For many complex diseases such as cancer and AIDS, the optimal individualized treatment rule or regime is a dynamical treatment process, involving a sequence of treatment decisions made at different time points throughout the disease evolving course.

Q-learning [@watkins1992q; @murphy2005generalization] and A-learning [@murphy2003optimal; @robins2004optimal] are two main approaches for finding optimal dynamic individualized treatment rules based on clinical trials or observational data. Q-learning is based on posing a regression model to estimate the conditional expectation of the outcome at each time point, and then applying a backward recursive procedure to fit the model. A-learning, on the other hand, only requires modeling the contrast function of the treatments at each time point, is therefore more flexible and robust to a model misspecification. See [@schulte2014q] for a complete review and comparison of these two methods under various scenarios, in terms of the parameter estimation accuracy and the estimation of expected outcomes. Q- and A-learning have good performance when model is correctly specified but are sensitive to model misspecification. To overcome this shortcoming, several “direct” methods have been proposed, which maximize value functions directly instead of modeling the conditional mean. See [@ZhaoYingQi2012OWL; @Zhang2013Robust] for example.

All existing methods for optimal individualized treatment rule estimation, including Q-learning and A-learning, belong to mean regression as they estimate the optimal estimator by maximizing expected outcomes. In the case of single decision point, Q-learning is equivalent to the least-squares regression. Least-squares estimates are optimal if the errors are i.i.d. normal random variables. However, skewed, heavy-tailed, heteroscedastic errors or outliers of the response are frequently encountered. In such situations, the efficiency of the least square estimates is impaired. One extreme example is that when the response takes i.i.d. Cauchy errors, neither Q-learning nor A-learning can consistently estimate the optimal ITR. For example, in AIDS Clinical Trials Group Protocol 175 (ACTG175) data [@hammer1996trial], HIV-infected subjects were randomized to four regimes with equal probabilities, and our objective is to find the optimal ITR for each patient based on their age, weight, race, gender and some other baseline measurements. The response CD4 count of the data follows a skewed, heteroscedastic errors, which weakens the efficiency of classical Q- and A-learning. A method to estimate optimal ITR which is robust against skewed, heavy-tailed, heteroscedastic errors or outliers is highly valuable. One possible solution is to construct the optimal decision rule based on the conditional median or quantiles of response given covariates than based on average effects.

In the following, we present a simple example where a quantile-based decision rule is more preferable than a mean-based decision rules. We use higher value of response $Y$ to indicate more favorable outcomes. Figure \[fig:plot1\] plots the conditional density of $Y$ under two treatments, $A$ and $B$, given a binary covariate $X$ which takes the value of male and female. Under the comparison based on conditional means, $A$ and $B$ are exactly equivalent. However, conditional quantiles provide us more insight. For the male group, the conditional distribution of response given treatment $B$ is a log-normal and skewed to the right. Therefore, treatment $B$ is less favorable when either 50% or 25% conditional quantile are considered. For the female group, the conditional distribution of response given treatment $A$ is a standard normal while a Cauchy distribution given treatment $B$. Therefore, if we make a comparison based on $25\%$ conditional quantile, treatment $A$ is more favorable.

![The distribution functions of the response $Y$, in a randomized clinical trial with two treatments, $A$ and $B$, for male (two panels on the left) and female (two panels on the right). The solid lines with triangle symbol, dashed line, and dotted lines are the conditional mean, $50\%$ quantile, and $25\%$ quantile functions of $Y$ given the gender and the treatment, respectively.[]{data-label="fig:plot1"}](plot1.eps){width="6in"}

In this paper, we propose a general framework for optimal individualized treatment rule estimation based on robust regression, including quantile regression and the regression based on Huber’s loss and $\epsilon$-insensitive loss. The proposed methodology has the following desired features. First, the new decision rule obtained by maximizing the conditional quantile, which is suitable for skewed, heavy-tailed errors or outliers. Second, the proposed estimator requires only modeling the contrast function between two treatments, and is therefore robust against misspecification of the baseline function. This property is shared by A-learning. Third, empirical results from our comprehensive numerical study suggest favorable performance of the new robust regression estimator.

The rest of the paper is organized as follows. In Section 2, we first review the classical Q- and A- learning methods. Then we propose the new procedure and method and discuss its connection with existing methods. In Section 3, we study and prove the asymptotic properties of the proposed method, including consistency and asymptotic normality. In Section 4, a comprehensive numerical study is conducted to assess finite sample performance of the new procedure. In Section 5, we apply the method to ACTG175 data. Concluding remarks are given in Section 6. Throughout the paper, we use upper case letters to denote random variables and lower case letters to denote their values.

New Optimal Treatment Estimation Framework: Robust Regression
=============================================================

Basic Notations and Assumptions
-------------------------------

For simplicity, we consider a single stage randomized clinical trial with two treatments. For each patient, the observed data is $({{\bm X}},A,Y)$, where ${{\bm X}}\in\mathcal{X}={\mathrm{I \! R} \mathit{^{p}}}$ denotes the baseline covariates, $A\in\mathcal{A}=\{0,1\}$ denotes the treatment assigned to the patient, and $Y$ is the real-valued response, which is coded so that higher values indicate more favorable clinical outcomes. An ITR $g$ is a function mapping from $\mathcal{X}$ to $\mathcal{A}$.

We first review the potential outcome framework [@neyman1923applications; @rubin1974estimating; @rubin1986comment]. The potential outcome $Y^*(a)$ is the outcome for an arbitrary individual has s/he received treatment $a$. In actuality, at most one of the potential outcomes can be observed for any individual. The optimal ITR under mean regression, which maximizes the expected outcome, is $g^{\mathrm{opt}}_\mu={\mathrm{argmax}}_{g\in\mathcal{G}}{\mbox{E}}[Y^*\{g({{\bm X}})\}]$. Define the propensity score $\pi({{\bm X}})\triangleq P(A=1|{{\bm X}})$. Following [@rubin1974estimating] and [@rubin1986comment], we can compute the expectation of the potential outcome under the following two key assumptions.

-   **Stable Unit Treatment Value Assumption (SUTVA):** a patient’s observed outcome is the same as the potential outcome for the treatment that s/he actually received. Based on [@rubin1986comment], the SUTVA assumption implies that the value of the potential outcome for a subject does not depend on what treatments other subject receive. Specifically, we can write the SUTVA assumption as $$Y_i=Y_i^*(1)A_i+Y_i^*(0)(1-A_i),\;i=1,\ldots,n.$$ This is also referred as consistency assumption.

-   **Strong Ignorability Assumption:** the treatment assignment $A$ for an individual is independent of the potential outcomes conditional on the covariates ${{\bm X}}$, i.e., $A\bot \{Y^*(a)\}_{a\in\mathcal{A}}|{{\bm X}}$. For a randomized clinical trial, this assumption is satisfied automatically. For an observational study, as clinicians make decisions based only on all past available information, this assumption essentially assumes no unmeasured confounders.

For consistent estimation of the optimal treatment rule, we also need to assume

-   **Positivity Assumption:** $0<\pi({{\bm x}})<1$, $\forall {{\bm x}}\in\mathcal{X}$.

Existing Learning Methods: Q-learning and A-learning
----------------------------------------------------

Define the Q-function $Q({{\bm x}},a)\triangleq{\mbox{E}}(Y|{{\bm x}},a)$. Under assumptions (C1)-(C2), one can show that $g^{\mathrm{opt}}_\mu({{\bm X}})={\mathrm{argmax}}_{a\in\mathcal{A}}Q({{\bm x}},a)={\mathrm{argmax}}_{a\in\mathcal{A}}{\mbox{E}}(Y|{{\bm X}},A=a)$. This suggests that, in order to find $g^{\mathrm{opt}}_\mu$, we only need to estimate the conditional expectation of $Y$ given $({{\bm X}},A)$. This result serves as the foundation of Q- and A-learning framework. We further define the value function $V_\mu(g)={\mbox{E}}_{{{\bm X}}}[Q\{{{\bm X}},g({{\bm X}})\}]$ which is simply the marginal mean outcome under the ITR $g$, and $g^{\mathrm{opt}}_\mu={\mathrm{argmax}}_gV_\mu(g)$.

Define the $\tau$-th conditional quantile of $Y$ given $({{\bm X}},A)$ as $Q_{\tau}({{\bm X}},A)\triangleq\inf\{y: F_{Y|{{\bm X}},A}(y)\geq\tau\}$. Then we define the value function based on the $\tau$-th conditional quantile as $V_{\tau-q}(g)={\mbox{E}}_{{{\bm X}}}[Q_{\tau}\{{{\bm X}},g({{\bm X}})\}]$, which is an analog to the definition of $V_\mu(g)$. The optimal ITR which maximizes the $\tau$-th conditional quantile is then defined as $$g^{\mathrm{opt}}_{\tau}({{\bm x}})=\underset{a\in\mathcal{A}}{\mathrm{argmax}}Q_{\tau}({{\bm x}},a),\;\tau\in[0,1],$$ and $g^{\mathrm{opt}}_\tau={\mathrm{argmax}}_gV_{\tau-q}(g)$.

Consider the general model ${\mbox{E}}(Y|{{\bm X}},A)=h_0({{\bm X}})+AC_0({{\bm X}})$, where $h_0({{\bm X}})$ represents the baseline effect, and $C_0({{\bm X}})$ denotes the contrast effect as $$C_0({{\bm X}})={\mbox{E}}(Y|{{\bm X}},A=1)-{\mbox{E}}(Y|{{\bm X}},A=0).$$ Therefore, $g^{\mathrm{opt}}_\mu({{\bm X}})={\mathrm{1}}\{C_0({{\bm X}})>0$}. In Q-learning, a parametric model is often employed as a working model, $${\mbox{E}}(Y|{{\bm X}},A)=h({{\bm X}}; {{\mbox{\boldmath $\gamma$}}})+AC({{\bm X}};{{\mbox{\boldmath $\beta$}}}),$$ where $h({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ are posited parametric models for $h_0({{\bm X}})$ and $C_0({{\bm X}})$ respectively. Commonly a linear model is assumed for simplicity and interpretability, i.e., $h({{\bm X}};{{\mbox{\boldmath $\gamma$}}})={{\mbox{\boldmath $\gamma$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})={{\mbox{\boldmath $\beta$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$, where $\tilde{{{\bm X}}}=({{\mathbf 1}},{{\bm X}}{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$. Given the observation $\{(Y_i,{{\bm X}}_i,A_i);\;i=1,\ldots,n\}$, the Q-learning procedure estimates the parameters $({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})$ by minimizing the squared error loss $$L_{1n}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})=\frac{1}{n}\sum_{i=1}^{n}\left\{Y_i-h({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})-A_iC({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})\right\}^2.$$ Denote the optimized point as $(\hat{{{\mbox{\boldmath $\beta$}}}}^Q,\hat{{{\mbox{\boldmath $\gamma$}}}}^Q)$. The estimated optimal ITR based on Q-learning is then $\hat{g}^{Q}({{\bm x}})\triangleq{\mathrm{1}}\{C({{\bm x}};\hat{{{\mbox{\boldmath $\beta$}}}}^Q)>0\}$, which is a consistent estimator of $g^{\mathrm{opt}}_\mu({{\bm x}})$ if both $h({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ are correctly specified.

A-learning is a semiparametric improvement of Q-learning by modeling only the contrast function $C_0({{\bm X}})$ rather than the full Q-function. This is reasonable based on the observation that the optimal ITR $g^{\mathrm{opt}}_\mu$ only depends on $C_0({{\bm X}})$. By positing $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ for the contrast function, in A-learning, one can estimate coefficients ${{\mbox{\boldmath $\beta$}}}$ by solving the following estimating equation $$\sum_{i=1}^{n}\lambda({{\bm X}}_i)\left\{A_i-\pi({{\bm X}}_i)\right\}\left\{Y_i-A_iC({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})-h({{\bm X}}_i)\right\}=0,
\label{eq:eeofAlearn}$$ where $\lambda({{\bm X}}_i)$ and $h({{\bm X}}_i)$ are arbitrary functions, and $\lambda({{\bm X}}_i)$ has the same dimension as ${{\mbox{\boldmath $\beta$}}}$. Denote the solution to by $\hat{{{\mbox{\boldmath $\beta$}}}}^{A}$. If $\mathrm{var}(Y|X)$ is constant and $C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})$ is correctly specified, the optimal choices of $\lambda(\cdot)$ and $h(\cdot)$ are $\lambda({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})=\partial/\partial {{\mbox{\boldmath $\beta$}}}C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})$ and $h({{\bm X}}_i)=h_0({{\bm X}}_i)$ [@robins2004optimal]. In practice, one may pose models, say $\pi({{\bm X}}_i;{{\mbox{\boldmath $\phi$}}})$ and $h({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})$ for $\pi({{\bm X}}_i)$ and $h({{\bm X}}_i)$ respectively, and take $\lambda({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})=\partial/\partial {{\mbox{\boldmath $\beta$}}}C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})$. Under randomized designs, the propensity score $\pi({{\bm X}}_i)$ is known. Otherwise, a logistic model can be proposed. Under the assumption that $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ is correctly specified, the double robustness property of A-learning states that as long as one of $\pi({{\bm X}};{{\mbox{\boldmath $\phi$}}})$ and $h({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$ is correctly specified, $\hat{g}^{A}({{\bm x}})\triangleq{\mathrm{1}}\{C({{\bm x}};\hat{{{\mbox{\boldmath $\beta$}}}}^A)>0\}$ is consistent estimator of $g^{\mathrm{opt}}_\mu({{\bm x}})$.

Recently, [-@lu2011variable] propose a variant of A-learning by a loss-based learning framework. Rewrite $$\begin{aligned}
{\mbox{E}}(Y|{{\bm X}},A) =& h_0({{\bm X}})+AC_0({{\bm X}})\\
=& \varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}C_0({{\bm X}}),\end{aligned}$$ where $\varphi_0({{\bm X}})=h_0({{\bm X}})+\pi({{\bm X}})C_0({{\bm X}})$. Based on the expression above, [@lu2011variable] propose to estimate $({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})$ by minimizing the following loss function $$L_{2n}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})=\frac{1}{n}\sum_{i=1}^{n}\left[Y_i-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})-\{A_i-\pi({{\bm X}}_i)\}C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})\right]^2,
\label{eq:A-loss}$$ where $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$, $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ are proposed models for $\varphi_0({{\bm X}})$ and $C_0({{\bm X}})$ respectively. Denote the minimizer of as $(\hat{{{\mbox{\boldmath $\beta$}}}}^A_{LS},\hat{{{\mbox{\boldmath $\gamma$}}}}^{A}_{LS})$. [@lu2011variable] show that $\hat{g}^{A}_{LS}({{\bm x}})\triangleq{\mathrm{1}}\{C({{\bm x}};\hat{{{\mbox{\boldmath $\beta$}}}}^A_{LS})>0\}$ is a consistent estimator of $g^{\mathrm{opt}}_\mu({{\bm x}})$ when the propensity score $\pi({{\bm X}})$ is known or can be consistently estimated from the data, and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ is correctly specified. We refer to this method as least square A-learning (lsA-learning).

One main advantage of the lsA-learning, compared to the classical A-learning, is its square loss, making the procedure easy to be coupled with penalized regression to achieve variable selection in high dimensional data. Specifically, [@lu2011variable] propose to identify important nonzero coefficients in ${{\mbox{\boldmath $\beta$}}}$ by applying an adaptive LASSO penalty to . Under some regularity conditions, both the selection consistency and asymptotic normality of the estimator are established in [@lu2011variable]. The downside of lsA-learning is that one direction of the double robustness property of the classical A-learning is lost, i.e., when $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$ is correctly specified, ${{\mbox{\boldmath $\beta$}}}$ may still not be consistent if the propensity score $\pi({{\bm X}})$ is not consistently estimated. Finally, it can be shown that lsA-learning and Q-learning are equivalent when $\pi({{\bm X}})$ is constant and both $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ take the linear form (with the space of $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ included in the space of $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$). Similar properties hold for A-learning and Q-learning [@schulte2014q].

New Proposal: Robust Regression
-------------------------------

Skewed, heavy-tailed, heteroscedastic errors or outliers of the response $Y$ are frequently encountered in clinical trials. It is well known that ordinary least square estimation fails to produce a reliable estimator in such situations. The immediate consequence is the efficiency loss in the estimators produced by Q-, A-, and lsA-learning. This motivates us to adopt robust regression techniques in optimal treatment regime estimation.

We consider the following additive model, $$Y_i=\varphi_0({{\bm X}}_i)+\{A_i-\pi({{\bm X}}_i)\}C({{\bm X}}_i;\beta_0)+\epsilon_i,\; i=1,\ldots,n,
\label{eq:model_nointeraction}$$ where $\varphi_0({{\bm X}})$ is the baseline function, $C({{\bm X}};\beta_0)$ is the contrast function, $\pi({{\bm X}})$ is the propensity score, and $\epsilon$ is the error term which satisfies the conditional independence assumption $\epsilon\perp A|{{\bm X}}$. We point out that the error term defined in can be very general. For example, we could take $\epsilon=\sum_{j=1}^{K}\sigma_j({{\bm X}})e_j$ for any $K\geq 1$ that allows the error distribution to change with ${{\bm X}}$, used to model heterogeneous errors, where $\sigma_j({{\bm X}})$ are arbitrary positive functions and $e_{j}\perp (A,{{\bm X}})$ for all $j=1,\ldots,K$. Throughout the paper, we assume $\{(Y_i,{{\bm X}}_i,A_i,\epsilon_i),i=1,\ldots,n\}$ are i.i.d random samples of the population.

We propose to estimate $({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})$ by minimizing $$L_{3n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})=\frac{1}{n}\sum_{i=1}^{n}M\left[Y_i-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})-\{A_i-\pi({{\bm X}}_i)\}C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})\right],
\label{eq:A-general-loss}$$ where ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$, ${{\mbox{\boldmath $\beta$}}}\in\mathcal{B}$ and $M:{\mathrm{I \! R} \mathit{^{\rightarrow}}} [0,\infty)$ is a convex function with minimum achieved at 0. Denote the minimizer of as $(\hat{{{\mbox{\boldmath $\beta$}}}}^R_{M},\hat{{{\mbox{\boldmath $\gamma$}}}}^{R}_{M})$, and the estimated ITR is then $\hat{g}^{R}_M({{\bm x}})\triangleq{\mathrm{1}}\{C({{\bm x}};\hat{{{\mbox{\boldmath $\beta$}}}}^R_M)>0\}$. In the following, we refer the robust regression with loss function $M(x)$ as RR(M)-learning. In this article, we consider the following three types of loss functions, i.e., the pinball loss $$M(x)=\rho_\tau(x)\triangleq
\begin{cases}
(\tau-1)x, &\text{if } x<0\\
\tau x, &\text{if } x\geq0
\end{cases}$$ where $0<\tau<1$, the Huber loss $$M(x)=H_\alpha(x)\triangleq
\begin{cases}
0.5x^2, &\text{if } |x|<\alpha\\
\alpha|x|-0.5\alpha^2, &\text{if } |x|\geq\alpha
\end{cases}$$ for some $\alpha>0$, and the $\epsilon$-insensitive loss $$M(x)=J_\epsilon(x)\triangleq\max(0, |x|-\epsilon)$$ for some $\epsilon>0$. The pinball loss are frequently applied for quantile regression [@koenker2005quantile], and the Huber losses and the $\epsilon$-insensitive are robust against heavy tailed errors or outliers. A dramatic difference of pinball loss, Huber loss and $\epsilon$-insensitive loss, compared with the square loss, is that they penalize large deviances linearly instead of quadratically. This property makes them more robust when dealing with responses with non-normal type of errors.

Asymptotic Properties {#section:asymptotic}
=====================

Consistency of Robust Regression: Pinball Loss
----------------------------------------------

Under the conditional independence assumption $\epsilon\perp A|{{\bm X}}$, we have $$\begin{aligned}
Q({{\bm X}},A)=& \varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}C({{\bm X}};\beta_0)+\mu_\epsilon({{\bm X}});\\
Q_\tau({{\bm X}},A)=& \varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}C({{\bm X}};\beta_0)+F^{-1}_{\epsilon}({{\bm X}};\tau).\end{aligned}$$ where $\mu_{\epsilon}({{\bm X}})$ and $F^{-1}_{\epsilon}({{\bm X}};\tau)$ denote the mean and the $\tau$-th quantile of $\epsilon$ conditional on ${{\bm X}}$ respectively. Therefore, in this situation, we have $g^{\mathrm{opt}}_\mu=g^{\mathrm{opt}}_{\tau}={\mathrm{1}}\{C({{\bm X}};{{\mbox{\boldmath $\beta$}}}_0)>0\}$. In other words, the underlying ITR which maximize the population mean and $\tau$-th quantile are equivalent. For a good ITR $\hat{g}={\mathrm{1}}\{C({{\bm X}};\hat{{{\mbox{\boldmath $\beta$}}}})>0\}$, it is reasonable to require $\hat{{{\mbox{\boldmath $\beta$}}}}$ to be a consistent estimator of ${{\mbox{\boldmath $\beta$}}}_0$. This consistency result is first shown for the robust regression with pinball loss, which is given in Theorem 1. We allocate all the proofs into the Appendix A.

Under regularity conditions (A1)-(A8) in the Appendix A, if the contrast function in is correctly specified and $\pi({{\bm x}})$ is known, then $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)}\inprob{{\mbox{\boldmath $\beta$}}}_0$ for all $\tau\in(0,1)$, where $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)}$ is the solution of when $M(x)=\rho_{\tau}(x)$.

**Remarks:**

1.  Theorem 1 doesn’t assume the finiteness of $E(Y)$. Therefore it can be applied to the cases when $\epsilon_i$ follows a Cauchy distribution.

2.  After fitting the model, the Assumption (A2), $\epsilon\perp A|{{\bm X}}$, can be verified by applying conditional independence test with $\hat{r}(\hat{{{\mbox{\boldmath $\beta$}}}},\hat{{{\mbox{\boldmath $\gamma$}}}})$ and $A$ given ${{\bm X}}$, where $\hat{r}(\hat{{{\mbox{\boldmath $\beta$}}}},\hat{{{\mbox{\boldmath $\gamma$}}}})$ is the estimated residual and $\hat{r}(\hat{{{\mbox{\boldmath $\beta$}}}},\hat{{{\mbox{\boldmath $\gamma$}}}})=Y-\varphi({{\bm X}};\hat{{{\mbox{\boldmath $\gamma$}}}})-\{A-\pi({{\bm X}})\}C({{\bm X}};\hat{{{\mbox{\boldmath $\beta$}}}})$. See [@lawrance1976conditional; @su2007consistent; @song2007testing; @huang2010testing; @zhang2012kernel] for more discussion of conditional independence hypothesis tests. In particular, we demonstrate the usefulness of the test by applying the Kernel-based conditional independence test (KCI-test, [@zhang2012kernel]) in Section 5. KCI-test doesn’t assume functional forms among variables and thus suits our need.

When the conditional independence assumption ($\epsilon\perp A|{{\bm X}}$) does not hold, $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)}$ may no longer be a consistent estimator of ${{\mbox{\boldmath $\beta$}}}_0$. This is intuitively reasonable as $\epsilon$ contains extra information with respect to $A$. In fact, a general result which can be derived in this case is that, $(\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)},\hat{{{\mbox{\boldmath $\gamma$}}}}^{R}_{\rho(\tau)})$ minimizes a weighed mean-square error loss function with specification error [@angrist2006quantile; @lee2013interpretation].

Instead of assuming response $Y$ takes an additive error term $\epsilon$ as in , we assume the conditional quantile function $Q_\tau({{\bm X}},A)=\varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}C({{\bm X}};{{\mbox{\boldmath $\beta$}}}_0(\tau))$, where we redundantly represent the baseline function and contrast function as $\varphi_0(\cdot)$ and $C(\cdot)$ respectively. Notice that we use ${{\mbox{\boldmath $\beta$}}}_0(\tau)$ instead of ${{\mbox{\boldmath $\beta$}}}_0$ to emphasize that the true ${{\mbox{\boldmath $\beta$}}}$ may vary with respect to $\tau$. The proposed model is $\hat{Q}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})=\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\{A-\pi({{\bm X}})\}C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ with $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ correctly specified. Define $$\label{eq:population_beta_tau}
\left({{\mbox{\boldmath $\beta$}}}(\tau),{{\mbox{\boldmath $\gamma$}}}(\tau)\right)=\underset{{{\mbox{\boldmath $\beta$}}}\in\mathcal{B},{{\mbox{\boldmath $\gamma$}}}\in\Gamma}{\mathrm{argmin}}{\mbox{E}}\left[\rho_{\tau}\{Y-\hat{Q}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\}
-\rho_{\tau}\{Y-\hat{Q}({{\mbox{\boldmath $\beta$}}}',{{\mbox{\boldmath $\gamma$}}}')\}\right]$$ where $({{\mbox{\boldmath $\beta$}}}',{{\mbox{\boldmath $\gamma$}}}')$ is any fixed point in $\mathcal{B}\times\Gamma$. Define the QR specification error as $\Delta_{\tau}({{\bm X}},A;{{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\triangleq\hat{Q}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})-Q_\tau({{\bm X}},A)$. Define the quantile-specific residual as $\epsilon_{\tau}\triangleq Y-Q_\tau({{\bm X}},A)$ with conditional density function $f_{\epsilon_{\tau}}(\cdot|{{\bm X}},A)$. Then we have the following approximation theorem. The proof of the theorem follows Theorem 1 of [@angrist2006quantile], and is omitted for brevity.

Suppose that (i) the conditional density $f_{Y}(y|{{\bm X}},A)$ exists a.s.; (ii)${\mbox{E}}[Q_{\tau}({{\bm X}},A)]$ and ${\mbox{E}}[\Delta^2_{\tau}({{\bm X}},A;{{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})]$ are finite; (iii) $\left({{\mbox{\boldmath $\beta$}}}(\tau),{{\mbox{\boldmath $\gamma$}}}(\tau)\right)$ uniquely solves . Then $$\left({{\mbox{\boldmath $\beta$}}}(\tau),{{\mbox{\boldmath $\gamma$}}}(\tau)\right)=\underset{{{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}}}{\mathrm{argmin}}{\mbox{E}}[w_{\tau}({{\bm X}},A;{{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\Delta^2_{\tau}({{\bm X}},A;{{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})]$$ where $$w_{\tau}({{\bm X}},A;{{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})=\int_{0}^{1}(1-u)f_{\epsilon_{\tau}}(u\Delta_{\tau}({{\bm X}},A;{{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})|{{\bm X}},A)du.$$ \[thm:approximation\_quantile\]

**Remarks:**

1.  Theorem 2 shows that $\hat{Q}\left({{\mbox{\boldmath $\beta$}}}(\tau),{{\mbox{\boldmath $\gamma$}}}(\tau)\right)$ is a weighted least square approximation to $Q_{\tau}({{\bm X}},A)$. In other word, $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}}(\tau))+\{A-\pi({{\bm X}})\}C({{\bm X}};{{\mbox{\boldmath $\beta$}}}(\tau))$ is close to $\varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}C({{\bm X}};{{\mbox{\boldmath $\beta$}}}_0(\tau))$. So even though it is not true that ${{\mbox{\boldmath $\beta$}}}(\tau)={{\mbox{\boldmath $\beta$}}}_0(\tau)$ holds exactly, the difference between them is small in general . This coupled with the fact that $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)}\inprob{{\mbox{\boldmath $\beta$}}}(\tau)$ (proved in Theorem \[thm:asymptotic\_normality\_pinball\]), leads to the conclusion that approximately ITR $\hat{g}^{R}_{\rho(\tau)}({{\bm x}})$ $(\triangleq{\mathrm{1}}\{C({{\bm x}};\hat{{{\mbox{\boldmath $\beta$}}}}^R_{\rho(\tau)})>0\})$ maximizes the $\tau$-th conditional quantile. This observation is justified numerically in Section 4.2.

2.  When there exists ${{\mbox{\boldmath $\gamma$}}}_0\in\Gamma$ such that $\varphi_0({{\bm X}})\equiv\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}}_0)$, then we have ${{\mbox{\boldmath $\beta$}}}(\tau)={{\mbox{\boldmath $\beta$}}}_0(\tau)$.

Consistency of Robust Regression: Other Losses
----------------------------------------------

Under model and the assumption $\epsilon\perp A|{{\bm X}}$, similar consistency results can be established for Huber loss and the $\epsilon$-insensitive loss, as stated in Theorem \[thm:consistency\_huber\].

\[thm:consistency\_huber\] Under regularity conditions (A1)-(A8), if the contrast function in is correctly specified and $\pi({{\bm x}})$ is known, then we have

(a) $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{H(\alpha)}\inprob{{\mbox{\boldmath $\beta$}}}_0$ for all $\alpha>0$, where $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{H(\alpha)}$ is the solution of when $M(x)=H_{\alpha}(x)$;

(b) $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{J(\epsilon)}\inprob{{\mbox{\boldmath $\beta$}}}_0$ for all $\epsilon>0$, where $\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{J(\epsilon)}$ is the solution of when $M(x)=J_{\epsilon}(x)$.

Asymptotic Normality: Pinball Loss
----------------------------------

Without loss of generality, in this section we assume both the $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$ take the linear form: $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})=\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})=\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\beta$}}}$, where $\tilde{{{\bm X}}}=(1,{{\bm X}}{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$. Denote $\hat{{{\mbox{\boldmath $\beta$}}}}(\tau)=\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)}$ and $\hat{{{\mbox{\boldmath $\gamma$}}}}(\tau)=\hat{{{\mbox{\boldmath $\gamma$}}}}^{R}_{\rho(\tau)}$. Denote ${{\bm W}}=(\{A-\pi({{\bm X}})\}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}},\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$, ${{\mbox{\boldmath $\theta$}}}(\tau)=({{\mbox{\boldmath $\beta$}}}(\tau){^{\mbox{\tiny {\sf T}}}},{{\mbox{\boldmath $\gamma$}}}(\tau){^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$, $\hat{{{\mbox{\boldmath $\theta$}}}}(\tau)=(\hat{{{\mbox{\boldmath $\beta$}}}}(\tau){^{\mbox{\tiny {\sf T}}}},\hat{{{\mbox{\boldmath $\gamma$}}}}(\tau){^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$ and $J(\tau)\triangleq{\mbox{E}}\left[f_Y({{\bm W}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\theta$}}}(\tau)|{{\bm X}},A){{\bm W}}{{\bm W}}{^{\mbox{\tiny {\sf T}}}}\right]$. Under the following regularity conditions, which is the same as the assumptions assumed in [@angrist2006quantile] and [@lee2013interpretation], we have the asymptotic normality of $\hat{{{\mbox{\boldmath $\theta$}}}}(\tau)$, which is given in Theorem \[thm:asymptotic\_normality\_pinball\].

-   $\{(Y_i,{{\bm X}}_i,A_i,\epsilon_i),i=1,\ldots,n\}$ are i.i.d random variables;

-   the conditional density $f_Y(y|{{\bm X}}={{\bm x}},A=a))$ exists, and is bounded and uniformly continuous in y, uniformly in ${{\bm x}}$ over the support of ${{\bm X}}$;

-   $J(\tau)$ is positive definite for all $\tau\in(0,1)$, where ${{\mbox{\boldmath $\theta$}}}(\tau)$ is uniquely defined in ;

-   ${\mbox{E}}\|{{\bm X}}\|^{2+\epsilon}$ for some $\epsilon>0$.

\[thm:asymptotic\_normality\_pinball\] If regularity conditions (B1)-(B4) are hold, we have

1.  (**Uniform Consistency**) $\sup_{\tau}\|\hat{{{\mbox{\boldmath $\theta$}}}}(\tau)-{{\mbox{\boldmath $\theta$}}}(\tau)\|=o_p(1)$;

2.  (**Asymptotic Normality**) $J(\cdot)\sqrt{n}(\hat{{{\mbox{\boldmath $\theta$}}}}(\cdot)-{{\mbox{\boldmath $\theta$}}}(\cdot))$ converge in distribution to a zero mean Gaussian process with covariance function $\Sigma(\tau,\tau')$ defined as $$\Sigma(\tau,\tau')={\mbox{E}}\left[\left(\tau-{\mathrm{1}}\{Y<{{\bm W}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\theta$}}}(\tau)\}\right)
            \left(\tau'-{\mathrm{1}}\{Y<{{\bm W}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\theta$}}}(\tau)\}\right){{\bm W}}{{\bm W}}{^{\mbox{\tiny {\sf T}}}}\right].$$

The proof is given in [@angrist2006quantile], and the asymptotic covariance matrix of $\hat{{{\mbox{\boldmath $\theta$}}}}(\tau)$ can be estimated by either a bootstrap procedure [@hahn1997bayesian] or a nonparametric kernel method [@angrist2006quantile]. We adopt the parametric bootstrap approach to estimate the asymptotic covariance matrix in Section 5. Under model the result of Theorem \[thm:asymptotic\_normality\_pinball\] can be further simplified, which is given in Theorem 5.

Under the condition of Theorem 4, if further we assume $Y=\varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\beta$}}}_0+\epsilon$, and $\epsilon\perp A|{{\bm X}}$, then

1.  $\sup_{\tau}\|\hat{{{\mbox{\boldmath $\beta$}}}}(\tau)-{{\mbox{\boldmath $\beta$}}}_0\|=o_p(1)$;

2.  $\sqrt{n}(\hat{{{\mbox{\boldmath $\beta$}}}}(\tau)-{{\mbox{\boldmath $\beta$}}}_0)\indist N({{\mathbf 0}}, J_{11}^{-1}(\tau)\Sigma_{11}(\tau,\tau)J_{11}^{-1}(\tau))$, where $$\begin{aligned}
            J_{11}(\tau)=& {\mbox{E}}\left[f_{\epsilon}\left(\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})|{{\bm X}}\right)
            \pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\right],\\
            \Sigma_{11}(\tau,\tau)=& {\mbox{E}}\left\{\left[\tau-{\mathrm{1}}\{\epsilon<\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})\}\right]^2
            \pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\right\}.
            \end{aligned}$$ Furthermore, we have $\Sigma_{11}(\tau,\tau)\leq\left(\tau^2+|1-2\tau|\right){\mbox{E}}\left[\pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\right]$.

Comparing the asymptotic normality of $\hat{{{\mbox{\boldmath $\beta$}}}}(\tau)$ with $\hat{{{\mbox{\boldmath $\beta$}}}}^A_{LS}$ yields interesting insights. Assuming that ${\mbox{E}}(Y|{{\bm X}},A)=\varphi_0({{\bm X}})+\{A-\pi({{\bm X}})\}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\beta$}}}_0$ and $({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}}^*)={\mathrm{argmin}}_{({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})}{\mbox{E}}[Y-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})-\{A-\pi({{\bm X}})\}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\beta$}}}]^2$, the asymptotic normality property of $\hat{{{\mbox{\boldmath $\beta$}}}}^{A}_{LS}$ can then be established, which is summarized in Theorem \[thm:lsA\]. Its proof has been omitted, and readers are referred to [@lu2011variable].

\[thm:lsA\] Under the regularity condition of A1-A4 of [@lu2011variable], $$\sqrt{n}(\hat{{{\mbox{\boldmath $\beta$}}}}^{A}_{LS}-{{\mbox{\boldmath $\beta$}}}_0)\indist N(0,U_{11}^{-1}\Omega_{11}U_{11}^{-1}),$$ where $U_{11}={\mbox{E}}\left[\pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\right]$ and $$\Omega_{11}={\mbox{E}}\left[\left\{\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}}^*)+\epsilon\right\}^2
    \pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\right]$$

**Remarks:**

1.  When the family of functions $\{\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}}),{{\mbox{\boldmath $\gamma$}}}\in\Gamma\}$ cannot well approximate the unknown baseline function $\varphi_0({{\bm X}})$, the $\Omega_{11}$ term in the asymptotic variance of $\hat{{{\mbox{\boldmath $\beta$}}}}^{A}_{LS}$ may explode, which makes $\hat{{{\mbox{\boldmath $\beta$}}}}^{A}_{LS}$ less efficient than $\hat{{{\mbox{\boldmath $\beta$}}}}(\tau)$.

2.  When $Y=\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\gamma_0+\{A-\pi({{\bm X}})\}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\beta$}}}_0+\epsilon$, $\epsilon\perp(A,{{\bm X}})$, $\pi({{\bm X}})\equiv0.5$ and $\epsilon\sim N(0,\sigma^2)$, the asymptotic variance of $\hat{{{\mbox{\boldmath $\beta$}}}}(\tau=0.5)$ is $2\pi\sigma^2{\mbox{E}}(\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}})^{-1}$, which is strictly larger than $4\sigma^2{\mbox{E}}(\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}})^{-1}$ (the asymptotic variance of $\hat{{{\mbox{\boldmath $\beta$}}}}^{A}_{LS}$).

Numerical Results: Simulation Studies {#section:simulation}
=====================================

To demonstrate finite sample performance of the proposed robust regression methods for optimal treatment rule estimation, we conduct two simulation studies: the errors independent with treatments, and the errors interactive with treatments, respectively.

Simulation Study I: error terms independent with treatment
----------------------------------------------------------

We consider the following two models with p=3,

-   Model I: $$Y_i=1+(X_{i1}-X_{i2})(X_{i1}+X_{i3})+\{A_i-\pi({{\bm X}}_{i})\}{{\mbox{\boldmath $\beta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\sigma({{\bm X}}_{i})\epsilon_i,$$ where ${{\bm X}}_{i}=(X_{i1},X_{i2},X_{i3}){^{\mbox{\tiny {\sf T}}}}$ are multivariate normal with mean 0, variance 1, and $\mathrm{Corr}(X_{ij},X_{ik})=0.5^{|j-k|}$, $\tilde{{{\bm X}}}_i=(1,{{\bm X}}_i{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$ and ${{\mbox{\boldmath $\beta$}}}_0=(0,1,-1,1){^{\mbox{\tiny {\sf T}}}}$.

-   Model II: $$Y_i={{\mbox{\boldmath $\gamma$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\{A_i-\pi({{\bm X}}_{i})\}{{\mbox{\boldmath $\beta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\sigma({{\bm X}}_{i})\epsilon_i,$$ where ${{\mbox{\boldmath $\gamma$}}}_0{^{\mbox{\tiny {\sf T}}}}=(0.5,4,1,-3)$, and ${{\bm X}}_{i}$, $\tilde{{{\bm X}}}_{i}$ and ${{\mbox{\boldmath $\beta$}}}_0$ are the same as Model I.

We take linear forms for both the baseline and the contrast functions, where $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})={{\mbox{\boldmath $\gamma$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})={{\mbox{\boldmath $\beta$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$. We assume the propensity scores $\pi(\cdot)$ are known, and we study both the constant case $(\pi({{\bm X}}_i)=0.5)$ and the non-constant case $(\pi({{\bm X}}_i)=\mathrm{logit}({{\bm X}}_{i1}-{{\bm X}}_{i2}))$. In addition, We consider two different $\sigma({{\bm X}}_{i})$ functions, i.e., the homogeneous case with $\sigma({{\bm X}}_{i})=1$, and the heterogenous case with $\sigma({{\bm X}}_{i})=0.5+(X_{i1}-X_{i2})^2$. The simulation results under constant and non-constant propensity scores are similar. Thus, for brevity, we only report the constant case and allocate the result of non-constant case to the Appendix B. The results of Model I and II with constant propensity score are given in Table \[table:modelI\_constant\] and \[table:modelII\_constant\] respectively.

[@ll ccc ccc ccc]{}\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 1.32 (0.040) & 80.7 & 1.06 & 2.36 (0.081) & 75.7 & 1.57 & & 58.4 & 3.75\
& P(0.5) & 1.44 (0.042) & 80.1 & 1.13 & 1.73 (0.051) & 78.0 & 1.31 & 2.69 (0.077) & 75.2 & 1.63\
& P(0.25) & 1.90 (0.057) & 78.3 & 1.34 & 1.63 (0.051) & 79.0 & 1.29 & 5.29 (0.168) & 70.4 & 2.25\
& Huber & 1.15 (0.034) & 81.9 & 0.93 & 1.45 (0.044) & 79.9 & 1.13 & 2.61 (0.072) & 74.9 & 1.66\
200 & LS & 0.68 (0.021) & 85.6 & 0.59 & 1.10 (0.033) & 82.0 & 0.91 & & 58.7 & 3.70\
& P(0.5) & 0.73 (0.021) & 85.3 & 0.62 & 0.78 (0.021) & 84.1 & 0.70 & 1.23 (0.037) & 81.3 & 0.99\
& P(0.25) & 0.92 (0.028) & 84.0 & 0.75 & 0.70 (0.023) & 86.0 & 0.59 & 2.48 (0.079) & 75.7 & 1.64\
& Huber & 0.58 (0.017) & 86.8 & 0.50 & 0.66 (0.018) & 85.5 & 0.58 & 1.24 (0.035) & 80.8 & 1.03\
400 & LS & 0.33 (0.009) & 90.3 & 0.26 & 0.56 (0.016) & 87.1 & 0.46 & & 59.2 & 3.61\
& P(0.5) & 0.35 (0.010) & 90.0 & 0.29 & 0.37 (0.010) & 89.0 & 0.34 & 0.56 (0.016) & 87.1 & 0.48\
& P(0.25) & 0.43 (0.013) & 89.1 & 0.34 & 0.33 (0.010) & 90.7 & 0.25 & 1.16 (0.037) & 82.9 & 0.86\
& Huber & 0.28 (0.008) & 91.1 & 0.22 & 0.31 (0.009) & 90.2 & 0.27 & 0.58 (0.017) & 86.7 & 0.49\
800 & LS & 0.17 (0.005) & 93.2 & 0.13 & 0.26 (0.008) & 90.9 & 0.23 & & 59.4 & 3.59\
& P(0.5) & 0.17 (0.005) & 93.1 & 0.13 & 0.19 (0.005) & 92.1 & 0.17 & 0.29 (0.009) & 90.7 & 0.24\
& P(0.25) & 0.22 (0.007) & 92.4 & 0.16 & 0.18 (0.006) & 93.6 & 0.12 & 0.59 (0.019) & 87.3 & 0.48\
& Huber & 0.14 (0.004) & 93.8 & 0.11 & 0.16 (0.005) & 93.1 & 0.14 & 0.29 (0.008) & 90.5 & 0.25\
\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 3.24 (0.110) & 74.7 & 1.70 & 8.98 (0.561) & 68.6 & 2.44 & & 56.2 & 4.05\
& P(0.5) & 1.70 (0.060) & 80.5 & 1.08 & 1.80 (0.064) & 80.1 & 1.08 & 3.45 (0.124) & 75.1 & 1.69\
& P(0.25) & 2.50 (0.085) & 77.4 & 1.42 & 2.51 (0.079) & 76.8 & 1.46 & 9.13 (0.341) & 67.2 & 2.66\
& Huber & 1.70 (0.057) & 80.4 & 1.10 & 1.87 (0.063) & 79.2 & 1.16 & 4.27 (0.155) & 72.8 & 1.93\
200 & LS & 1.54 (0.050) & 80.6 & 1.06 & 4.71 (0.244) & 73.4 & 1.85 & & 55.2 & 4.17\
& P(0.5) & 0.78 (0.028) & 86.7 & 0.53 & 0.90 (0.032) & 85.3 & 0.63 & 1.49 (0.052) & 81.9 & 0.95\
& P(0.25) & 1.16 (0.039) & 83.5 & 0.81 & 1.23 (0.039) & 82.0 & 0.91 & 3.95 (0.150) & 73.2 & 1.90\
& Huber & 0.77 (0.025) & 86.4 & 0.55 & 0.94 (0.032) & 84.5 & 0.69 & 1.94 (0.071) & 79.3 & 1.19\
400 & LS & 0.80 (0.026) & 86.0 & 0.58 & 2.69 (0.136) & 77.8 & 1.34 & & 54.7 & 4.26\
& P(0.5) & 0.39 (0.013) & 90.5 & 0.27 & 0.44 (0.017) & 89.6 & 0.32 & 0.71 (0.024) & 86.9 & 0.50\
& P(0.25) & 0.56 (0.019) & 88.8 & 0.37 & 0.66 (0.020) & 86.9 & 0.50 & 1.70 (0.055) & 79.6 & 1.17\
& Huber & 0.38 (0.012) & 90.4 & 0.27 & 0.48 (0.017) & 88.8 & 0.36 & 0.91 (0.029) & 84.9 & 0.65\
800 & LS & 0.41 (0.013) & 89.9 & 0.29 & 1.35 (0.150) & 83.1 & 0.82 & & 56.5 & 4.00\
& P(0.5) & 0.18 (0.006) & 93.6 & 0.12 & 0.20 (0.007) & 92.6 & 0.16 & 0.36 (0.013) & 91.0 & 0.25\
& P(0.25) & 0.28 (0.009) & 92.2 & 0.18 & 0.31 (0.010) & 90.8 & 0.24 & 0.89 (0.031) & 85.8 & 0.60\
& Huber & 0.19 (0.006) & 93.3 & 0.13 & 0.22 (0.007) & 92.1 & 0.18 & 0.47 (0.017) & 89.2 & 0.34\

[@ll ccc ccc ccc]{}\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 0.24 (0.006) & 91.1 & 0.21 & 1.23 (0.061) & 82.4 & 0.87 & & 58.6 & 3.73\
& P(0.5) & 0.36 (0.010) & 89.0 & 0.32 & 0.39 (0.012) & 88.8 & 0.34 & 0.80 (0.024) & 84.2 & 0.69\
& P(0.25) & 0.45 (0.012) & 87.8 & 0.40 & 0.13 (0.004) & 93.4 & 0.12 & 2.37 (0.083) & 76.0 & 1.49\
& Huber & 0.25 (0.007) & 90.8 & 0.22 & 0.31 (0.010) & 90.3 & 0.26 & 0.99 (0.029) & 82.4 & 0.84\
200 & LS & 0.11 (0.003) & 93.7 & 0.10 & 0.52 (0.018) & 87.3 & 0.45 & & 58.7 & 3.69\
& P(0.5) & 0.17 (0.005) & 92.4 & 0.16 & 0.17 (0.005) & 92.4 & 0.15 & 0.32 (0.009) & 89.5 & 0.30\
& P(0.25) & 0.20 (0.005) & 91.8 & 0.18 & 0.06 (0.002) & 95.6 & 0.05 & 1.03 (0.033) & 82.1 & 0.88\
& Huber & 0.12 (0.003) & 93.6 & 0.11 & 0.13 (0.003) & 93.5 & 0.12 & 0.43 (0.013) & 87.9 & 0.40\
400 & LS & 0.05 (0.001) & 95.7 & 0.05 & 0.26 (0.008) & 90.7 & 0.23 & & 59.4 & 3.60\
& P(0.5) & 0.09 (0.002) & 94.5 & 0.08 & 0.09 (0.002) & 94.5 & 0.08 & 0.15 (0.004) & 92.8 & 0.14\
& P(0.25) & 0.10 (0.002) & 94.2 & 0.09 & 0.03 (0.001) & 96.9 & 0.02 & 0.44 (0.012) & 87.9 & 0.39\
& Huber & 0.06 (0.001) & 95.5 & 0.05 & 0.06 (0.002) & 95.4 & 0.06 & 0.21 (0.006) & 91.6 & 0.19\
800 & LS & 0.03 (0.001) & 96.9 & 0.03 & 0.13 (0.004) & 93.5 & 0.11 & & 59.4 & 3.58\
& P(0.5) & 0.04 (0.001) & 96.1 & 0.04 & 0.04 (0.001) & 96.2 & 0.04 & 0.07 (0.002) & 95.1 & 0.06\
& P(0.25) & 0.05 (0.001) & 95.8 & 0.05 & 0.01 (0.000) & 97.9 & 0.01 & 0.20 (0.005) & 91.5 & 0.19\
& Huber & 0.03 (0.001) & 96.8 & 0.03 & 0.03 (0.001) & 96.8 & 0.03 & 0.10 (0.002) & 94.2 & 0.09\
\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 1.97 (0.072) & 79.8 & 1.13 & 7.75 (0.514) & 70.4 & 2.22 & & 56.4 & 4.02\
& P(0.5) & 0.84 (0.029) & 86.1 & 0.55 & 1.21 (0.045) & 84.3 & 0.74 & 1.82 (0.071) & 80.5 & 1.07\
& P(0.25) & 1.37 (0.049) & 82.1 & 0.90 & 1.56 (0.051) & 80.5 & 1.04 & 6.20 (0.261) & 69.8 & 2.25\
& Huber & 0.84 (0.031) & 85.9 & 0.57 & 1.33 (0.046) & 82.8 & 0.85 & 2.69 (0.106) & 77.0 & 1.42\
200 & LS & 0.99 (0.035) & 84.7 & 0.66 & 4.16 (0.237) & 75.2 & 1.62 & & 55.1 & 4.19\
& P(0.5) & 0.41 (0.014) & 90.2 & 0.28 & 0.58 (0.024) & 89.4 & 0.37 & 0.79 (0.030) & 86.7 & 0.52\
& P(0.25) & 0.64 (0.021) & 87.4 & 0.45 & 0.74 (0.024) & 86.1 & 0.54 & 2.48 (0.096) & 76.9 & 1.40\
& Huber & 0.39 (0.013) & 90.3 & 0.27 & 0.69 (0.027) & 87.7 & 0.45 & 1.17 (0.044) & 83.4 & 0.78\
400 & LS & 0.51 (0.018) & 89.0 & 0.35 & 2.48 (0.133) & 79.3 & 1.20 & & 54.7 & 4.25\
& P(0.5) & 0.20 (0.007) & 93.2 & 0.14 & 0.29 (0.011) & 92.6 & 0.17 & 0.32 (0.011) & 91.2 & 0.22\
& P(0.25) & 0.30 (0.009) & 91.3 & 0.22 & 0.39 (0.012) & 89.9 & 0.28 & 0.99 (0.030) & 83.0 & 0.78\
& Huber & 0.20 (0.007) & 93.2 & 0.14 & 0.34 (0.012) & 91.4 & 0.22 & 0.53 (0.016) & 88.4 & 0.37\
800 & LS & 0.25 (0.008) & 92.2 & 0.17 & 1.25 (0.159) & 84.2 & 0.73 & & 56.4 & 4.00\
& P(0.5) & 0.10 (0.004) & 95.3 & 0.07 & 0.14 (0.006) & 94.7 & 0.09 & 0.16 (0.006) & 93.9 & 0.11\
& P(0.25) & 0.14 (0.005) & 94.0 & 0.10 & 0.18 (0.006) & 92.9 & 0.14 & 0.49 (0.015) & 88.0 & 0.39\
& Huber & 0.09 (0.004) & 95.3 & 0.06 & 0.17 (0.006) & 93.9 & 0.11 & 0.26 (0.009) & 91.8 & 0.19\

Comparison is made among four methods. They are: lsA-learning, robust regression with $\rho_{0.5}$ (RR($\rho_{0.5}$)), robust regression with $\rho_{0.25}$ (RR($\rho_{0.25}$)), and robust regression with Huber loss (RR(H)). The error terms $\epsilon_i$ are taken as standard i.i.d. normal, log-normal or Cauchy distribution, and independent with both $A$ and ${{\bm X}}$. It is easy to check that the conditional independence assumption $\epsilon\perp A|{{\bm X}}$ is satisfied, and $g^{\mathrm{opt}}_{\mu}=g^{\mathrm{opt}}_{\tau}={\mathrm{1}}\{{{\mbox{\boldmath $\beta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i>0\}$. We consider four different sample sizes 100, 200, 400 and 800. To evaluate the performance of each method, we compare three groups of criteria: (1) the mean squared error $\|\hat{{{\mbox{\boldmath $\beta$}}}}-{{\mbox{\boldmath $\beta$}}}_{0}\|^2_2$ (mse), which measures the distance between estimated parameters and the true parameter ${{\mbox{\boldmath $\beta$}}}_0$; (2) the percentage of making correct decisions (PCD), which are calculated based on a validation set with 10000 observations. Specifically, we take the formula $100*\left(1-\sum_{i=1}^{N_T}|{\mathrm{1}}\{\hat{{{\mbox{\boldmath $\beta$}}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i>0\}-{\mathrm{1}}\{{{\mbox{\boldmath $\beta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i>0\}|/N_T\right)$ with $N_T=10000$; (3) the differences of $V_\mu(g)$ and $V_{0.5-q}(g)$ between the optimal ITR and the estimated ITR, where $\delta_{\mu}=V_{\mu}(g_{\mu}^{\mathrm{opt}})-V_{\mu}(\hat{g})$ and $\delta_{\tau}=V_{\tau-q}(g_{\mu}^{\mathrm{opt}})-V_{\tau-q}(\hat{g})$, $\forall\tau\in(0,1)$. $V_{\mu}(g)$ and $V_{\tau-q}(g)$ (defined in Section 2.1) are estimated from the validation set as well, and they evaluate the overall performance of an ITR $g$, where the former one focuses on the response’s mean and the latter one focuses on the response’s conditional $\tau$-th quantile. Under our setting, $\delta_{\mu}=\delta_{0.5}$ when they both exists. Thus, only $\delta_{0.5}$ is reported. For each scenario, we take 1000 replications. All numbers in the tables are based on the sample average of all replications. We further report the standard errors of mse to evaluate the variability of the corresponding statistics.

When the propensity score is constant, lsA-learning is equivalent to both Q- and A-learning under our setting. If we compare the performance of the methods under homogeneous and heterogeneous errors, the first thing we find is that lsA-learning works much worse under the heterogeneous errors, while all other methods are generally less affected by the heterogeneity of the errors. When the baseline function is misspecified as in Model I, under the homogeneous normal errors, RR(H) works slightly better than lsA-learning, while $\mathrm{RR}(\rho_{0.25})$ works the worst. However, the difference in general is small. For the homogeneous log-normal errors, again RR(H) works the best, while $\mathrm{RR}(\rho_{0.5})$ and $\mathrm{RR}(\rho_{0.25})$ have similar performance, and lsA-learning works the worst. Under the homogeneous Cauchy errors, $\mathrm{RR}(\rho_{0.5})$ works the best and RR(H) has a close performance. The lsA-learning is no longer consistent, and its mse explodes. The actual numbers are too large and thus leave as blank in Table \[table:modelI\_constant\] and \[table:modelII\_constant\]. Furthermore, with the Cauchy errors, the PCD of lsA-learning are less than 60% under all scenarios, while other methods’ PCD can be as high as 90%. When baseline function is correctly specified as in Model II, under homogeneous normal errors, lsA-learning performs the best. However, in this case RR(H) also has a very close performance, and thus makes no difference from a practical point of view to choose between these two methods. The results of Model II under other cases draw similar conclusion as Model I. To sum up, the overall conclusion is that, under the conditional independence assumption, the proposed robust regression method RR(M) is more efficient than Q-, A- and lsA-learning in the circumstances when observations have skewed, heterogeneous or heavy-tailed errors. On the other hand, when the error terms indeed follows i.i.d. normal distribution, the loss of efficiency of RR(M) is not significant. This is especially true when Huber loss is applied.

Simulation Study II: error terms interactive with treatment
-----------------------------------------------------------

We consider the following model with p=2, $$Y_i=1 + 0.5\sin[\pi(X_{i1}-X_{i2})]+
0.25(1+X_{i1}+2X_{i2})^2+(A_i-\pi({{\bm X}}_{i})){{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\sigma({{\bm X}}_{i},A_i)\epsilon_i,$$ where ${{\bm X}}_{i}=(X_{i1},X_{i2}){^{\mbox{\tiny {\sf T}}}}$, $\tilde{{{\bm X}}}_i=(1,{{\bm X}}_i{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$, $\sigma({{\bm X}}_{i},A_i)=1+A_i d_0 X_{i1}^2$, ${{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}}=(0.5,2,-1)$ and $X_{ik}$ are i.i.d. Uniform\[-1,1\].

Similar as Section 4.1, we take linear forms for both the baseline and the contrast functions, where $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})={{\mbox{\boldmath $\gamma$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$, $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})={{\mbox{\boldmath $\beta$}}}{^{\mbox{\tiny {\sf T}}}}{{\bm W}}$ and ${{\bm W}}=(\tilde{{{\bm X}}},X_{1}^2,X_{2}^2,X_{1}X_{2})$. $d_0=5$, 10 or 15. The error terms $\epsilon_i$ follows i.i.d. N(0,1) or Gamma(1,1)-1 distribution. The propensity scores $\pi(\cdot)$ are known, and we consider both the constant case $\pi({{\bm X}}_i)=0.5$ and the non-constant case $\pi({{\bm X}}_i)=\mathrm{logit}({{\bm X}}_{i1}-{{\bm X}}_{i2})$. We report only the result of the constant case (Table \[table:interacted\_constant\_ps\]), and allocate the non-constant case to Appendix B.

  ---------------------------------------------------------- ------- ----- ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------
                                                                                                                                                                                                                                                                                    
  (r)[4-6]{} (lr)[7-9]{} (lr)[10-12]{} (lr)[13-15]{} Error   $d_0$   n       $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$   $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$   $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$   $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$
  Normal                                                     5       100                0.16              0.16               0.31              0.16              0.16               0.27              0.25              0.25               0.17              0.14              0.14               0.26
                                                                     200                0.09              0.09               0.24              0.10              0.10               0.19              0.18              0.18               0.09              0.08              0.08               0.19
                                                                     400                0.05              0.05               0.18              0.07              0.07               0.12              0.15              0.15               0.05              0.05              0.05               0.13
                                                                     800                0.02              0.02               0.14              0.05              0.05               0.09              0.14              0.14               0.04              0.03              0.03               0.09
                                                             10      100                0.28              0.28               0.92              0.22              0.22               0.81              0.39              0.39               0.40              0.21              0.21               0.82
                                                                     200                0.19              0.19               0.85              0.15              0.15               0.71              0.33              0.33               0.28              0.13              0.13               0.72
                                                                     400                0.12              0.12               0.79              0.10              0.10               0.60              0.30              0.30               0.23              0.09              0.09               0.63
                                                                     800                0.06              0.06               0.73              0.07              0.07               0.50              0.27              0.27               0.22              0.06              0.06               0.54
                                                             15      100                0.35              0.35               1.55              0.25              0.25               1.40              0.47              0.47               0.62              0.26              0.26               1.43
                                                                     200                0.27              0.27               1.48              0.18              0.18               1.31              0.45              0.45               0.45              0.18              0.18               1.34
                                                                     400                0.19              0.19               1.47              0.13              0.13               1.17              0.44              0.44               0.37              0.12              0.12               1.23
                                                                     800                0.12              0.12               1.39              0.09              0.09               1.03              0.41              0.41               0.35              0.08              0.08               1.07
  Gamma                                                      5       100                0.15              0.18               0.31              0.15              0.11               0.16              0.22              0.12               0.09              0.12              0.09               0.15
                                                                     200                0.09              0.12               0.26              0.10              0.06               0.10              0.18              0.08               0.05              0.08              0.05               0.09
                                                                     400                0.05              0.07               0.21              0.08              0.03               0.07              0.16              0.06               0.04              0.06              0.02               0.07
                                                                     800                0.02              0.04               0.17              0.07              0.03               0.06              0.15              0.06               0.03              0.05              0.02               0.07
                                                             10      100                0.26              0.33               0.90              0.22              0.16               0.54              0.39              0.13               0.27              0.22              0.14               0.50
                                                                     200                0.19              0.29               0.88              0.17              0.08               0.44              0.37              0.10               0.22              0.17              0.07               0.41
                                                                     400                0.12              0.24               0.87              0.13              0.04               0.39              0.35              0.08               0.19              0.14              0.03               0.36
                                                                     800                0.06              0.17               0.78              0.12              0.03               0.37              0.33              0.07               0.19              0.13              0.02               0.35
                                                             15      100                0.36              0.57               1.52              0.30              0.31               0.98              0.53              0.19               0.40              0.32              0.28               0.89
                                                                     200                0.28              0.53               1.51              0.22              0.19               0.81              0.55              0.16               0.29              0.24              0.16               0.71
                                                                     400                0.19              0.47               1.50              0.17              0.13               0.73              0.57              0.15               0.26              0.21              0.11               0.63
                                                                     800                0.11              0.43               1.50              0.15              0.11               0.71              0.58              0.15               0.24              0.18              0.09               0.62
  ---------------------------------------------------------- ------- ----- ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------

  : Summary results with constant propensity scores when errors interacted with treatment. Least square stands for lsA-learning. Pinball(0.5) stands for robust regression with pinball loss and parameter $\tau=0.5$. Pinball(0.25) stands for robust regression with pinball loss and parameter $\tau=0.25$. Huber stands for robust regression with Huber loss, where parameter $\alpha$ is tuned automatically with R function rlm.[]{data-label="table:interacted_constant_ps"}

We compare the performance of four methods: lsA-learning, robust regression with $\rho_{0.5}$ ($\mathrm{RR}(\rho_{0.5})$), robust regression with $\rho_{0.25}$ ($\mathrm{RR}(\rho_{0.25})$) and robust regression with Huber loss ($\mathrm{RR}(H)$). We consider four different sample sizes 100, 200, 400 and 800. For each scenario, we again simulate 1000 replications. When error terms are interactive with treatment, the true ${{\mbox{\boldmath $\beta$}}}_0$ associated with $g_{\mu}^{\mathrm{opt}}$ and $g_{\tau}^{\mathrm{opt}}$ are different. Specifically, under our model, ${{\mbox{\boldmath $\beta$}}}_0=({{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}},0,0,0){^{\mbox{\tiny {\sf T}}}}$ for $g_{\mu}^{\mathrm{opt}}$, ${{\mbox{\boldmath $\beta$}}}_0=({{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}},d_0 F^{-1}_{\epsilon}(0.5),0,0){^{\mbox{\tiny {\sf T}}}}$ for $g_{0.5}^{\mathrm{opt}}$ and ${{\mbox{\boldmath $\beta$}}}_0=({{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}},d_0 F^{-1}_{\epsilon}(0.25),0,0){^{\mbox{\tiny {\sf T}}}}$ for $g_{0.25}^{\mathrm{opt}}$. Thus, the two criteria, mse and PCD used in simulation study I, are no longer meaningful. So we evaluate the performance of methods in this simulation study based on value differences $\delta_{\mu}$, $\delta_{0.5}$ and $\delta_{0.25}$.

Based on Theorem \[thm:lsA\], we can prove that $\hat{g}^{A}_{LS}({{\bm x}})$ is consistent which converges to $g^{\mathrm{opt}}_{\mu}$ as sample size goes to infinity. This is shown in Table \[table:interacted\_constant\_ps\] such that the $\delta_{\mu}$ column for the lsA-learning method converges to 0 as sample size increases. We also know under Normal error terms, $\delta_{0.5}=\delta_{\mu}$. Thus, the $\delta_{0.5}$ column for the lsA-learning method also converges to 0. However, all other columns in Table \[table:interacted\_constant\_ps\] converge to a positive constant instead of 0 as sample size goes to infinity.

Another observation we discover from Table \[table:interacted\_constant\_ps\] is $\mathrm{RR}(H)$ and $\mathrm{RR}(\rho_{0.5})$ perform similarly. One additional observation we have is even though lsA-learning outperform all other methods in $\delta_{\mu}$ when sample size is large. It may be worse than $\mathrm{RR}(\rho_{0.5})$ and $\mathrm{RR}(H)$ when sample size is small. This is due to the fact that lsA-learning is inefficient under the heteroscedastic or skewed errors. The last observation we have is overall lsA-learning, $\mathrm{RR}(\rho_{0.5})$ and $\mathrm{RR}(\rho_{0.25})$ perform best at the columns $\delta_{\mu}$, $\delta_{0.5}$ and $\delta_{0.25}$ accordingly. The reason is given in the Remark under Theorem \[thm:approximation\_quantile\], which shows that $\hat{g}^{R}_{\rho(\tau)}$ $(\triangleq{\mathrm{1}}\{C({{\bm x}};\hat{{{\mbox{\boldmath $\beta$}}}}^R_{\rho(\tau)})>0\})$ in general approximates the unknown optimal ITR $g^{\mathrm{opt}}_{\tau}$ even when the conditional independence assumption $\epsilon\perp A|{{\bm X}}$ does not hold.

Application to AIDS study {#section:aids}
=========================

We illustrate the proposed robust regression method to data from AIDS Clinical Trials Group Protocol 175 (ACTG175), which has been previously studied by various authors [@leon2003semiparametric; @tsiatis2008covariate; @zhang2008improving; @lu2011variable]. In the study, 2139 HIV-infected subjects were randomized to four different treatment groups in equal proportions, and the treatment groups are zidovudine (ZDV) monotherapy, ZDV + didanosine (ddI), ZDV + zalcitabine, and ddI monotherapy. Following [@lu2011variable], we choose CD4 count $(\mathrm{cells/mm}^3)$ at $20\pm5$ weeks post-baseline as the primary continuous outcome $Y$, and include five continuous covariates and seven binary covariates as our covariates. They are: 1. age (years), 2. weight (kg), 3. karnof=Karnofsky score (scale of 0-100), 4. cd40=CD4 count $(\mathrm{cells/mm}^3)$ at baseline, 5. cd80=CD8 count $(\mathrm{cells/mm}^3)$ at baseline, 6. hemophilia=hemophilia (0=no, 1=yes), 7. homosexuality=homosexual activity (0=no, 1=yes), 8. drugs=history of intravenous drug use (0=no, 1=yes), 9. race (0=white, 1=non-white), 10. gender (0=female, 1=male), 11. str2= antiretroviral history (0=naive, 1=experienced), and 12. sympton=symptomatic status (0=asymptomatic, 1=symptomatic). For brevity, we only compare the treatment ZDV + didanosine (ddI) $(A=1)$ and ZDV + zalcitabine $(A=0)$, and restrict our samples to subjects receiving these two treatments. Thus, the propensity scores $\pi({{\bm X}}_i)\equiv 0.5$ in our restricted samples as the patients are assigned into one of two treatments with equal probability.

In our analysis, we assume linear models for both the baseline and the contrast functions. For interpretability, we keep the response $Y$ (the CD4 count) at its original scale, which is also consistent with the way clinicians think about the outcome in practice [@tsiatis2008covariate]. We plot the scatter plot of response Y against age. It shows some skewness and heterogeneity. With some preliminary analysis (fitting full model with lsA-learning and RR(M)), we find that only covariates age, homosexuality and race may possibly interact with the treatment. So in our final model, only these three covariates are included in the contrast function, while at the same time we still keep all twelve covariates in the baseline function. The estimated coefficients associated with their corresponding standard errors and p-values are given in Table \[table:aids\], where standard errors are estimated with 1000 bootstrap samples (parametric bootstrap) and p-values are calculated with normal approximation. Only coefficients included in the contrast function are shown.

  ---------------------------------------------------------- -------- ------- ----------- -------- ------- ----------- -------- ------- ----------- -------- ------- -----------
                                                                                                                                                                     
  (r)[2-4]{} (lr)[5-7]{} (l)[8-10]{} (l)[11-13]{} Variable       Est.      SE          PV     Est.      SE          PV     Est.      SE          PV     Est.      SE          PV
  intercept                                                    -42.61   32.93       0.196   -33.45   37.32       0.370   -35.77   39.17       0.361   -42.76   31.40       0.173
  age                                                            3.13    0.85   **0.000**     2.62    0.97   **0.007**     2.46    1.06   **0.020**     2.80    0.79   **0.000**
  homosexuality                                                -40.66   16.73   **0.015**   -33.18   17.68       0.061   -35.38   18.28       0.053   -27.33   15.19       0.072
  race                                                         -25.70   17.69       0.146   -33.56   18.12       0.064   -34.21   18.32       0.062   -25.29   16.08       0.116
  ---------------------------------------------------------- -------- ------- ----------- -------- ------- ----------- -------- ------- ----------- -------- ------- -----------

  : Analysis results for AIDS data. Est. stands for estimate; SE stands for standard error; PV stands for p-value. All p-values which are significant at level 0.1 are highlighted.[]{data-label="table:aids"}

From Tables \[table:aids\], we make the following observations. First, lsA-learning (equivalent to Q- and A-learning with this model setting) and robust regression with pinball loss and Huber loss all have estimates with the exact same signs. Second, the estimated coefficients are distinguishable across different methods. Third, the covairiate homosexuality is significant under lsA-learning, but it is not significant under robust regression with either pinball losses or Huber loss, when the significant level $\alpha$ is set to 0.05.

We could further estimate the values $(V_{\mu}(\hat{g}))$ associated with each method by either the inverse probability weighted estimator (IPWE) [@robin2000marginal] or the augmented inverse probability weighted estimator (AIPWE) [@robins1994estimation], where $$\begin{aligned}
\hat{V}^{\mathrm{IPWE}}_{\mu}(\hat{g})=&\frac{\sum_{i=1}^n{\mathrm{1}}{\{A_i=\hat{g}({{\bm X}}_i)\}}Y_i/p(A_i|{{\bm X}}_i)}
{\sum_{i=1}^n{\mathrm{1}}{\{A_i=\hat{g}({{\bm X}}_i)\}}/p(A_i|{{\bm X}}_i)},\\
\hat{V}^{\mathrm{AIPWE}}_{\mu}(\hat{g})=&\frac{1}{n}\sum_{i=1}^{n}\hat{{\mbox{E}}}(Y_i|{{\bm X}}_i,\hat{g}({{\bm X}}_i))
+\frac{1}{n}\sum_{i=1}^{n}\frac{{\mathrm{1}}{\{A_i=\hat{g}({{\bm X}}_i)\}}}{p(A_i|{{\bm X}}_i)}\left[Y_i-\hat{{\mbox{E}}}(Y_i|{{\bm X}}_i,A_i)\right],\end{aligned}$$ $\hat{{\mbox{E}}}(Y_i|{{\bm X}}_i,A_i))=\varphi({{\bm X}}_i;\hat{{{\mbox{\boldmath $\gamma$}}}})+\left\{A_i-p(A_i|{{\bm X}}_i)\right\}C({{\bm X}}_i;\hat{{{\mbox{\boldmath $\beta$}}}})$, and $p(A_i|{{\bm X}}_i)\equiv0.5$. Both $\hat{V}^{\mathrm{IPWE}}_{\mu}(\hat{g})$ and $\hat{V}^{\mathrm{AIPWE}}_{\mu}(\hat{g})$ are consistent estimator of value $V_{\mu}(\hat{g})$, and their asymptotic covariance matrix can also be consistently estimated from the data [@zhang2012robust; @mckeague2014estimation]. The estimates of $(V_{\mu}(\hat{g}))$ and their corresponding 95% confidence interval of four methods based on both IPWE and AIPWE are given in Table \[table:aids\_value\].

  Estimator   method           Value     SE           CI
  ----------- --------------- -------- ------ ------------------
  IPWE        Least Square     405.05   6.72   (391.88, 418.22)
              Pinball(0.5)     406.77   6.71   (393.63, 419.92)
              Pinball(0.25)    406.07   6.73   (392.87, 419.26)
              Huber            407.03   6.71   (393.87, 420.18)
  AIPWE       Least Square     404.39   6.12   (392.40, 416.38)
              Pinball(0.5)     405.93   6.13   (393.92, 417.94)
              Pinball(0.25)    403.60   6.62   (390.62, 416.58)
              Huber            406.00   6.15   (393.95, 418.04)

  : Result of estimated values and their corresponding 95% confidence interval for four methods based on IPWE and AIPWE.SE stands for standard error. CI stands for 95% confidence interval.[]{data-label="table:aids_value"}

From Table \[table:aids\_value\], robust regression with $\rho_{0.5}$ and Huber loss perform slightly better than lsA-learning, while robust regression with $\rho_{0.25}$ performs worse than lsA-learning when the values $(V_{\mu}(\hat{g}))$ is estimated based on AIPWE. We conduct KCI-test to check the conditional independence assumption $\epsilon\perp A|{{\bm X}}$. For $\mathrm{RR}(\rho(0.5))$, $\mathrm{RR}(\rho(0.25))$ and RR(H), their p-values associated with KCI-test are 0.060, 0.002 and 0.083 respectively. The conditional independence assumption holds at the significance level of 0.05 for $\mathrm{RR}(\rho(0.5))$ and RR(H), so the estimated ITR can be thought to maximize $V_{\mu}(g)$. On the other hand, this assumption doesn’t hold for $\mathrm{RR}(\rho(0.25))$, and its estimated ITR doesn’t maximize $V_{\mu}(g)$, instead it approximately maximizes $V_{0.25-q}(g)$. This partly explains the relatively bad performance of RR($\rho_{0.25}$) in Table \[table:aids\_value\]. Again, as $\mathrm{RR}(\rho(0.5))$ and RR(H) are more robust against heterogeneous, right skewed errors comparing with the least square method, they slightly outperform lsA-learning in term of $V_{\mu}(g)$.

Discussion
==========

In this article, we propose a new general loss based robust regression framework for estimating the optimal individualized treatment rules. This new method has the desired property to be robust against skewed, heterogeneous, heavy-tailed errors and outliers. And similar as A-learning, it produces consistent estimates of the optimal ITR even when the baseline function is misspecified. However, the consistency of the proposed method does require the key conditional independence assumption $\epsilon\perp A|{{\bm X}}$, which is somewhat stronger than the condition needed for the consistency of Q- and A-learning $({\mbox{E}}(\epsilon|{{\bm X}},A)=0)$. So there are situations when the classical Q- and A-learning are more appropriate to apply. Furthermore, we also point out in the article that when pinball loss $\rho_{\tau}$ is chosen and the assumption $\epsilon\perp A|{{\bm X}}$ doesn’t hold, the estimated ITR approximately maximize the conditional $\tau$-th quantile and thus maximize $V_{\tau-q}(g)$. From a practice point of view, there are situations when maximizing $V_{\tau-q}(g)$ is a much more reasonable approach comparing with maximizing $V_{\mu}(g)$, especially when the conditional distribution of response $Y$ is highly skewed to one side.

In practice, there are cases when multiple treatment groups need to be compared simultaneously. For brevity, we have limited our discussion to two treatment groups. However, the proposed method can be readily extended to multiple cases by just replacing equation with the following more complex form, $$L_{3n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})=\frac{1}{n}\sum_{i=1}^{n}M\left[Y_i-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})-\sum_{k=1}^{K-1}(I(A_i=k)-\pi_k({{\bm X}}_i))C_k({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}}_k)\right],
\label{eq:A-general-loss2}$$ where $\mathcal{A}=\{1,\ldots,K\}$, $K$-th treatment is the baseline treatment, $\pi_k({{\bm X}}_i)=\Pr(A_i=k|{{\bm X}}_i)$ and $C_k({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}}_k)$ denotes the contrast function comparing $k$-th treatment and the baseline treatment. All Theorems can be easily extended to this multiple treatments setting as well.

When the dimension of prognostic variables is high, regularized regression is needed in order to produce parsimonious yet interpretable individualized treatment rules. Essentially this is a variable selection problem in the context of M-estimator, which has been previously studied in [-@wu2009variable; -@li2011nonconcave], etc. This is an interesting topic that needs further investigation. Another interesting direction is to extend the current method to the multi-stage setting, where sequential decisions are made along the time line.

Appendix A: Proof of Asymptotic Properties {#appendix-a-proof-of-asymptotic-properties .unnumbered}
==========================================

We consider the following additive model, $$Y_i=\varphi_0({{\bm X}}_i)+\{A_i-\pi({{\bm X}}_i)\}C({{\bm X}}_i;\beta_0)+\epsilon_i,\; i=1,\ldots,n,$$ where $\varphi_0({{\bm X}})$ is the baseline function, $C({{\bm X}};\beta_0)$ is the contrast function, $\pi({{\bm X}})$ is the propensity score, and $\epsilon$ is the error term. We estimate $({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})$ by minimizing $$L_{3n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})=\frac{1}{n}\sum_{i=1}^{n}M\left[Y_i-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})-\{A_i-\pi({{\bm X}}_i)\}C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})\right],
\label{eq:A-general-loss}$$ where ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$, ${{\mbox{\boldmath $\beta$}}}\in\mathcal{B}$ and $M:{\mathrm{I \! R} \mathit{^{\rightarrow}}} [0,\infty)$ is a convex function with minimum achieved at 0. We consider the following three types of loss functions, i.e., the pinball loss $$M(x)=\rho_\tau(x)\triangleq
\begin{cases}
(\tau-1)x, &\text{if } x<0\\
\tau x, &\text{if } x\geq0
\end{cases}$$ where $0<\tau<1$, the Huber loss $$M(x)=H_\alpha(x)\triangleq
\begin{cases}
0.5x^2, &\text{if } |x|<\alpha\\
\alpha|x|-0.5\alpha^2, &\text{if } |x|\geq\alpha
\end{cases}$$ for some $\alpha>0$, and the $\epsilon$-insensitive loss $$M(x)=J_\epsilon(x)\triangleq\max(0, |x|-\epsilon)$$ for some $\epsilon>0$. Define $\Delta C({{\bm x}};{{\mbox{\boldmath $\beta$}}})=C({{\bm x}};{{\mbox{\boldmath $\beta$}}})-C({{\bm x}};{{\mbox{\boldmath $\beta$}}}_0)$. Assume ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$, ${{\mbox{\boldmath $\beta$}}}\in\mathcal{B}$ and ${{\mbox{\boldmath $\gamma$}}}'$ is any arbitrary fix point in $\Gamma$.

**Regularity conditions A:**

-   $\{(Y_i,{{\bm X}}_i,A_i,\epsilon_i),i=1,\ldots,n\}$ are i.i.d random variables.

-   $\epsilon_i\perp A_i|{{\bm X}}_i$ $\forall i=1,\ldots,n$.

-   ${\mbox{E}}|\Delta C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})|<\infty$ $\forall{{\mbox{\boldmath $\beta$}}}\in\mathcal{B}$.

-   $\Pr\{{{\bm x}}\in\mathcal{X}:\;\Delta C({{\bm x}};{{\mbox{\boldmath $\beta$}}})\neq 0\}>0$ for all ${{\mbox{\boldmath $\beta$}}}\neq{{\mbox{\boldmath $\beta$}}}_0$.

-   ${\mbox{E}}|\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})|<\infty$ $\forall{{\mbox{\boldmath $\gamma$}}}\in\Gamma$.

-   $G_2({{\mbox{\boldmath $\gamma$}}})$ has unique minimizer ${{\mbox{\boldmath $\gamma$}}}^*$, where $G_2({{\mbox{\boldmath $\gamma$}}})$ is the pointwise limit of $L_{3n}({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}})-L_{3n}({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}}')$ in probability.

-   $L_{3n}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})$ is strictly convex with respect to $({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})$.

-   $\epsilon|{{\bm X}}={{\bm x}}$ has nonzero density on $\mathbb{R}$ for almost all ${{\bm x}}\in\mathcal{X}$.

$\left|\rho_\tau(x-y)-\rho_\tau(x)\right|\leq |y|$, for all $\tau\in(0,1)$.

$$\begin{aligned}
    \left|\rho_\tau(x-y)-\rho_\tau(x)\right| &= \left|\tau\left\{(x-y)_{+}-x_{+}\right\}+(1-\tau)\left\{(x-y)_{-}-x_{-}\right\}\right|\\
    &\leq|(x-y)_{+}-x_{+}|+|(x-y)_{-}-x_{-}|=|y|
    \end{aligned}$$

$$\begin{aligned}
    \rho_\tau(x-y)-\rho_\tau(x)=&-\tau y{\mathrm{1}}\{x\geq 0\}+(1-\tau)y{\mathrm{1}}\{x< 0\}+(y-x){\mathrm{1}}\{x\geq 0\}{\mathrm{1}}\{y>x\}\\
    &+(x-y){\mathrm{1}}\{x< 0\}{\mathrm{1}}\{y< x\},
    \end{aligned}$$

for all $\tau\in(0,1)$.

Denote $D=\rho_\tau(x-y)-\rho_\tau(x)$.

1.  If $x\geq0$, $y\leq0$ $\Rightarrow$ $D=-\tau y$;

2.  If $x\geq0$, $y>0$, $|x|\geq|y|$ $\Rightarrow$ $D=-\tau y$;

3.  If $x\geq0$, $y>0$, $|x|<|y|$ $\Rightarrow$ $D=-\tau y+(y-x)$;

4.  If $x<0$, $y\geq0$ $\Rightarrow$ $D=(1-\tau)y$;

5.  If $x<0$, $y<0$, $|x|\geq|y|$ $\Rightarrow$ $D=(1-\tau)y$;

6.  If $x<0$, $y<0$, $|x|<|y|$ $\Rightarrow$ $D=(1-\tau)y+(x-y)$;

Combining the above 6 cases, Lemma 2 is proved.

**Proof of Theorem 1.**

Recall that the loss function defined in takes the form $$L_{3n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})=\frac{1}{n}\sum_{i=1}^{n}\rho_\tau\left[\varphi_0({{\bm X}}_i)-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})+\epsilon_i
    -(A_i-\pi({{\bm X}}_i))\Delta C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})\right].$$ By definition, $$\begin{aligned}
    (\hat{{{\mbox{\boldmath $\beta$}}}}^{R}_{\rho(\tau)},\hat{{{\mbox{\boldmath $\gamma$}}}}^{R}_{\rho(\tau)}) =& {\mathrm{argmin}}_{({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})}L_{3n}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})-L_{3n}({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}}')\\
    =& {\mathrm{argmin}}_{({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})}\left[L_{3n}({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})-L_{3n}({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}})\right]+
    \left[L_{3n}({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}})-L_{3n}({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}}')\right],
    \end{aligned}$$ Define $$\begin{aligned}
    S_{1n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}}) =& L_{3n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})-L_{3n}({{\mbox{\boldmath $\beta$}}}_0, {{\mbox{\boldmath $\gamma$}}})=1/n\sum_{i=1}^n d_{1i};\\
    S_{2n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}}) =& L_{3n}({{\mbox{\boldmath $\beta$}}}_0, {{\mbox{\boldmath $\gamma$}}})-L_{3n}({{\mbox{\boldmath $\beta$}}}_0, {{\mbox{\boldmath $\gamma$}}}')=1/n\sum_{i=1}^n d_{2i}
    \end{aligned}$$ where $$\begin{aligned}
    d_{1i} =& \rho_\tau\left[\varphi_0({{\bm X}}_i)-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})+\epsilon_i
    -(A_i-\pi({{\bm X}}_i))\Delta C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})\right]-\rho_\tau\left[\varphi_0({{\bm X}}_i)-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})+\epsilon_i\right],\\
    d_{2i} =& \rho_\tau\left[\varphi_0({{\bm X}}_i)-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})+\epsilon_i\right]-\rho_\tau\left[\varphi_0({{\bm X}}_i)-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}}')+\epsilon_i\right].
    \end{aligned}$$ By Lemma 1, A3 and A5, ${\mbox{E}}|d_{1i}|\leq{\mbox{E}}|(A_i-\pi({{\bm X}}_i))\Delta C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})|\leq{\mbox{E}}|\Delta C({{\bm X}}_i;{{\mbox{\boldmath $\beta$}}})|<\infty$ and ${\mbox{E}}|d_{2i}|\leq{\mbox{E}}|\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})-\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}}')|\leq{\mbox{E}}|\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}})|+{\mbox{E}}|\varphi({{\bm X}}_i;{{\mbox{\boldmath $\gamma$}}}')|
    <\infty$. Then, by Law of Large Number, $\forall\;{{\mbox{\boldmath $\beta$}}}\in\mathcal{B}$, ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$, we have $S_{1n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})\inprob G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\triangleq{\mbox{E}}(D)$, and $S_{2n}({{\mbox{\boldmath $\beta$}}}, {{\mbox{\boldmath $\gamma$}}})\inprob G_2({{\mbox{\boldmath $\gamma$}}})$, where $$\begin{aligned}
    D=&\rho_\tau\left[\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\epsilon
    -\{A-\pi({{\bm X}})\}\Delta C({{\bm X}};{{\mbox{\boldmath $\beta$}}})\right]-\rho_\tau\left[\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\epsilon\right].
    \end{aligned}$$ Below we show that a) $({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}}^*)$ is the minimizer of $G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})+G_2({{\mbox{\boldmath $\gamma$}}})$, b) $({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}}^*)$ is the unique minimizer. The consistency then follows from the argmax continuous mapping theorem under Assumption (A7).

Denote $K_1=\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\epsilon$, $K_2=\{A-\pi({{\bm X}})\}\Delta C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$. By Lemma 2, $$\begin{aligned}
    D =& -\tau K_2{\mathrm{1}}\{K_1\geq0\}+(1-\tau) K_2{\mathrm{1}}\{K_1<0\}+(K_2-K_1){\mathrm{1}}\{K_1\geq0\}{\mathrm{1}}\{K_2>K_1\}\\
    &+(K_1-K_2){\mathrm{1}}\{K_1<0\}{\mathrm{1}}\{K_2<K_1\}.
    \end{aligned}$$ Since $\epsilon\perp A|{{\bm X}}$ and $\Pr(A|{{\bm X}})=\pi({{\bm X}})$, applying double expectation rule with ${{\bm X}}$, we have ${\mbox{E}}[-\tau K_2{\mathrm{1}}\{K_1\geq0\}]={\mbox{E}}[(1-\tau) K_2{\mathrm{1}}\{K_1<0\}]=0$. Thus, $$G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})={\mbox{E}}[(K_2-K_1){\mathrm{1}}\{K_1\geq0\}{\mathrm{1}}\{K_2>K_1\}]+{\mbox{E}}[(K_1-K_2){\mathrm{1}}\{K_1<0\}{\mathrm{1}}\{K_2<K_1\}].
    \label{eq:Gfunction}$$ It is easy to check $G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\geq0$ and achieves minimal value 0 at point $({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}})$ for all ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$. In addition, by A6, we know $G_2({{\mbox{\boldmath $\gamma$}}})$ has unique minimizer ${{\mbox{\boldmath $\gamma$}}}^*$. Combining the above two facts, a) is proved.

Combining A4, A8 and , we could prove $G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})>0$ for all ${{\mbox{\boldmath $\beta$}}}\neq{{\mbox{\boldmath $\beta$}}}_0$ and ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$. So b) holds.

**Proof of Theorem 3.**

\(a) When $M(x)=H_{\alpha}(x)$, the proof follows similar steps as Theorem 1. The only difference is that $G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})$ takes a different expression now and we need to redo the proof of 1) $G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})>0$ $\forall{{\mbox{\boldmath $\beta$}}}\neq{{\mbox{\boldmath $\beta$}}}_0$, ${{\mbox{\boldmath $\gamma$}}}\in\Gamma$, and 2) $G_1({{\mbox{\boldmath $\beta$}}}_0,{{\mbox{\boldmath $\gamma$}}})=0$ $\forall{{\mbox{\boldmath $\gamma$}}}\in\Gamma$. By definition, $G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\triangleq{\mbox{E}}(D)$, where $$\begin{aligned}
    D=H_\alpha\left[\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\epsilon
    -\{A-\pi({{\bm X}})\}\Delta C({{\bm X}};{{\mbox{\boldmath $\beta$}}})\right]-H_\alpha\left[\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\epsilon\right].
    \end{aligned}$$ Then, 2) holds immediately. Denote $K_1=\varphi_0({{\bm X}})-\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})+\epsilon$, $K_2=\{A-\pi({{\bm X}})\}\Delta C({{\bm X}};{{\mbox{\boldmath $\beta$}}})$. We have the following four cases:

1.  If $K_1>\alpha$ then $H_\alpha(K_1-K_2)\geq\alpha(K_1-K_2)-0.5\alpha^2$. Thus, $D\geq-\alpha K_2$;

2.  If $K_1<-\alpha$ then $H_\alpha(K_1-K_2)\geq\alpha(K_2-K_1)-0.5\alpha^2$. Thus, $D\geq\alpha K_2$;

3.  If $K_1\in[-\alpha,\alpha]$ and $K_1-K_2\in[-\alpha,\alpha]$ then $D=1/2(K_1-K_2)^2-1/2K_1^2=-K_1K_2+1/2K_2^2$;

4.  If $K_1\in[-\alpha,\alpha]$ and $K_1-K_2\not\in[-\alpha,\alpha]$ then $H_\alpha(K_1-K_2)\geq1/2(K_1-K_2)^2-\left\{1/2(\alpha+|K_2|)^2-\left[\alpha(\alpha+|K_2|)-1/2\alpha^2\right]\right\}=1/2(K_1-K_2)^2-1/2K_1^2$. Thus, $D\geq1/2(K_1-K_2)^2-1/2K_1^2-1/2K_2^2=-K_1K_2$.

Combining the above four equalities and inequalities, $$\begin{aligned}
    G_1({{\mbox{\boldmath $\beta$}}},{{\mbox{\boldmath $\gamma$}}})\geq& {\mbox{E}}[-\alpha K_2{\mathrm{1}}\{K_1>\alpha\}] + {\mbox{E}}[\alpha K_2{\mathrm{1}}\{K_1<-\alpha\}] +
    {\mbox{E}}[-K_1K_2{\mathrm{1}}\{K_1\in[-\alpha,\alpha]\}]\\
    &+ {\mbox{E}}\left[1/2K_2^2{\mathrm{1}}\left(\{K_1\in[-\alpha,\alpha]\}\cup\{K_1-K_2\in[-\alpha,\alpha]\}\right)\right]
    \end{aligned}$$ Since $\epsilon\perp A|{{\bm X}}$ and $\Pr(A|{{\bm X}})=\pi({{\bm X}})$, applying double expectation rule with ${{\bm X}}$, we have ${\mbox{E}}[-\alpha K_2{\mathrm{1}}\{K_1>\alpha\}]={\mbox{E}}[\alpha K_2{\mathrm{1}}\{K_1<-\alpha\}]={\mbox{E}}[-K_1K_2{\mathrm{1}}\{K_1\in[-\alpha,\alpha]\}]=0$. Thus, $$G_1({{\mbox{\boldmath $\beta$}}};{{\mbox{\boldmath $\gamma$}}})\geq{\mbox{E}}\left[1/2K_2^2{\mathrm{1}}\left(\{K_1\in[-\alpha,\alpha]\}\cup\{K_1-K_2\in[-\alpha,\alpha]\}\right)\right].
    \label{eq:Gfunction_thm2}$$ Combining , A4 and A8, we can check that 1) holds. Thus, part (a) is proved.

\(b) When $M(x)=J_{\epsilon}(x)$, similarly $D=J_\epsilon\left(K_1-K_2\right)-J_\epsilon\left(K_1\right)$. Notice that we have the following three cases:

1.  If $K_1>\epsilon$ then $D\geq -K_2$;

2.  If $K_1<-\epsilon$ then $D\geq K_2$;

3.  If $K_1\in[-\epsilon,\epsilon]$ then $D\geq 0$;

The rest of the proof follows similar steps as part (a).

**Proof of Theorem 5.**

From Theorem 1, ${{\mbox{\boldmath $\beta$}}}_{\tau}={{\mbox{\boldmath $\beta$}}}_0$. Plugging this into Theorem 4 and applying double expectation rules, we have $$J(\tau)={\mbox{E}}\left[f_{\epsilon}\left(\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})|{{\bm X}}\right)
    \left(\begin{array}{cc}
    \pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}& {{\mathbf 0}}\\
    {{\mathbf 0}}& \tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\end{array}
    \right)\right]$$ and $$\Sigma(\tau,\tau)={\mbox{E}}\left\{\left[\tau-{\mathrm{1}}\left\{\epsilon<\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})\right\}\right]^2
    \left(\begin{array}{cc}
    \pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}& {{\mathbf 0}}\\
    {{\mathbf 0}}& \tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\end{array}
    \right)\right\}.$$ Thus, $\sqrt{n}(\hat{{{\mbox{\boldmath $\beta$}}}}(\tau)-{{\mbox{\boldmath $\beta$}}}_0)\indist N({{\mathbf 0}}, J_{11}^{-1}(\tau)\Sigma_{11}(\tau,\tau)J_{11}^{-1}(\tau))$, where $J_{11}^{-1}(\tau)$ and $\Sigma_{11}(\tau,\tau)$ are defined as in Theorem 5. Conditional on ${{\bm X}}$, ${\mathrm{1}}\left\{\epsilon<\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})\right\}$ is a binomial random variable with $p=\Pr\left(\epsilon<\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})\right)$. Then, ${\mbox{E}}\left\{\left[\tau-{\mathrm{1}}\{\epsilon<\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}{{\mbox{\boldmath $\gamma$}}}(\tau)-\varphi_0({{\bm X}})\}\right]^2|{{\bm X}}\right\}=(p-\tau)^2+p(1-p)\leq \tau^2+|1-2\tau|$. Thus, $\Sigma_{11}(\tau,\tau)\leq\left(\tau^2+|1-2\tau|\right){\mbox{E}}\left[\pi({{\bm X}})\{1-\pi({{\bm X}})\}\tilde{{{\bm X}}}\tilde{{{\bm X}}}{^{\mbox{\tiny {\sf T}}}}\right]$.

Appendix B: Additional Simulation Results {#appendix-b-additional-simulation-results .unnumbered}
=========================================

We conducted additional simulations with non-constant propensity scores. Specifically, we considered the following examples.

Examples with error terms independent with treatment {#examples-with-error-terms-independent-with-treatment .unnumbered}
----------------------------------------------------

We consider the following two models with p=3,

-   Model I: $$Y_i=1+(X_{i1}-X_{i2})(X_{i1}+X_{i3})+\{A_i-\pi({{\bm X}}_{i})\}{{\mbox{\boldmath $\beta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\sigma({{\bm X}}_{i})\epsilon_i,$$ where ${{\bm X}}_{i}=(X_{i1},X_{i2},X_{i3}){^{\mbox{\tiny {\sf T}}}}$ are multivariate normal with mean 0, variance 1, and $\mathrm{Corr}(X_{ij},X_{ik})=0.5^{|j-k|}$, $\tilde{{{\bm X}}}_i=(1,{{\bm X}}_i{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$ and ${{\mbox{\boldmath $\beta$}}}_0=(0,1,-1,1){^{\mbox{\tiny {\sf T}}}}$.

-   Model II: $$Y_i={{\mbox{\boldmath $\gamma$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\{A_i-\pi({{\bm X}}_{i})\}{{\mbox{\boldmath $\beta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\sigma({{\bm X}}_{i})\epsilon_i,$$ where ${{\mbox{\boldmath $\gamma$}}}_0{^{\mbox{\tiny {\sf T}}}}=(0.5,4,1,-3)$, and ${{\bm X}}_{i}$, $\tilde{{{\bm X}}}_{i}$ and ${{\mbox{\boldmath $\beta$}}}_0$ are the same as Model I.

We take linear forms for both the baseline and the contrast functions, where $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})={{\mbox{\boldmath $\gamma$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$ and $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})={{\mbox{\boldmath $\beta$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$. We assume the propensity scores $\pi(\cdot)$ are known, and we study the non-constant case $(\pi({{\bm X}}_i)=\mathrm{logit}({{\bm X}}_{i1}-{{\bm X}}_{i2}))$ here. In addition, We consider two different $\sigma({{\bm X}}_{i})$ functions, i.e., the homogeneous case with $\sigma({{\bm X}}_{i})=1$, and the heterogenous case with $\sigma({{\bm X}}_{i})=0.5+(X_{i1}-X_{i2})^2$. The simulation results are given in Table \[table:modelI\_nonconstant\] and Table \[table:modelII\_nonconstant\].

[@ll ccc ccc ccc]{}\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 1.70 (0.061) & 81.9 & 0.91 & 2.90 (0.114) & 77.6 & 1.34 & & 59.3 & 3.61\
& P(0.5) & 1.90 (0.069) & 80.1 & 1.09 & 2.13 (0.073) & 78.3 & 1.25 & 3.54 (0.128) & 75.7 & 1.57\
& P(0.25) & 2.35 (0.080) & 78.2 & 1.33 & 1.95 (0.076) & 80.4 & 1.08 & 8.45 (0.431) & 69.8 & 2.28\
& Huber & 1.51 (0.053) & 82.1 & 0.89 & 1.77 (0.065) & 80.6 & 1.02 & 3.67 (0.127) & 75.4 & 1.60\
200 & LS & 0.77 (0.026) & 86.8 & 0.50 & 1.35 (0.045) & 82.2 & 0.91 & & 59.2 & 3.63\
& P(0.5) & 0.88 (0.028) & 85.5 & 0.60 & 1.00 (0.029) & 83.0 & 0.79 & 1.54 (0.050) & 81.1 & 1.00\
& P(0.25) & 1.06 (0.035) & 84.5 & 0.68 & 0.83 (0.027) & 85.9 & 0.59 & 3.61 (0.143) & 74.7 & 1.70\
& Huber & 0.68 (0.022) & 87.3 & 0.46 & 0.81 (0.025) & 85.2 & 0.62 & 1.58 (0.052) & 80.7 & 1.03\
400 & LS & 0.39 (0.012) & 90.2 & 0.28 & 0.65 (0.020) & 86.9 & 0.48 & & 58.0 & 3.79\
& P(0.5) & 0.43 (0.013) & 89.3 & 0.32 & 0.47 (0.014) & 88.4 & 0.38 & 0.73 (0.022) & 86.5 & 0.51\
& P(0.25) & 0.53 (0.016) & 88.5 & 0.38 & 0.41 (0.013) & 90.5 & 0.27 & 1.50 (0.049) & 81.7 & 0.96\
& Huber & 0.34 (0.010) & 90.6 & 0.25 & 0.39 (0.012) & 89.6 & 0.30 & 0.72 (0.022) & 86.3 & 0.53\
800 & LS & 0.18 (0.006) & 93.3 & 0.13 & 0.32 (0.010) & 90.2 & 0.27 & & 58.3 & 3.75\
& P(0.5) & 0.21 (0.007) & 92.7 & 0.15 & 0.24 (0.007) & 91.5 & 0.20 & 0.36 (0.011) & 90.3 & 0.27\
& P(0.25) & 0.28 (0.009) & 92.4 & 0.17 & 0.21 (0.007) & 93.4 & 0.13 & 0.78 (0.026) & 86.9 & 0.50\
& Huber & 0.16 (0.005) & 93.7 & 0.11 & 0.19 (0.006) & 92.6 & 0.15 & 0.37 (0.010) & 89.9 & 0.28\
\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 2.84 (0.111) & 78.2 & 1.33 & 9.96 (0.773) & 72.0 & 2.06 & & 55.2 & 4.18\
& P(0.5) & 2.01 (0.082) & 80.6 & 1.09 & 2.18 (0.080) & 79.2 & 1.21 & 4.18 (0.189) & 74.1 & 1.81\
& P(0.25) & 2.91 (0.110) & 76.7 & 1.52 & 3.22 (0.105) & 74.2 & 1.76 & 10.62 (0.475) & 65.3 & 2.87\
& Huber & 1.90 (0.074) & 80.9 & 1.06 & 2.38 (0.090) & 78.1 & 1.32 & 5.06 (0.230) & 71.9 & 2.04\
200 & LS & 1.46 (0.053) & 83.1 & 0.83 & 4.47 (0.371) & 76.8 & 1.51 & & 56.3 & 4.04\
& P(0.5) & 0.92 (0.033) & 86.4 & 0.55 & 0.98 (0.035) & 85.3 & 0.64 & 1.69 (0.065) & 81.5 & 0.98\
& P(0.25) & 1.35 (0.049) & 83.3 & 0.81 & 1.47 (0.049) & 81.6 & 0.97 & 4.73 (0.241) & 71.9 & 2.05\
& Huber & 0.86 (0.030) & 86.6 & 0.52 & 1.02 (0.036) & 84.7 & 0.68 & 2.11 (0.079) & 79.3 & 1.18\
400 & LS & 0.74 (0.029) & 87.4 & 0.47 & 2.65 (0.402) & 81.4 & 1.04 & & 56.2 & 4.06\
& P(0.5) & 0.45 (0.016) & 90.2 & 0.29 & 0.44 (0.017) & 89.5 & 0.34 & 0.79 (0.029) & 87.2 & 0.49\
& P(0.25) & 0.66 (0.025) & 88.3 & 0.41 & 0.70 (0.023) & 86.9 & 0.50 & 2.12 (0.091) & 79.5 & 1.19\
& Huber & 0.43 (0.016) & 90.2 & 0.28 & 0.48 (0.018) & 89.0 & 0.36 & 1.01 (0.036) & 85.0 & 0.65\
800 & LS & 0.36 (0.013) & 90.8 & 0.25 & 1.09 (0.066) & 85.0 & 0.69 & & 56.3 & 4.02\
& P(0.5) & 0.21 (0.008) & 93.2 & 0.14 & 0.24 (0.009) & 92.3 & 0.19 & 0.39 (0.014) & 90.5 & 0.27\
& P(0.25) & 0.33 (0.013) & 91.7 & 0.21 & 0.36 (0.012) & 90.8 & 0.25 & 1.01 (0.034) & 84.9 & 0.65\
& Huber & 0.20 (0.008) & 93.2 & 0.14 & 0.25 (0.009) & 92.1 & 0.19 & 0.49 (0.016) & 89.1 & 0.34\

[@ll ccc ccc ccc]{}\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 0.36 (0.011) & 89.8 & 0.29 & 1.65 (0.085) & 80.8 & 1.06 & & 58.7 & 3.69\
& P(0.5) & 0.57 (0.017) & 86.9 & 0.46 & 0.61 (0.026) & 86.4 & 0.55 & 1.31 (0.045) & 81.7 & 0.93\
& P(0.25) & 0.65 (0.020) & 86.2 & 0.52 & 0.22 (0.008) & 91.7 & 0.20 & 4.67 (0.312) & 74.7 & 1.64\
& Huber & 0.38 (0.012) & 89.5 & 0.30 & 0.45 (0.018) & 88.3 & 0.40 & 1.70 (0.060) & 79.5 & 1.14\
200 & LS & 0.16 (0.004) & 92.9 & 0.14 & 0.74 (0.030) & 85.6 & 0.61 & & 59.1 & 3.64\
& P(0.5) & 0.25 (0.007) & 91.2 & 0.21 & 0.26 (0.008) & 90.7 & 0.24 & 0.52 (0.017) & 87.8 & 0.41\
& P(0.25) & 0.30 (0.008) & 90.3 & 0.26 & 0.09 (0.003) & 94.8 & 0.08 & 1.69 (0.074) & 81.3 & 0.92\
& Huber & 0.17 (0.005) & 92.8 & 0.14 & 0.19 (0.006) & 92.2 & 0.17 & 0.70 (0.022) & 86.2 & 0.53\
400 & LS & 0.08 (0.002) & 95.1 & 0.06 & 0.36 (0.013) & 89.7 & 0.30 & & 58.0 & 3.79\
& P(0.5) & 0.12 (0.003) & 93.8 & 0.10 & 0.12 (0.003) & 93.8 & 0.10 & 0.22 (0.006) & 91.6 & 0.19\
& P(0.25) & 0.14 (0.004) & 93.3 & 0.12 & 0.04 (0.001) & 96.5 & 0.03 & 0.63 (0.021) & 86.5 & 0.49\
& Huber & 0.08 (0.002) & 95.0 & 0.07 & 0.09 (0.002) & 94.8 & 0.07 & 0.30 (0.009) & 90.3 & 0.26\
800 & LS & 0.04 (0.001) & 96.5 & 0.03 & 0.18 (0.006) & 92.3 & 0.16 & & 58.2 & 3.76\
& P(0.5) & 0.06 (0.002) & 95.6 & 0.05 & 0.06 (0.002) & 95.6 & 0.05 & 0.10 (0.003) & 94.4 & 0.09\
& P(0.25) & 0.07 (0.002) & 95.3 & 0.06 & 0.02 (0.001) & 97.5 & 0.02 & 0.29 (0.009) & 90.6 & 0.23\
& Huber & 0.04 (0.001) & 96.4 & 0.03 & 0.04 (0.001) & 96.3 & 0.04 & 0.14 (0.004) & 93.2 & 0.12\
\
& & & &\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} n & method & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$ & mse & PCD & $\delta_{0.5}$\
(r)[3-5]{} (lr)[6-8]{} (l)[9-11]{} 100 & LS & 1.45 (0.059) & 82.9 & 0.85 & 8.53 (0.784) & 72.4 & 2.01 & & 54.9 & 4.22\
& P(0.5) & 0.94 (0.034) & 85.6 & 0.61 & 1.29 (0.058) & 83.3 & 0.86 & 2.27 (0.132) & 78.9 & 1.24\
& P(0.25) & 1.46 (0.051) & 81.5 & 0.96 & 1.78 (0.071) & 78.2 & 1.30 & 7.88 (0.422) & 68.1 & 2.46\
& Huber & 0.89 (0.034) & 86.1 & 0.57 & 1.46 (0.067) & 81.7 & 0.99 & 3.28 (0.157) & 75.1 & 1.65\
200 & LS & 0.84 (0.035) & 86.6 & 0.53 & 3.85 (0.358) & 77.6 & 1.43 & & 55.9 & 4.09\
& P(0.5) & 0.44 (0.016) & 90.0 & 0.29 & 0.60 (0.024) & 89.0 & 0.39 & 0.87 (0.034) & 86.3 & 0.56\
& P(0.25) & 0.69 (0.025) & 87.0 & 0.49 & 0.75 (0.024) & 85.5 & 0.59 & 3.08 (0.179) & 75.3 & 1.58\
& Huber & 0.43 (0.016) & 90.3 & 0.28 & 0.66 (0.025) & 87.7 & 0.47 & 1.32 (0.050) & 82.4 & 0.87\
400 & LS & 0.44 (0.020) & 90.3 & 0.28 & 2.34 (0.393) & 82.4 & 0.95 & & 55.9 & 4.09\
& P(0.5) & 0.23 (0.009) & 92.9 & 0.16 & 0.28 (0.011) & 92.5 & 0.19 & 0.39 (0.015) & 90.8 & 0.26\
& P(0.25) & 0.33 (0.011) & 91.0 & 0.23 & 0.36 (0.012) & 90.1 & 0.27 & 1.25 (0.048) & 82.8 & 0.82\
& Huber & 0.22 (0.008) & 93.1 & 0.15 & 0.31 (0.012) & 91.7 & 0.21 & 0.60 (0.022) & 88.0 & 0.43\
800 & LS & 0.23 (0.009) & 93.0 & 0.15 & 0.90 (0.057) & 86.2 & 0.60 & & 56.3 & 4.03\
& P(0.5) & 0.11 (0.004) & 95.0 & 0.07 & 0.14 (0.005) & 94.8 & 0.09 & 0.18 (0.006) & 93.6 & 0.12\
& P(0.25) & 0.17 (0.006) & 93.7 & 0.12 & 0.18 (0.006) & 93.0 & 0.14 & 0.59 (0.017) & 87.3 & 0.44\
& Huber & 0.10 (0.004) & 95.1 & 0.07 & 0.15 (0.006) & 94.2 & 0.11 & 0.29 (0.010) & 91.4 & 0.21\

We firstly notice that lsA-learning works much worse under the heterogeneous errors, while all other methods are generally less affected by the heterogeneity of the errors. When the baseline function is misspecified as in Model I, under the homogeneous normal errors, RR(H) works slightly better than lsA-learning, while $\mathrm{RR}(\rho_{0.25})$ works the worst. The difference in general is small. For the homogeneous log-normal errors, again RR(H) works the best, while $\mathrm{RR}(\rho_{0.5})$ and $\mathrm{RR}(\rho_{0.25})$ work slightly worse. Here lsA-learning has the worst performance. Under the homogeneous Cauchy errors, the lsA-learning is no longer consistent and work the worst. Both $\mathrm{RR}(\rho_{0.5})$ and RR(H) have good performance under the homogeneous Cauchy errors. When baseline function is correctly specified as in Model II, under homogeneous normal errors, lsA-learning performs the best. However, in this case RR(H) also has a very close performance. Under homogeneous log-normal errors, $\mathrm{RR}(\rho_{0.25})$ work the best and lsA-learning work the worst. Under homogeneous Cauchy errors, $\mathrm{RR}(\rho_{0.5})$ has the best performance and RR(H) has a close performance. lsA-learning is again not consistent.

Examples with error terms interacted with treatment {#examples-with-error-terms-interacted-with-treatment .unnumbered}
---------------------------------------------------

We consider the following model with p=2, $$Y_i=1 + 0.5\sin[\pi(X_{i1}-X_{i2})]+
0.25(1+X_{i1}+2X_{i2})^2+(A_i-\pi({{\bm X}}_{i})){{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}_i+\sigma({{\bm X}}_{i},A_i)\epsilon_i,$$ where ${{\bm X}}_{i}=(X_{i1},X_{i2}){^{\mbox{\tiny {\sf T}}}}$, $\tilde{{{\bm X}}}_i=(1,{{\bm X}}_i{^{\mbox{\tiny {\sf T}}}}){^{\mbox{\tiny {\sf T}}}}$, $\sigma({{\bm X}}_{i},A_i)=1+A_i d_0 X_{i1}^2$, ${{\mbox{\boldmath $\theta$}}}_0{^{\mbox{\tiny {\sf T}}}}=(0.5,2,-1)$ and $X_{ik}$ are i.i.d. Uniform\[-1,1\]. We take linear forms for both the baseline and the contrast functions, where $\varphi({{\bm X}};{{\mbox{\boldmath $\gamma$}}})={{\mbox{\boldmath $\gamma$}}}{^{\mbox{\tiny {\sf T}}}}\tilde{{{\bm X}}}$, $C({{\bm X}};{{\mbox{\boldmath $\beta$}}})={{\mbox{\boldmath $\beta$}}}{^{\mbox{\tiny {\sf T}}}}{{\bm W}}$ and ${{\bm W}}=(\tilde{{{\bm X}}},X_{1}^2,X_{2}^2,X_{1}X_{2})$. $d_0=5$, 10 or 15. The error terms $\epsilon_i$ follows i.i.d. N(0,1) or Gamma(1,1)-1 distribution. The propensity scores $\pi(\cdot)$ are known, and we consider the non-constant case ($\pi({{\bm X}}_i)=\mathrm{logit}({{\bm X}}_{i1}-{{\bm X}}_{i2})$) here. The simulation results are given in Table \[table:interacted\_nonconstant\_ps\].

  ---------------------------------------------------------- ------- ----- ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------
                                                                                                                                                                                                                                                                                    
  (r)[4-6]{} (lr)[7-9]{} (lr)[10-12]{} (lr)[13-15]{} Error   $d_0$   n       $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$   $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$   $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$   $ \delta_{\mu}$   $ \delta_{0.5}$   $ \delta_{0.25}$
  Normal                                                     5       100                0.19              0.19               0.36              0.20              0.20               0.34              0.30              0.30               0.23              0.17              0.17               0.33
                                                                     200                0.11              0.11               0.28              0.13              0.13               0.25              0.21              0.21               0.12              0.11              0.11               0.23
                                                                     400                0.06              0.06               0.21              0.08              0.08               0.17              0.17              0.17               0.06              0.06              0.06               0.16
                                                                     800                0.03              0.03               0.16              0.06              0.06               0.11              0.15              0.15               0.04              0.04              0.04               0.11
                                                             10      100                0.29              0.29               0.93              0.24              0.24               0.88              0.44              0.44               0.50              0.24              0.24               0.88
                                                                     200                0.21              0.21               0.92              0.18              0.18               0.84              0.37              0.37               0.34              0.17              0.17               0.83
                                                                     400                0.13              0.13               0.87              0.14              0.14               0.75              0.32              0.32               0.25              0.12              0.12               0.75
                                                                     800                0.08              0.08               0.80              0.11              0.11               0.64              0.28              0.28               0.21              0.08              0.08               0.64
                                                             15      100                0.35              0.35               1.58              0.27              0.27               1.51              0.53              0.53               0.72              0.26              0.26               1.51
                                                                     200                0.29              0.29               1.56              0.21              0.21               1.47              0.50              0.50               0.54              0.20              0.20               1.47
                                                                     400                0.21              0.21               1.58              0.17              0.17               1.37              0.48              0.48               0.39              0.15              0.15               1.38
                                                                     800                0.14              0.14               1.52              0.14              0.14               1.26              0.45              0.45               0.31              0.12              0.12               1.27
  Gamma                                                      5       100                0.18              0.21               0.34              0.20              0.17               0.24              0.28              0.18               0.14              0.18              0.15               0.21
                                                                     200                0.10              0.14               0.29              0.13              0.10               0.15              0.21              0.11               0.07              0.11              0.07               0.13
                                                                     400                0.06              0.09               0.23              0.10              0.05               0.10              0.18              0.07               0.04              0.07              0.03               0.08
                                                                     800                0.03              0.06               0.19              0.08              0.03               0.06              0.16              0.06               0.03              0.06              0.02               0.07
                                                             10      100                0.27              0.34               0.90              0.28              0.25               0.67              0.46              0.21               0.33              0.28              0.22               0.62
                                                                     200                0.20              0.32               0.94              0.21              0.16               0.57              0.43              0.14               0.24              0.21              0.13               0.49
                                                                     400                0.13              0.27               0.92              0.16              0.09               0.46              0.38              0.10               0.18              0.15              0.06               0.39
                                                                     800                0.08              0.21               0.85              0.13              0.05               0.40              0.35              0.07               0.16              0.13              0.03               0.35
                                                             15      100                0.34              0.55               1.49              0.33              0.37               1.09              0.59              0.25               0.46              0.33              0.33               0.99
                                                                     200                0.27              0.54               1.57              0.26              0.29               1.00              0.60              0.19               0.31              0.27              0.23               0.85
                                                                     400                0.19              0.50               1.56              0.20              0.21               0.88              0.61              0.15               0.21              0.22              0.14               0.70
                                                                     800                0.12              0.47               1.58              0.17              0.14               0.76              0.62              0.15               0.18              0.19              0.09               0.63
  ---------------------------------------------------------- ------- ----- ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------ ----------------- ----------------- ------------------

  : Summary results with non-constant propensity scores when errors interacted with treatment. Least square stands for lsA-learning. Pinball(0.5) stands for robust regression with pinball loss and parameter $\tau=0.5$. Pinball(0.25) stands for robust regression with pinball loss and parameter $\tau=0.25$. Huber stands for robust regression with Huber loss, where parameter $\alpha$ is tuned automatically with R function rlm.[]{data-label="table:interacted_nonconstant_ps"}

Based on Theorem 6 of the main paper, $\delta_{\mu}$ column for the lsA-learning method in Table \[table:interacted\_nonconstant\_ps\] converges to 0 as sample size increases. Under Normal error terms, we have $\delta_{0.5}=\delta_{\mu}$. Thus, the $\delta_{0.5}$ column for the lsA-learning method under Normal error also converges to 0. All other columns in Table \[table:interacted\_nonconstant\_ps\] converge to a positive constant instead of 0 as sample size goes to infinity. $\mathrm{RR}(H)$ and $\mathrm{RR}(\rho_{0.5})$ perform similarly in Table \[table:interacted\_nonconstant\_ps\]. We also find even though lsA-learning outperform all other methods in $\delta_{\mu}$ when sample size is large. It may be worse than $\mathrm{RR}(\rho_{0.5})$ and $\mathrm{RR}(H)$ when sample size is small due to the fact that lsA-learning is inefficient under the heteroscedastic or skewed errors. Last, we find that lsA-learning, $\mathrm{RR}(\rho_{0.5})$ and $\mathrm{RR}(\rho_{0.25})$ perform best at the columns $\delta_{\mu}$, $\delta_{0.5}$ and $\delta_{0.25}$ accordingly. The reason is given in the Remark under Theorem 2 of the main paper.
