Smoothed Concordance-Assisted Learning for Optimal Treatment Decision in High Dimensional Data

In order to give precised treatment decisions, we have to take patients’ heterogeneity into account. Decision rules based on patients’ own features are called individualized treatment rules. In a binary treatment decisions setting, $\mathcal{A}=\{1,0\}$ are the treatment indicators and

are the i.i.d. observations, where $Y_{i}$ and $X_{i}$ are the outcome and features of subject $i$. Suppose $Y_{i}^{*}(A_{i})$ is the expected outcome for subject $i$ with features $X_{i}$ and treatment $A_{i}$. If $Y_{i}^{*}(1)>Y_{i}^{*}(0)$, treatment $1$ is more favorable to subject $i$ compared to treatment $0$. Assuming the linear decision rule $d(X)=\mathbbm{1}(\beta^{T}X>c)$, a treatment regime is determined by the coefficients $\beta$ and the threshold $c$. Mean treatment outcomes can be modeled as

An optimal treatment regime is the treatment regime which yields the most favorable mean treatment outcomes. Optimal treatment decision with multiple decision timepoints is called optimal dynamic treatment regimes. There are currently two leading dynamic treatment learning approaches. The first one is Q-learning by Watkins, (1989). Watkins and Dayan, (1992) proved the convergence property of Q-learning. Song et al., (2011) extended Q-learning into penalized Q-learning (PQ-learning). The second one is Advantage Learning (A-learning) by Murphy, (2003). Blatt et al., (2004) showed that A-learning is likely to have smaller bias than Q-learning.

Due to the rapid accumulation of patients’ heterogeneity data, people started to consider the high dimensional setting. Zhu et al., (2019) extended Q-learning to high dimensional Q-learning, with a focus on the two-stage dynamic treatment regimes. Shi et al., (2018) proposed Penalized A-learning (PAL), which applied the Dantzig selector under the A-learning framework, with a focus on two or three-stage settings.

For a single timepoint treatment decision, Zhang et al., (2012) used a doubly robust augmented inverse probability weighted estimator to handle the possible model misspecification issue. Zhao et al., (2012) adopted the support vector machine framework and proposed outcome weighted learning (OWL). Song et al., (2015) modified OWL to be penalized outcome weighted learning (POWL). Lu et al., (2013) introduced a robust estimator which doesn’t require estimating the baseline function and is easy to perform variable selection in large dimensions. CAL by Fan et al., (2017) used the contrast between two treatments to construct a concordance function, and then used the concordance function to estimate the coefficients $\beta^{*}$. CAL has a fast convergence rate, and it does not assume that the contrast function is a linear combination of the features. They assumed the stable unit treatment value assumption (SUTVA) in Rubin, (1980),

and the no-unmeasured-confounders condition,

The proposed concordance function $C(\beta)$ is

where the contrast function $D(X_{i})=E[Y_{i}|A_{i}=1,X_{i}]-E[Y_{i}|A_{i}=0,X_{i}]$. Within the concordance function, an important assumption is that $\beta^{*T}X_{i}$ is concord with $D(X_{i})$, which means

where $Q(x)$ is an unknown monotone increasing function. Using the propensity score $\pi(X_{i})=P(A_{i}=1|X_{i})$, their $D(X_{i})$ is estimated by its unbiased estimator

where $v(X_{i,}\hat{\theta})$ can be any arbitrary function independent of $A_{i}$, while it is usually chosen as the mean response of those patients who received treatment $0$. The true $\beta^{*}$ and the true threshold $c_{0}$ in the linear decision rule are estimated by

and

In the high dimensional setting, however, due to the indicator function inside the concordance function, the optimization of CAL is somewhat difficult. In order to deal with this optimization issue of CAL, Liang et al., (2018) proposed SCAL, which used hinge loss to replace $\mathbbm{1}(x)$, and they added the $L_{1}$ penalty for variable selection purpose. They estimate $\beta^{*}$ and $c_{0}$ by

and

It achieved $n^{1/2}$-rate but induced a relatively large bias due to the difference between the hinge loss and the original 0-1 loss.

Our way to overcome the optimization difficulties in CAL under high dimensional data is based on the smoothing method proposed by Han et al., (2017). A family of sigmoid functions are used to substitute the indicator function in CAL. The $L_{1}$ penalty is also added. We then employ the coordinate descent algorithm for optimization. Our method is called Smoothed Concordance-assisted Learning (SMCAL).

Compared to SCAL, our SMCAL has a slower convergence rate but achieves a much smaller bias. Numerical comparisons with POWL, SCAL and PAL demonstrated the advantage of our method, especially when the sample size is relatively large. Other than the continuous covariates cases in the SCAL settings, we expand our SMCAL to the discrete covariates cases.

There are three main differences between our work and Han et al., (2017). Firstly, they focused on the Maximum Rank Correlation (MRC) proposed in Han, (1987), but our smoothing is applied on the concordance function. Secondly, their paper studied Monotone Transformation model, but we focus on Monotone Single Index Model, and we have a faster convergence rate. Thirdly, we assume weaker assumptions instead of their normality assumption posted on the covariates, and we expand our method to discrete cases.

The rest of this paper is organized as follows: Section 2 introduces SMCAL and describes the coordinate descent algorithm. Section 3 displays the $L_{2}$-error rate in continuous cases and the $L_{1}$-error rate in discrete cases, which are the main results of our paper. In Section 4, we conduct numerical comparisons with SCAL, PAL and POWL. A real data application to STAR*D dataset is provided in Section 5. The proofs of lemmas and theorems in Section 3 are left to the supplementary material.

In this section, we establish the error rates for $\hat{\beta}_{r(\gamma),\alpha_{n}}-\beta^{*}$ under continuous and then discrete covariates cases. Let $\lambda_{n}\geq 2||\nabla\hat{L}_{n}(\beta_{r,\alpha_{n}}^{*})||_{\infty}$. Assume the sparsity of the real coefficients $\beta^{*}$. Let $S$ represent the nonzero indexes in $\beta^{*}$, $s=|S|$, and define $S^{c}=\{1,2,...d\}\setminus S$. Let $X_{i,I}$ represent all the entries of $X_{i}$ whose indexes belong to set $I$. The notation $\lesssim$ means that if $a_{n}\lesssim b_{n}$, or equivalently $b_{n}\gtrsim a_{n}$, then there exist $C>0$ and $N_{0}>0$, such that for all $n>N_{0}$, we have $|a_{n}|\leq C|b_{n}|$. If $a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$, then we write $a_{n}\asymp b_{n}$.

Four methods are compared in this section: POWL in Song et al., (2015), PAL in Shi et al., (2018), SCAL in Liang et al., (2018) and SMCAL. POWL method minimizes

where $g(X_{i})=\beta^{T}\cdot X_{i}+\beta_{0}$ is the decision function and the tuning parameter $\lambda$ is selected by the largest IPW estimator. PAL is a dynamic treatment regime approach but can be used in the single-stage problem because PAL estimates the latest stage at first. PAL perform variable selection using penalized A-learning, then use unpenalized A-learning to solve the coefficients. We used the R package provided by the authors of PAL to get its results.

To evaluate the estimated coefficients, we report the MSE after normalization. To evaluate the variable selection accuracy, we report the Incorrect Zeros (true coefficient is zero but the estimation is nonzero) and Correct Zeros (true coefficient is zero and the estimation is zero). To evaluate the estimated treatment regime, we report the Percentage of Correct Decision (PCD) and the mean response if we follow the estimated treatment regime (Estimated Value). Estimated Value is calculated by drawing 1000 new samples. Standard errors of the PCDs and Estimated Values are in the parenthesis.