Learning Individualized Treatment Rules with Many Treatments: A Supervised Clustering Approach Using Adaptive Fusion
Authors: Haixu Ma
Donglin Zeng
Yufeng Liu
Venue: NeurIPS 2022 Accept (2022)
Keywords: Fusion Penalty
High-dimensional Regression
Individualized Treatment Rule
Precision Medicine

Abstract: Learning an optimal Individualized Treatment Rule (ITR) is a very important problem in precision medicine. This paper is concerned with the challenge when the number of treatment arms is large, and some groups of treatments in the large treatment space may work similarly for the patients. Motivated by the recent development of supervised clustering, we propose a novel adaptive fusion based method to cluster the treatments with similar treatment effects together and estimate the optimal ITR simultaneously through a single convex optimization. The problem is formulated as balancing \textit{loss}$+$\textit{penalty} terms with a tuning parameter, which allows the entire solution path of the treatment clustering process to be clearly visualized hierarchically. For computation, we propose an efficient  algorithm based on  accelerated proximal gradient and further conduct a novel group-lasso based algorithm for variable selection to boost the performance. Moreover, we demonstrate the theoretical guarantee of recovering the underlying true clustering structure of the treatments for our method. Finally, we demonstrate the superior performance of our method via both simulations and a real data application on cancer treatment, which may assist the decision making process for doctors.

Learning Individualized Treatment Rules with Many Treatments: A Supervised Clustering Approach Using Adaptive Fusion
Haixu Ma
Donglin Zeng
Department of Statisties and Operations Research University of North Carolina at Chapel Hill Chapel Hill, NC 27516 haixuma@1 live. unc. edu
Department of Biostatistics University of North Carolina at Chapel Hill
Chapel Hill, NC 27516 dzeng@email, unc. edu
Yufeng Liu
Department of Statistics and Operations Research
Department of Genetics
Department of Biostatistics University of North Carolina Chapel Hill
Chapel Hill, NC 27516 yfliu@email. unc. edu
Abstract
Learning optimal Individualized Treatment Rule (ITR) is very important prob- lem in precision medicine. This paper is concerned with the challenge when the number of treatment arms is large, and some groups of treatments in the large treat- ment space may work similarly for the patients. Motivated by the recent develop- ment of supervised clustering, we propose novel adaptive fusion based method to cluster the treatments with similar treatment effects together and estimate the optimal ITR simultaneously through single convex optimization. The problem is formulated as balancing loss+ penalty terms with a tuning parameter, which allows the entire solution path of the treatment clustering process to be clearly visualized hierarchically. For computation, we propose an efficient algorithm based on accel- erated proximal gradient and further conduct novel group-lasso based algorithm for variable selection to boost the performance. Moreover, we demonstrate the theoretical guarantee of recovering the underlying true clustering structure of the treatments for our method. Finally, we demonstrate the superior performance of our method via both simulations and real data application on cancer treatment, which may assist the decision making process for doctors.
Introduction
Data-driven individualized decision making problems have received lot of attentions in precision medicine [MI. 2. C. H]. Due to possible significant heterogeneity of treatment effects among individ- uals, it is necessary for decision makers precisely tailor the treatment decision rules to different subgroups of individuals rather than simply implement the traditional "one size fits all" medical procedures. The main goal of the paper is to recommend the optimal Individualized Treatment Rule (ITR), mapping from the covariate space based on patients' characteristics to the treatment assign- ment, that optimizes the expected value of specified reward, known the value function [5].
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

There are many existing machine learning based approaches for estimating the optimal ITR. Some methods estimate the treatment effect related functions under some prespecified outcome mean mod- els. The estimated ITR is induced by maximizing the estimated treatment effects conditional on patients' covariates. These methods include Q-learning [5]], A-learning [B. Z] weighted ordinary least square based method [8] and D-learning [B]. Some other methods circumvent modeling the treatment effect functions by directly estimating the ITR that maximizes the Inverse Probability Weighting (IPW) version of the value function [TTA I. T2, IK [El]. To overcome the potential mis- specification problems for these methods, [US.. 16, 17, [8] proposed the doubly robust IPW methods.
Problem description. This paper is concerned with the following two questions: (1) Given large number of treatment options but limited observations of some specific treatments, how can we effectively estimate the optimal ITR? (2) Some treatments in the large treatment space may work similarly for patients. How can we identify this homogeneous treatment structure and cluster the treatments with similar treatment effects together to reduce the dimension of the treatment space?
For question (I), the estimated ITR using either the model-based methods or the IPW IPW-based methods can become inaccurate with large variability and numerical instability due to the insufficient number of observations for certain treatments. For question (2), it is common that the drugs are developed based on intervening the same disease symptoms and mechanisms but may be produced by different pharmaceutical companies. [[[G] and [[9] directly combined the treatments with similar treatment effects into certain classes of treatments based some prior knowledge. In general, it is desirable to propose data-driven methods to identify the homogeneous treatment structure automatically. How- ever, to our best of knowledge, few existing methods deal with clustering treatments because most literature only consider binary or a moderate number of treatments.
Related work and our method. Our goal is to cluster the heterogeneous treatment effects and estimate the ITR under the Q-learning framework [20] 5]] Some methods have been developed to address the problem of identifying subgroups with homogeneous relationships. In supervised
learning, with the special focus on exploring the homogeneous relations between the covariates and
response, [2I]]. [22] and [23] estimated the group structure for the regression coefficients of covari- ates with various fusion penalties. To identify subgroups from heterogeneous population, [[4]] and [25]] assumed the population comes from a mixture of subgroups with their own distributions and utilized the mixture model analysis to classify the observations. However, the mixture model may cause misspecification problems in practice since it needs to specify the underlying distribution for the data and the number of mixture components in the population. Moreover, these supervised learning methods only utilized the information of response and covariates without considering the treatment assignment. Hence, they are not applicable to the ITR problem. In unsupervised learning, clustering analysis is a popular tool. It is usually used to detect the similarities of observations using predefined distance measure. In partieular, popular visualization method is to draw dendrogram of the hierarchical clustering using a "bottom up" approach [26]..
Regarding to our problem of clustering the treatment effects, we would like to adapt the idea of latent supervised clustering [2Z 28). We are interested clustering the relationship characterized by the regression of "response covariates treatments" where "x" refers to the interactions be- tween covariates and treatments. To achieve this goal, we model the treatment effeets as the regres- sion problem with rreatment-speciff coefficient vectors. We implement an adaptive pairwise fusion penalty of the l fi-distance between the treattent-specific coefficient vectors in order to merge these coefficient vectors with similar values together. This is equivalent to clustering the treatments with similar treatment effects. We formulate the clustering process as a convex minimization problem in- volving loss + fusion penalty, with a tuning parameter balancing these two terms. With this convex formulation, we achieve to maximize the goodness of fit for estimating the heterogeneous treatment effects, while at the same time clustering the treatments without the need of prior knowledge.
Contribution. The main contributions of our paper can be summarized as follows. (1) We propose the Supervised Clustering approach using the Adaptive Fusion (SCAF) to identify possible homo- geneous group structure within the large treatment space and estimate the ITR in a more effective way. (2) Compared with the two-step method [29], we simultaneously cluster the treatments with similar treatment effects and estimate the optimal ITR. Our method combines the unsupervised and supervised learning together within a single convex optimization problem. (3) We do not need to specify the number of treatment clusters and the entire solution path of the clustering process can

---

be visualized with new dendrogram by using a "bottom up" approach. (4) Although [27] and [28]] also adopted the supervised clustering technique, they focused on identifying the heterogeneous subpopulations and aimed to cluster the subject-specific regression coefficients. In contrast, SCAF only needs to estimate the treatment-speciffc coefficients and the dimension of the estimated param- eters is dramatically reduced. (5) For computation, we propose an effective algorithm to solve the convex minimization problem based on the accelerated proximal gradient algorithm [[O]] Specif- ically, in order to boost the performance of our fusion algorithm, we propose novel group-lasso based algorithm to classify the covariates into homogeneous variables and heterogeneous variables respectively. Thus, the fusion penalty only needs to be imposed to the heterogeneous variables that have interaction with treatments. (6) For the theoretical study, we prove the consistency of the re- gression coefficients and demonstrate that SCAF is able recover the true underlying clustering structure with probability 1. (7) We conduct both simulation studies and a real data analysis on cancer treatment to illustrate the superior performance of SCAF.
2
Methodology
We first introduce the ITR problem under the regression bassed framework. Consider the training data (zi, as 344 for i = ...... n as i.i.d. realizations from the joint distribution of ( (Z, A, Y), where 2, 2 C. Iº' denotes the patient's prognostic variables, a, A = (1.2...., ,M} is the treatment assignment, and y/ IR is the observed reward for each patient i. Let (Y(a)) aE E H be the potential outcome. In addition, define the propensity scores of treatment p(a[z) := P(A a|Z z) for a E A. An ITR D E D is a map from the covariate space Z to the treatment space A. Here, D C. A is prespecified ITR class. The value function of an ITR is defined V(D) B[Y ]- Assuming that larger reward is better, our goal is to find the optimal ITR D* e 'D that maximizes the value function, i.e., D* E arg maXDeD V(D).
We assume the following identifiability assumptions are satisfied [ (1) consistency: Y Eag = alY(a); (2) no unmeasured confounders: for each a E A,Y(() 1 A Z: (3) positivity: p(alz)
0 for any z Z. Note that the no unmeasured confounders as- sumption can be relaxed. If there are possible confounding issues as in the observational stud- ies, proposed method can be further generalized deal with such issues by using propen- sity scores. Based on the three identifiability assumptions, the value function can be written as V(D) JE [ [Z) from D*) E arg 1 Y Z A a] for each Z. This motivates us to estimate the conditional treatment effect E Y Z. A]. Specifically, we consider the following regression model:
a]IE[Y| A al]. Therefore, the optimal ITR D can be derived
Y Mo(Z) + E I[A a a)T(Z; So) E, s. >I T(Z: G.) = 0; E[c Z, A] = 0; varl var[c[2.A] < +00,
(I)
where the redundant function Mo(Z) is the main effect of treatments, and T(Z: Ca) is the inter- action effect between treatment and the covariates. We assume that the interaction effect for each treatment a A can be characterized by the treatment-specific parameter Co E R2 that has the same dimension as the covariates Z. Here, we consider that the covariates 2 contain the intercept term so that the treatment-specinee effect is included in T(Z: Ca). Note that a sum-to- zero constraint for the interaction terms is assumed for identifiability of the regression function. Given an estimated parameter Ç = (ŚI .., "T, T, the estimated ITR can be equivalently written as D(z) E arg max.s.4 T(z; E.) for each = E 2.
With Model (I) in place, we propose our SCAF method to cluster the treatments with similar treat- ment effects into a combined treatment group. To achieve this, we aim to estimate and fuse the Ça's 's into several groups. The group structure of the Ça's can equivalently represent the cluster- ing results for the treatments. Note that the main effect function My (Z) can be estimated from the weighted parametric the main effect can be found in Appendix A. of supplementary materials. Denote the estimation
nonparametric regression models [62]] More details about estimating
of Mu as My. and let Y be the residual Y - Mo ZZ) In this article, we focus on linear interac- tion effects, i.e., T(z, Sa) ITaa. Let Ç = (CT .C....., (CT Denote the true value of Ç to be - = (SP" ,..., ( ) R.M, where ch E RR is the true parameter for treatment a I,...; M. In order to estimate and cluster Ça 's, we consider the following optimization problem by imposing
3

---

pairwise fusion penalty to each pair of the seratment-speciffe parameters:
min [ E.
E. A a]ZT EiSICEEM PA. (lls: Sell2)
(2)
where .)..)) is the prespecified loss function to characterize the goodness of fit, I-ll1 is the & norm of a vector, PA,. is the penalty function to encourage fusing Sa's into groups, and A, is the tuning parameter. Here, E, denotes the empirical mean of the training data. Specifically, the optimization problem (2) can be interpreted by maximizing the goodness of fit, while at the same time minimiz- ing the heterogeneity among treatments. When the true treatment effects C. s have homogeneous clustering structure, the estimated parameter Ça s are expected to recover this structure.
In practice, one may assume that only certain elements of Z Rp (denoted as X R-) contribute to the interaction term T(Z; Sa) while others (denoted as V E R) only show up in the main effects Mo (Z) [28]] Without loss of generality, assume the first d components of Z are the hetero- geneous variables and the remaining p variables are the homogeneous variables v.c. Z = VT)T. Let the coefficients for treatment a be C. = ((Ga.I...... Sa.a), (Sa.d+1 .... Sap))) = (8)-))) where Sa. the regression coefficient of covariate Z& =....., .......)) p). and B, and % are the coefficients for X and V respectively. By removing the main effect Mo (Z) and consider the residual Y. the true regression coeffecient of V for each treatment should be 0, i.e., 7 = 0 for all a E A. In this case, the penalty term should only be imposed on the coefficients of X, i.e. 3. If prior knowledge of whether the covariate belongs to X or V is unknown, we propose the group-lasso based algorithm in Section B to classify Z into homogeneous variables V and heteroge- neous variables X. Denote the training data (æ;, Vi, 0i, Jk) for i I,...,74, ni, where æ; E IR and We rewrite Model (II) by decomposing T [2 Ca) into two parts:
v; E RR-dd
Y = Mo(X, V) + L..A I[A = a) + T(V:7)) +
= MX V) + E.EA I[A = a]XTBa+ + c,
s.t. E.EA XTB, 0; E[e[X, V. V] 0; var V, A] < <00.
After distinguishing X and V from the whole variables Z, the optimization problem in (L) becomes
((Am)) arg min B-(3] (T)
(o.(0:A.) (.: B. [C(, Σ =|= a]xT
Σ 1slciEM
ps., (13 - -
(4)
where (B; An.) is the optimization function, and (An)) is the estimated coefficients based on the training data. Note that if we do not detect the homogeneous variables X from Z, then the pairwise penalty in (2) would be applied to the p-dimensional vector C. 's. For comparison, in (H), only the subvector, i.e., the d-dimensional vectors Ba 's are included in the penalty term. Our empirical results show that our fusing algorithm becomes more effective to solve the optimization problem since the low dimensional vector easier to be merged together.
The minimization problem (H)) can be viewed as the supervised clustering process. In the regression- based Q-learning framework, supervised clustering can be interpreted by clustering the treatment effects (relationships described by the regression model: reward Y ~ covariates x treatments A) with the penalized regression problem in (H)), which is formulated by balaneing loss + penalty terms with the tuning parameter A,. The pairwise fusion penalty encourages B, to merge together so that the treatments with similar treatment effects can be clustered into the same group.
3 Algorithms for SCAF
Group-lasso based algorithm to classify Z into X and V. After obtaining the residual Y by subtracting M2( from Y. we propose the following group-lasso based algorithm to identify X and V. Recall that (CT.... STy) and C.. (Sa,1, TT for treatment
A. Note that
Sa,k is the regression coefficient of covariate Ze (k 1,.... p) for treatment a E A. We rewrite Ç by sorting it based the covariate order, i.e., E (ET ..... TT)T, where &, (Gi,k..... (M,,,)) for 1, P. is the coefficient vector of Z₁ for all treatments. By the definition of the homogeneous variables V E R the true regression coefficient Ca, 0 for all k d + 1,..., p and all a A. Hence, the whole true coefficient vector can be written by the following sparse group structure,

---

ie, 5-' = (Ef" ......2 "'' .'. where the coefficients of each covariate are considered as group. This motivates us to implement group-lasso penalty to shrink the whole group vector & that corresponds to covariate 2 to 0 to achieve group sparsity. If & = O, then Z should be classified as V. otherwise it should belong to X. In particular, consider the following optimization problem:
f = arg min En E. IA- = a](Y - Ek=1 Sa,kZa) + Aglasso E I1Eallo
(5)
where 1-l1, is given base norm for each &s and Aglasoo is the tuning parameter. Typical choices of base norms include the ||:| or ll-llcc norm, so that the group-lasso penalty term becomes the block f//l2 E1/La norm, i.e., the norm on all groups with the lg or lx norm for each group.
Adattive proximal gradient algorithm. Next we introduce our main algorithm to solve problem (=)) We consider the Ey-loss function for C(..)) and choose the adaptive pairwise fusion penalty weighted by w for (.). Let B = [B]. B BM) Specifically, the optimization problem (i) can be written in the following compact form
B(An)) arg ming (Qn(3:Am) == (.)
(6)
where
(B) = JEm - - = x[XT3,)3] Pa( = An EiSl<IEM Wt, 1B Bella
(7)
Here, W1,1 is the weight of penalizing the di-dittance between B1 and 3.. We iteratively solve the problem coming from (65) and (0)) In practice, we find that W1, = min (Bw, 1/1l 3elli3 can be good option where G (ØT ...... BT,) T is estimate for ,3 from the last iteration. The initial estimate for 3 can be obtained from (5) in the group-lasso step. Moreover, B. is a prespecified constant to upper bound W1. to handle small distance of the pair (Bi, Be). The adaptive weight helps to adjust for the potential bias created by the f1-penalty. Intuitively, large weight expected to be assigned to treatments within the same group, otherwise, small weights are implemented.
Note that both T. (3) and (33) are convex functions for B, and (3) has Lipschitz gradient. Thus, we can utilize the accelerated proximal gradient algorithm [GO]] to iteratively solve the prob- lem coming from (E) and (α.). Denote the gradient vector of (5) as (L) and the Lipschitz constant of L, (3) as ln- The computational cost of our algorithm mainly includes two follow- ing parts: (1) calculate in and evaluate the gradient (L) for the updated B; (2) compute the proximal operator of P.. which is defined as prox P. (0) := arg ming (3) 2s. 3 (1 for any updated (9 and the step size 8n O. For fixed An. the main steps of the proposed algo- rithm for SCAF are summarized in Algorithm 1. The tuning parameter Ay can be tuned by cross validation based on the IPWE of the value function performed on the validation data. For z Z, the final estimated ITR (22) obtained from randomly sampling one treatment from the optimal estimated treatment group with probability proportional to the propensity score. The tuning crite- ria is E, [((2)))))
p(Alz) YJ/E,, [I which is larger the better [33 [0]. Here, E, denotes the p(A[Z)
empirical mean for the validation data.
4 Theoretical properties
In this section, we establish the theoretical guarantee for consistency of the estimated regression coefficients. For Model (L)) consider the M treatments can be partitioned into K disjoint treatment groups { where the effects of treatments within the same treatment group are identical. Here, the group number K is typically unknown in practice. Let 3 = (3) ,39 ...., 330)"" treatment-speciffic regression coefficients in (L). The values of (00s from the same treatment group 9% for l 1,..., K are the same. In particular, we aim to develop the convergence of ((\,)) g as n 00, which equivalently demonstrates consistency of recovering the true group structure.
be the true
Let X E I nx be the design matrix of X. Denote X,, E Re.XP xP the submatrix of X where the ob-
servations X, are assigned to treatment a. Here, n, is the numher of patients receiving treatment and hence ΣM group structure, let a (a. aa2 .....
n. Further denote U diag( X, X..... XM) E Rp, Under the true
k T E IR Ky be the distinct values of 3 where af R
5

---

Algorithm 1: SCAF
Step I: Sort the observations based on the assigned treatment order.
Step 2: Remove the main effect MZ(Z) and calculate residual y.
Step 3: Implement group-lasso to identify heterogencous variables X from Z.
Step 4: Use adaptive fast proximal gradient algorithm to solve problem (F:
(1) Obtain the initial point ( from Step and set the desired tolerance co >
(2) Compute the Lipschitz constant in = A (U U) and set the step-size 8n = 1/h,, to = 1;
Let ((0)= == 3 and set wi, == min Bw, 1/1
3( [13 for p. (3) (1, A);
(4) For i 0, l...., keax. do:
a. Compute g((+1) & prox P() (3 8n VLn (3(1)) [28];
b. Update f++1 (1 +
482 1
C. Perform FISTA [6_]] with 3
g((11)
+ 4=1 (g(111) - =3(1);;
d. If (13(++11)
3 & co. then end the loop:
c. Update w min Bw,, 1/13((+1))
3 3(1+1) |1} for (()+)(()) (1. A);
(5) End of the main loop.
Step 5: Obtain the estimated ITR ((x) arg a x 33 for x X.
the true treatment group-specific coefficients for k 1...., K. Denote N₂. to be the number of patients whose assigned treatments belong treatment group GR for k I,..., K. Then we have 2k-1 Nu = n. In addition, use [Ga] denote the number of treatments in the k-th group where Ek=1 [9kl = .. For simplicity of notations, we assume that the treatment index has been well sorted based on the true group structure. Then, under the true group memberships G1, $1..... Gk. let H2 e RMxx xp be the submatrix of X where the observations in H, receive treatments in group Gn for k I,..., K. Further denote H diag(H, .2....,. x) RRnxKy In particular, under the true group structure, define the oracle estimator for a to be 2 :- arg min, j H Equivalently, a (HTH)-
The corresponding oracle estimator for B, denoted for, is obtained by expanding aer based on the true group structure.
Hereafter, we allow the number of treatments M, the true treatment group number K, and the covariates' dimension p to grow as the sample size n increases. We use Mn, Kn, and Pn de- note them respectively. In addition, for any given matrix G (G denote |G|
maxigies 21 G.yyl For vector g (g1...., T R*, denote || x maxisisa [gi]. We use 4, >> k, to denote k, o(1).
We first establish the concentration bound for consistency of the oracle estimator when the obser- vation numbers in each treatment group employ suitable distribution pattern. Specifically, let Nmi := mink=11.... N₂ be the minimum number of observations receiving the treatments from the same treatment group. We make the following regularity assumptions in the high-dimensional statistical literature [27. 35];
Assumption 1. Assume |X,12 for j 1.... ,r. and |UTU| i0 & BiMapn where B₁ > 0.
Assumption 2. The minimum eigen value Amin! (HTH) C,Niiin where C₁ 0.
Assumption 3. Let E = (EI..... e R" be the error vector in model (G). Suppose e has sub- Gaussian tail bound, i.c., there exists constant e1 +00) such that, for any vector E RR" and any > O, we have Pr([aTe] > ||a]] 2 t) < 2exp(-c,t").
Theorem 1 (Consistency of oracle estimator). Suppose Assumptions 1-3 hold. If Kn,n = o(n), and Nmin >> (Kapu)) log n, then with probability at least 2K,pa/n, we have |300 30 =
|@ - a & Φn, where on C1 -1/2 VKrPhNN Vn logn.
Since / (KaPh)n log & Nmin & n/ n/Kn, we have kny(Knn)) = o( n// logn). Hence, Kn must satisfy K,, = (log3/((ogn))1))) treatment groups. If we set Nuin = ón/Ks Kn where E (0, 1), then the bound on = 11/2 10-8-1 / n/nn/n.
It indicates that we cannot have too many

---

Our next step is to show that the oracle estimator 3 is local minimizer of the object function (n: in (6) with probability by giving a lower bound of the minimum signal difference between treatment groups. Then the convergence of ((,,.) - 3o can be established. To achieve that, we need assume the penalty function PA... (-) in (3) has sharp derivative around Ó to adjust the bias of the fusion penalty. Some popular penalties considered in [GG]] [67] and [G8] satisfy similar properties as the following assumption.
Assumption 4. Suppose px.. (-) is symmetric around 0 for each X.. and there exists a constant a > 0 such that px.. (1) is a constant for all & gdm. Assume Px.. (0) 0, and is differentiable around 0 with its derivative satisfying PAA.. (t) > OVv'logn) infi [Gu] for all < 2
Let b, inf ,ee A 133 B3112 infe* | ag I be the minimal difference of the regression coefficients between two groups. Intuitively, we need the distances between two different groups to be large enough so that we can recover the group membership. In particular, we make the following assumption. A similar assumption was also used in [27]] and [39].
Assumption 5. Suppose bn/Pn. adn and Xn >> On.
Theorem 2. Suppose Assumptions 1-5 hold, then there exists a local minimizer (((a) of the objec- rive function Qn Xn) in (B) such that Pr Pr((()n) = 3rr)
Theorem 1I together with Theorem D establish the convergence property of 3(Am) - 3o. For our implementation in Section EL we iteratively solve the minimization problem (8I with the adaptive fusion penalty shown in (Z)). We use 300 = (3(a) ,.......,,,,,, (()) I to denote the estimation of /3 in the i-th iteration. For 1, E A, the weight the i-th iteration (1)) specified as min( 1/11((
3 3(t
|1.. Hence, the derivative of pA. (.) evaluated at |/(() - -3(")), [|, which equals to X,
will become extremely large when | / 1 is close to 0, and it will be relatively small if the pair (8("). ((") are not fused together. Consequently, if 3 is close to 8 then the derivative of PA., (-) will approximately satisfy the properties in Assumption E within a neighbourhood around g Therefore, based Theorem A. the oracle estimator 30r can still be the local minimizer of
@n Xn) within a neighbourhood around 30 when our adaptive fusion penalty is utilized.
Experiments
We evaluate the finite-sample performance of our method using simulation and a real data applica- tion to the Patient-Derived Xenograft (PDX) study. We compare our method with the sfate-of-art methods for the ITR problem: (a) RF (Random Forest, [+]]; (b) PLS (l, Penalizzze Least Square method, [S)) (c) AD (Multi-armed Angle-based Direct learning method, [.)). However, none of these existing methods consider the treatment group structure in large treatment space.
Data generation for simulation. We generate 10-dimensional independent feature variables Zi..... Z1o. following U[-1, The outcome Y is normally distributed with E[Y|Z, A] = 1+ 22, + Z2 + 0.52, + To (Z,(Z)) A) and standard deviation 1, where To A) reflects the inter- action between the treatment and the feature variables. Specifically, we conduct the following three scenarios where the treatment effects have the homogeneous grouping structure G:
Scenario M = 10, To(Z, A) L. Zi 2a) (-I) 9.) 1):
2.g = (9, = (11,2,3,4,55 (6. 7, 7,8, 1033. = (1, 71, 72), and
Scenario 2. .M = 15, K = 3,9 = (91, = (11,2,33,4,5,, 10), (11 12, 13, 14, 15)). X = (1, 21, 22). and (2(,,4) = (2)(00--+++)))
+2, 22) e Ga] + (-0.2+
[Z)) E Gs]):
Scenario 3. M = 20, - (1,.,66666_|111222,,,666778,,8
(13, 14,
15, 16), (17, 18, 19, 20 (1 Zo, 7,), and = + +1.5% + 27g +
1.5 +1 27% 7 37g- A E (-1+-1 27, + 7a -ZJTTA
Ga] (-2 7 7 -Zs + Z 94] +1.5 + 1.5% + +552.)] 9a1).
As the number of treatments M increases from 10, 15 to 20, the above three scenarios correspond to the treatments having 2,3 and 5 treatment groups respectively. For each scenario, the true decision
7

---

boundaries are linear under the group domain. Take Scenario for example. For any x R'o, the optimal ITR (*) would recommend any treatment in G1 when 0.2 2 22 < 0, and any treatment in G2 when 0.2 - 21 22 3 0. The treatments in all the scenarios are assigned in the unbalanced structure, which means the propensity score of specific treatments can be small. The tuning parameter X, is selected by the 5-fold voossvvaliatiion. For each scenario, the train- ing sample sizes vary from 200, 400 to 600 and we replicate the simulations for 200 times. On an independently generated testing data of size 10, 000, we evaluate the above methods using (a) empirical value function; (b) group-based misclassification rate between the estimated group deci- sion rules and the optimal group decision rules. The empirical value function is calculated using E., p(a]2) A]] [ [53, I0]] where E., denotes the empirical mean for the testing data. Since the treatment effects are identical within the same treatment group, the misclassification rate
p(A]z)
under the group domain is equivalent to that under the individual treatment domain.
Results for simulation. The boxplots of the empirical value and the misclassification rate eval- uated on the independent testing data are shown in Figure and Figure C.I in the supplementary materials. In most cases, compared with RF, AD and PLS, our proposed SCAF method has supe- rior performance with higher empirical values, smaller misclassification rates, and extremely lower variabilities for both evaluation criteria. When the treatment space becomes large, RF, AD and PLS suffer from inaccurate estimation of the individual treatment effect due to the small amount of obser- vations for some specific treatments. Moreover, these methods do not consider the possible group structure in the treatment space and hence employ relatively large variance. As comparison, SCAF is able to recover the treatment group structure and reduce the dimension of treatment space by clus- tering the treatments with similar treatment effects into the same treatment group. We also show the ratio of exactly recovering the true treatment group structure of SCAF among the 200 replications in Figure C.2 in the supplementary materials. Furthermore, we use Scenario to demonstrate the advantage of our group-lasso step for finding the heterogeneous variables X from Z. By implement- ing the group-lasso before the fusion step, as shown in Figure C.3 in the supplementary materials, SCAF with the group-lasso step has better performance than that without the group-lasso step since we only need to fuse the lower dimensional subvectors Ba 's rather than Ca's.
Scenario 1
Scenario 2
Scenario 3
0.80
1.3
Log(Value)
(og(Value)
1.2
0.75
2.3 Segivaa
1.1
1.0
0.70
2.2
0.9
200
400
600
200
400
600
200
400
600
Sample.size
Sample.size
Sample.size
Methods
RF
PLS
SCAF
Figure 1: Boxplots of empirical value based on the testing data in simulations. Red dashed lines demonstrate the oracle values.
Scenario
Scenario 2
Scenario
LogiA11000)
100/A110000
LogiA-100)
Treatment
Treatment
Treatment
Figure 2: Solution path of the estimated treatment group structure as À increases. The true treatment group memberships are demonstrated with different colors. The red dotted horizontal lines show the best tuned À using cross-vale coossvalidation.. The estimated groups based on the best À are framed with grey dotted rectangles.
8

---

Figure 2 demonstrates the solution path of the estimated treatment group structure as the tuning parameter À increases based on one of the replications for the three scenarios. When À 0. SCAF is equivalent to RF and PLS since the fusion penalty is not imposed and does not show any clustering pattern. As À increases, the fusion penalty encourages the treatments to merge together and cluster them with the expected group structure. When À becomes large enough, all of the treatments will be merged together and only one treatment group will be formed. With the best tuned X, 2, and 5 groups are correctly selected for the three scenarios respectively.
Application Patient- Deriived- Xenograft (PDX) study. Due to the complexity of human cancer, there's 's significant heterogeneity of treatment effects among the high throughput patients' genomic biomarkers such as RNA and DNA sequencing [- Utilizing the transferred tumor pieces from patients to mice, the PDX study [22, 93S. [ aims to personalize the optimal cancer treatment for five types of cancer among a large number of FDA-approved preclinical cancer therapies, given the assigned treatments, the observed genomic features and responses. In particular, among five types of cancer, we focus on the ColoRectal Cancer (CRC), which includes 43 PDX lines with complete genomic information. For CRC, 847 mice expanded from different PDX lines were treated with unique treatments among M = 20 single or combination treatments. The response Y is measured by the scaled maximum observed tumor shrinkage from the baseline time, where the larger value is preferable. After the preprocessing and supervised screening steps shown in [29]] we select 93 significant genomic biomarkers among all features.
Our goal is to identify the possible homo- geneous group structure among these 20 treatments and cluster the treatments with
Boxplot of Response for Treatments in PDX Data
Response Scaled
Single Treatment
Combination Treatment
similar treatment effects in order to boost
the performance the estimated ITR. The
boxplot of the observed response for the 20 treatments in PDX data is shown in Fig- ure D. We can see that, with the highest mean and median of the response, the com-
1.
UC
ENC
LKA
CGM
CET
HDM
LEE
BYL
STU
CLR
BIN
Ciok"
CET CET-ENC
DYL4ENC
BYL BYLICET-ENE
BYL -CE
(ice-l C
BYL LLL M
BYL LBNN
bination treatment BYL+BIN is superior than other treatments. For all methods, we randomly split the data into six folds with 200 replications where five folds are used to train the model in each replication. We use the remaining fold to compute the empirical value based on the estimated ITR. The tuning pa- rameters are selected based on the three-fold cross validation.
Treatment
Figure 3: Observed response in PDX study.
Results for the PDX study. We first include all the single and combination treatments in the model. As shown in the left panel of Figure gl with the highest empirical value, SCAF can provide more effective ITR than other comparison methods. On the overall dataset, we plot the solution path of the hierarchical clustering for all treatments in the left panel of Figure 5L With the best tuned X, it can be seen that the 20 treatments can be clustered into three groups. Most of the combination treat- ments belong to the same group (the green one). One possible explanation is that they include BYL as a common component. It is interesting to point out that the combination treatment BYL+BIN forms treatment group itself. Moreover, SCAF recommends BYL+BIN as the optimal treatment to 92.1% patients among the PDX data, which is consistent with the outstanding performance of BYL+BIN shown in Figure EL This also explains why BYL+BIN itself forms a group because its recommendation proportion nearly dominates the other treatments.
To better examine the group structure of the treatment effect, exclude the combination treatments and only consider the 13 single treatments due to the superior performance of BYL+BIN. Based on the right panel of Figure E the overall values of all four methods decrease compared with the left panel because we drop the combination treatments and reduce the treatment space, while our pro- posed SCAF still has the best performance. Similarly, we draw the solution path of the hierarchical clustering for the single treatments in the right panel of Figure EL The best tuned À further suggests that the single treatments can be clustered into three groups where the second treatment group (the green one) is recommended as the optimal to 94.1% patients in the PDX data. The outstanding performance of the second group matches with Figure because the treatments HDM and LKA in the first group (the pink one) and CGM, ENC and LJC in the third group (the blue one) are the treatments with smallest mean but highest variance among the 13 single treatments. In particular,

---

Treatments
Without Combination Treatments
Response.Scaled
0.5
Response Response.ccaled Scaled
0.5
0.0
0.0
-0.5
-0.5
-1.0
Methods
RF
AD
PLS
SCAF
Figure 4: Buxplots of value the testing data in PDX study.
Figure C.4 in the supplementary materials shows the path of empirical value function of the testing data for both cases. As becomes larger, the empirical value reaches the maximum point due to the helpful clustering structure and finally dramatically decreases because of the excessive merging pattern imposed by the large A. The results in Figure C.4 consistent with Figure
With All Treatments
Without Combination Treatments
0.008
0.015
0.000
0.00
LEE
BYL Livee
BYL PrLICCET
BYL-ENCE
BYL hee
CETHE CETHENC
cole
CLie
Ehice
BYL IaLI
100
SF
D00
BICK
DYL
C
CUR
CG
BI04
Figure 5: Solution paths for the PDX study as shown in Figure A
6 Discussion
In this paper, we consider the challenging ITR problem with a large number of treatments which may have homogeneous group structure. We propose the SCAF method to identify this structure and cluster the treatments in order to learn the optimal ITR more effectively. In particular, we adapt the idea of supervised clustering and formulate the problem convex optimization that consists of penalyy We are able to cluster the treatments and estimate the optimal ITR simultaneously within single optimization problem. The whole clustering process can be intuitively visualized with dendrogram plot similar to the hicrarchical clustering. An efficient algorithm is proposed to solve the problem and the numerical results demonstrate the superior performance SCAF.
The proposed method can be extended several aspects. First, we consider the continuous outcome in this paper. One interesting future direction is to extend SCAF to deal with various types of outcomes, such as discrete outcome and survival outcome. Second, our method can be extended to learn the group structures of treatments for multi-stage dynamic treatment regimes [6, [T]. We leave these interesting directions for future research.
Acknowledgments and Disclosure of Funding
The authors would like to thank the helpful and constructive comments from the reviewers which led to much improved presentation of this paper. Yufeng Liu was supported in part by NSF grant DMS 2100729 and NIH grant R01GM BDGGM225500 26550.. Donglin Zeng was partially supported by NIH R01GM124104, 201MH123499 and RO1NS073671. Any opinions, findings, and conclusions or recommendations expressed in this material are those the authors and do not necessarily reflect the views of National Science Foundation National Institute of Health.
10

---

References
[1] Bibhas Chakraborty and EE Moodie. Statistical methods for dynamic treatment regimes. Springer-Verlag. doi, 10:978-1, 2013.
[2] Francis S Collins and Harold Varmus. A new initiative on precision medicine. New England journal of medicine, 372(9))993-75,, 2015.
[3] Dimitris Bertsimas, Nathan Kallus, Alexander M Weinstein, and Ying Daisy Zhuo Person- alized diabetes management using electronic medical records. Diabetes care, 40(2):210-217, 2017.
[4] Michael R Kosorok and Eric B Laber. Precision medicine. Annual review of statistics and its application, 6:263-286, 2019.
[5] Min Qian and Susan A Murphy, Performance guarantees for individualized treatment rules. Annals of statistics, 39(2): 180, 2011.
[6] Susan A Murphy. Optimal dynamic treatment regimes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65(2):331-335 2003.
[7] Chengehun Shi, Alin Fan, Rui Song, and Wenbin Lu. High-dimensional a-learning for optimal dynamic treatment regimes. Annals of statistics, 46(3):925, 2018.
[8] Michael P Wallace and Erica EM Moodie. Doubly-robust dynamic treatment regimen estima- tion via weighted least squares. Biometrics, 71(3):636-644 644, 2015.
[9] Zhengling Qi, Dacheng Liu, Haoda Fu, and Yufeng Liu. Multi-armed angle-based direct learn- ing for estimating optimal individualized treatment rules with various outcomes. Journal of the American Statistical Association, 115(530)678-691,, 2020.
[10] Yingqi Zhao, Donglin Zeng, A John Rush, and Michael R Kosorok. Estimating individual- ized treatment rules using outcome weighted learning. Journal of the American Statistical Association, 107(499): 100--1118 2012.
[11] Ying-Qi ning-Qiiz Zhao, Donglin Zeng, Eric B Laber, and Michael R Kosorok. New statistical learning methods for estimating optimal dynamic treatment regimes. Journal the Association, 110(510):583-5555 598, 2015.
the American Statistical
[12] Xin Zhou, Nicole Mayer-Hamblett, Umer Khan, and Michael R Kosorok. Residual weighted learning for estimating individualized treatment rules. Journal of the American Statistical Association, 11961177 2017.
[13] Toru Kitagawa and Aleksey Tetenov. Who should be treated? empirical welfare maximization methods for treatment choice. Econometrica, 86(2):591-616, 2018.
[14] Chong Zhang, Jingxiang Chen, Haoda Fu, Xuanyao He, Ying-Qi Zhao, and Yufeng Liu. Mul- ticategory outcome weighted margin-based learning for estimating individualized treatment rules. Statistica sinica, 30:1857, 2020.
[15] Baqun Zhang, Anastasios A Tsiatis, Eric B Laber, and Marie Davidian. A robust method for estimating optimal treatment regimes. Biometrics, 6814)1000- 2012.
[16] Ying Liu, Yuanjia Wang, Michael R Kosorok, Yingqi Zhao, and Donglin Zeng. Augmented outcome-weighted learning for estimating optimal dynamic treatment regimens. Statistics in medicine, 37766)3377--3788, 2018.
[17] Nathan Kallus. Balanced policy evaluation and learning. Advances in neural information processing systems, 31, 2018.
[18] Susan Athey and Stefan Wager. Policy learning with observational data. 89(1):133-1611 2021.
Econometrica,
[19] Yinghao Pan and Ving-Qi Zhao. Improved doubly robust estimation in learning optimal indi- vidualized treatment rules. Journal of the American Statistical Association, 116(533):283-294 2021.
11

---

[20] Christopher John Cornish Hellaby Watkins. Learning from delayed rewards. PhD thesis, Uni- versity of Cambridge, UK, 1989.
[21] Robert Tibshirani, Michael Saunders, Saharon Rosset, Ji Zhu, and Keith Knight. Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society; Series B (Statistical Methodology), (7(11:911000) 2005.
[22] Xiaotong Shen and Hsin-Cheng Huang. Grouping pursuit through a regularization solution surface. Journal of the American Statistical Association, 105(490):27-733, 2010.
[23] Zheng Tracy Ke, Jianqing Fan, and Yichao Wu. Homogeneity pursuit. Journal of the American Statistical Association, (100(509)1155-9909 2015.
[24] Susan Wei and Michael R Kosorok. Latent supervised learning. Journal of the American Statistical Association, 108(503):557-970, 2013.
[25] Juan Shen and Xuming He. Inference for subgroup analysis with structured logistic-norma mixture model. Journal of the American Statistical Association, 110(509):303-333 312, 2015.
[26] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
[27] Shujie Ma and Jian Huang. A concave pairwise fusion approach subgroup analysis. Journal of the American Statistical Association, 1126517:410-4233 2017.
[28] Jingxiang Chen, Quoc Tran-Dinh, Michael R Kosorok, and Yufeng Liu. Identifying heteroge- neous effect using latent supervised clustering with adaptive fusion. Journal of Computational and Graphical Statistics, 30(1):43-55) 2021.
[29] Naim U Rashid, Daniel J Luckett, Jingxiang Chen, Michael T Lawson, Longshaokan Wang, Yunshu Zhang, Eric B Laber, Yufeng Liu, Jen Jen Yeh, Donglin Zeng. et al. digh-dimensional precision medicine from patient-derived xenografts. Journal of the American Statistical Asso- ciation, 1 (6((55)) 11135)):144-1155 140-1154, 2020.
[30] Yu Nesterov. Gradient methods for minimizing composite functions. Mathematical program- ming, 140(1))12--161, 2013.
[31] Donald B Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of educational Psychology, 66(5):688, 1974.
[32] Weibin Mo and Yufeng Liu. Efficient learning of optimal individualized treatment rules for heteroscedastic or misspecified reatment-free effect models. arXiv preprint arXiv:2109.0255 029002570 2021.
[33] Susan A Murphy, Mark J van der Laan, James M Robins, and Conduct Problems Preven- tion Research Group. Marginal mean models for dynamic regimes. Journal of the American Statistical Association, 166455): 2001.
[34] Amir Beck and Marc Teboulle. A fast iterative strinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183-202, 2009.
[35] Martin J Wainwright. High-ddmenssonal statistics: A non-asymptotic viewpoint, volume 48. Cambridge University Press, 2019.
[36] Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American statistical Association, 93(456) 2001.
[37] Hui Zou and Runze Li. One-step sparse estimates in nonconcave penalized likelihood models. Annals of statistics, 36(4):1509, 2008.
[38] Cun-Hui Zhang. Nearly unbiased variable selection under minimax concave penalty. The Annals of statistics, 38(2):894 942, 2010.
[39] Wen Wang, Shihao Wu, Ziwei Zhu, Ling Zhou, and Peter X-K Song. Supervised homogeneity fusion: a combinatorial approach. arXiv preprint arXiv:2201.01036 2022.
12

---

[40] Stefan Wager and Susan Athey, Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association. 113(523):1228-11 12222 2018.
[41] Kornelia Polyak et al. Heterogeneity in breast cancer. The Journal of clinical investigation, 121(10):3786-35) 3788, 2011.
[42] Despina Siolas and Gregory Hanmon. Patient-derived tumor xenografts: transforming clinical samples into mouse models. Cancer research, (1117):3315-3119) 2013.
[43] Manuel Hidalgo, Frederic Amant, Andrew V Biankin, Eva Budinská, Annette T Bymee Carlos Caldas, Robert B Clarke, Steven de Jong, Jos Jonkers, Gunhild Mari Mællandsmo, et al. Patient- derived xenograft models: an emerging platform for translational cancer research. Cancer discovery, (99:999-11113, 2014.
[44] Hui Gao, Joshua M Korn, Stéphane Ferretti, John E Monahan, Youzhen Wang, Mallika Singh, Chao Zhang, Christian Schnell, Guizhi Yang, Yun Zhang, et al. High-throughput screening using patient-derived tumor xenografts to predict clinical trial drug response. Nature medicine, 21(11):1318-1325, 2015.
13

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you descrihe the limitations of your work? [Yes] See Section f
(c) Did you discuss any potential negative societal impacts of your work? [Yes] See Section
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendix B in supplementary materials
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [Yes] See Section and supplementary materials
(b) Did you specify all the training details (c.g., data splits, hyperparameters, how they were chosen)? [Yes] See Sections 2 and B
(c) Did you report error bars (e.g.. with respect to the random seed after running experi- ments multiple times)? [Yes] See Section
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Appendix A in supplementary materials
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes] See Section
(b) Did you mention the license of the assets? [Yes] See Section R
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
(d) Did you discuss whether and how consent was obtained from people whose data you're re using/curating? [Yes] See Section B
(e) Did you discuss whether the data you using/curating contains personally identifi- able information or offensive content? [Yes] See Section I
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
14


Reviews:
[{"review_id": "vGgzTaa1pZi", "reviewer": "Reviewer_cwpD", "summary": "The authors consider the problem of finding optimal individualised treatment\nrules – an important problem in precision medicine – that assigns treatments based\non covariates which may change over time. The particular setting of interest\nis one where many treatments exist and where treatments have a structure induced\nby shared mechanisms of action (e.g., different drugs that target the same pathway).\n\nThe proposed approach uses a fusion penalty term that encourages clustering\nbetween treatments. This incorporates the structure assumption in a tunable way,\nand allows sharing of information between infrequently observed treatments and\nmore frequently observed ones that are related.\n", "questions": "- The PDX study employes a number of preprocessing steps, some of which are\n  supervised. How was bias avoided with the supervised preprocessing? Was\n  separate data used?\n\n- The treatments employed are surely employing drugs with known targets. Do the\n  targets/targeted pathways align with the inferred grouping structure?\n", "limitations": "- Given the strong patient focus in the intro, the confounding limitation should\n  be discussed. The PDX study used in experiments is free of confounding, but\n  typical observational patient data is not.\n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- The literature is well covered in their citations, with the exception of\n  confounding. There are works considering ITRs in the context of confounding\n  through the use of instrumental variables. This is not strictly a problem\n  since the authors aren't considering confounding.\n\n- The work builds on existing with the novelty being the fusion term.\n\n\n- The formulation is straightforward and easy to understand, as is the writing\n  and presentation.\n\n- There are two parts to the method, first there's the group lasso based\n  classification into homogeneous and heterogenous variables, and second is the\n  optimisation of the main objective function with the fusion penalty term. Most\n  of the work concentrates on the second part: the theoretical properties assume\n  the hetrogeneous variables are given and the experiments do not specifically\n  investigate the identification of homo/hetrogenous variables. There is a short\n  note on lines 150–152 stating it has been empirically observed to be more\n  effective, but these are not shown anywhere. Furthermore, it's unclear if the\n  group lasso classification was used in experiments or not, and if it was how\n  the hyperparameter was chosen.\n\n- The theoretical results are relevant and interesting but are limited to the\n  fusion model only.\n\n- The PDX study is not sufficiently described, it is unclear exactly which\n  dataset is used as no reference is provided. There is a citation to another\n  methods paper, cited as the authors have used their preprocessing steps, but\n  no citation to the main data reference.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "h0tDYJjkh6V", "reviewer": "Reviewer_NJYx", "summary": "This paper proposed a method for learning the optimal individualized treatment rule (ITR). It focuses on a situation of having extensive treatment options but limited observations for only a small number of specific treatments. This work designs an algorithm to merge similar treatments and provide optimal ITR in a reduced treatment space. The model is trained based on a designed convex optimization problem with an adaptive proximal gradient algorithm. The empirical study reveals its ability to detect the structure of the treatment space and merge similar ones to enable the optimal ITR from many treatment options. A theoretical analysis is also provided.", "questions": "Referred to Strengths and Weaknesses", "limitations": "Will such a treatment algorithm introduce any bias to the patients?", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The algorithm focuses on a more practical situation of optimal ITR in precision medicine. The author models the problem as a convex minimization problem and solves it with an accelerated gradient method. This algorithm can provide a hierarchical structure of the relationship between treatments. By merging similar treatments, users can give optimal ITR on a reduced treatment space which is more feasible in practice.\n\nI have the following concerns:\n1. It looks like the final result automatically generates the hierarchical structure of the treatments. However, I didn't understand how such a hierarchical structure could happen by just solving the convex minimization problem. Your algorithm seems can only provide the proximity/dissimilarity matrix on treatments. Did you later apply any linkage-based hierarchical clustering algorithms on the proximity/dissimilarity matrix to generate the structure? Please specify.\n2. Is there an analysis of your algorithm's time/space complexity?\n3. Did you compare with other methods that can also provide a dendrogram? \n4. In lines 170-171, should it be if \\hat{\\xi}_k = 0, then Z_k is classified as V? \n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "P09UwO0zaG9", "reviewer": "Reviewer_BQsP", "summary": "The work proposes a machine learning approach for precision medicine applications, where one needs to choose from a set of suitable treatments the best one for a given patient. Further, the proposed approach allows clustering similarly behaving treatments together, and the authors give theoretical guarantees that the underlying clustering will be recovered.\n\nThe proposed approach combines together a number of different approaches including:\n- Basic regression model, where the outcome is modelled as a sum of a main effect M(z), and interaction effect T(z,a) \n- After solving a standard regression problem to recover M(z), group Lasso is used to find the subset of features that may contribute to the interaction effect\n- Proposed adaptive proximal gradient algorithm is used to predict the remaining residual with T(z,a), and a pairwise regularizer is used to divide the treatments into clusters\n\nExperiments on simulated data show that the approach works on finding the true underlying structure on simple simulated problems, and evaluation on a cancer treatment data set shows both competitive performance against basic ML regression approaches, as well as the treatment structure recovered by the method.\n\n", "questions": "- How does the proposed method compare on real data to the baseline that always chooses the treatment that worked best on average on training data, and predicts the mean effect for that treatment? I find it a bit difficult to interpret from the metrics provided, how much is actually learned from the data compared to just predicting the mean.\n- What does the \"response scaled\" -metric mean, what units, scaled how? Maybe something to add to the supplementary?\n- Regarding comparison to the baseline methods, it would be good to tell in the supplementary materials the used parameter grids, and also what kind of feature representations these methods were supplied with. For example, if the goal is to learn both M(z) and T(z,a), in what form was the identity of the treatment a supplied to the baseline methods? The code is provided, but it would be good to be able to tell these details already from reading paper + supplementaries, in order to understand what the comparison means.", "limitations": "Limitations and negative societal impact not considered. The one possible limitation that comes to my mind, is that the method should be much more comprehensively validated on more realistic data sets and settings, before it could be ever considered for making treatment decisions for real-world patients.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper addresses an important problem in precision medicine, and provides ideas that appear to be novel about how to, in addition to predicting, simultaneously cluster together similar treatments. The writing is fairly clear, the work and contributions well motivated, and I found the mathematical presentation to be rigorous yet accessible.\n\nThe proposed learning approach seems reasonable, though to me some of the choices made to make divide this into a sequence of solvable (convex) optimization problems feel a bit like not so well justified \"hacks\" (solving for M(z) and T(z,a) separately, the assumption that Z decomposes nicely to X and V that can be recovered with group Lasso...). Theoretical analysis seems convincing, but I cannot verify that all the details are correct.\n\nThe experimental analysis on the PDX study is interesting, though from results on a single data set it is difficult to ascertain yet how well the method can be in general expected to perform. The one thing that was left a bit unclear to me was does the proposed approach perform more accurate predictions, than the simple approach of doing non-personalized predictions (i.e. always predicting the one treatment that works best on training data, on the benchmark data most likely usually the BYL+BIN treatment). Based on zooming on Figures 3 and 4 this might be the case, but this could be more clearly discussed in the text. \n\nMinor comments:\n\nEquation (1): I was a bit surprised there was not in addition to M(z) and T(z,a) a U(a) type of term that would depend only on the treatment, since in addition to personalized effects one would except some of the treatments to be overall superior to others. E.g. the BYL+BIN treatment in the PDX study. Though I suppose such effect can be encoded into T(z,a) with constant covariate for z, and identity encoding covariate for a...\n\n\"In practice, one may assume that only certain elements of Z...[divide into X and V]\" - this is not intuitively obvious to me, as in many cases one could imagine same covariates contributing both to M and Z? Any plausible real-world examples of this?\n\nI did wonder if some of the reinforcement learning type of terminology used in the beginning (rewards, value functions) were necessary, since in the end the considered setup is regression where training data has been gathered in advance, and gathering feedback and considering exploration/exploitation types of tradeoffs is not considered. The presentation could be simpler by leaving these connections out, but perhaps there is also value in making these connections.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "6PRfNw-4rRG", "reviewer": "Reviewer_DGyw", "summary": "This submission deals with treatment selection among a large collection of candidate ones for a specific individual or a group of individuals which share common response to the treatment. The authors propose to adopt a treatment clustering strategy by gathering treatments with similar effects. \n\nThe proposed approach follows the line of thought of some prior arts in which a statistical regression model is defined for the reward viewed as a function comprising a general and common effect incurred from all treatments plus some treatment specific effects and an additive noise. The common and treatment specific parts of the reward are functions of the patient whose information is embedded in a fixed-size vector. The goal of the authors is to learn the regression coefficients based on supervision coming from data in the form triplets (patient, treatment, reward). They leverage a group promoting penalty to encourage treatments with similar effects to be clustered together. \nTo rule out some of the covariates inside the vectors of patient features, the authors first use a group lasso (on the same reward regression) as a preliminary step. Finally, the main regression problem has properties allowing it to be solved by an accelerated proximal gradient algorithm. \n\nThe paper is technically sound and quite well written. A theorem proves the consistency of the estimator under technical (but not unusual) assumptions.", "questions": "Major : \n\nHave the authors considered investigating more flexible models by using for instance a kernelized version of linear regression ? \n\nI think that the interpretation of random variable $A$ should be clarified at the beginning of section 2. Indeed, should it be understood as the identity of the best treatment, or is it the result of some sampling distribution (e.g. uniform) from which the data is obtained ? Because it is difficult to imagine obtaining data from the former, I suspect the correct interpretation is the latter. The action sampling process in the first experiment of section 5 should be documented. \n\nHow is the term $M_0(.)$ estimated in the real data experiment in section 5 ?\n\nMinor : \n\nFor better readability, I suggest not to use bold letters for scalar variables such as $\\zeta_{a,k}$.", "limitations": "Not adapted to discrete rewards and limited learning capacity.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nPros : \n- the paper incorporates a novel aspect (treatment grouping) in a standard model used for personalized treatment selection\n- the proposed method has a theoretical group-consistency property\n- the efficiency of the approach is experimentally validated on synthetic and real data\n\nCons : \n- the model is designed for continuous outputs whereas discrete ones may be observed for several treatments/pathologies\n- the model backbone is a linear regression and may underfit the data in various situations", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept