Optimal personalized treatment selection with multivariate outcome measures in a multiple treatment case

Abstract In this work we propose a novel method for individualized treatment selection when there are correlated multiple treatment responses. For the K treatment ( ) scenario, we compare quantities that are suitable indexes based on outcome variables for each treatment conditional on patient-specific scores constructed from collected covariate measurements. Our method covers any number of treatments and outcome variables, and it can be applied for a broad set of models. The proposed method uses a rank aggregation technique that takes into account possible correlations among ranked lists to estimate an ordering of treatments based on treatment performance measures such as the smooth conditional mean. The method has the flexibility to incorporate patient and clinician preferences into the optimal treatment decision on an individual case basis. A simulation study demonstrates the performance of the proposed method in finite samples. We also present data analyses using HIV clinical trial data to show the applicability of the proposed procedure for real data.


Introduction
The concept of personalized medicine is fairly old, but the idea advanced dramatically after the introduction of randomized controlled clinical trials that also collected additional patient information.The primary aim of clinical trials is to make population-level decisions which do not necessarily are optimal at an individual patient level.Such population-level decisions do not always account for patient heterogeneity.But the increasing availability of vast amounts of additional patient data from such studies has increased the awareness of heterogeneity in both patient characteristics and outcomes and lead to new evidence-based medicine concepts.Over the last two decades, the statistical methodology in personalized medicine has led to new methodologies and insights, mostly owing to advancement in computational power, bioinformatics discoveries, and access to electronic health data (Murphy 2003;van't Veer and Bernards 2008;Zhao et al. 2012;Vazquez 2013;Kosorok and Moodie 2015).The goal of personalized medicine is to use data to improve decision making in health care to provide the "best" outcome for a patient based on his/ her individualized features.
In real-life situations, the success of treatment may not be fully reflected through a single outcome as a variety of factors may compel both patients and clinicians to consider cure/recovery in a rather broad view.For example, in the treatment for Type-2 Diabetes, control in HbA1c, systolic blood pressure, low-density lipoproteins, cholesterol levels, and prevention from hypoglycemia and weight gain have been suggested as therapeutic goals to address the net clinical utility of a treatment (Unnikrishnan et al. 2013).In cancer studies, although the overall survival is considered the most important, a variety of other factors such as a reduction in tumor size or eradication of cancerous cells, are also considered to be meaningful outcomes (Johnson et al. 2015).Similarly, in situations where the disease is a life-threatening condition, while the time-to-event outcomes (e.g., overall survival) are commonly considered the best outcome, there could be other factors with secondary importance, such as those relating to the quality of life and economic impact.Also, it is common to use a collection of surrogate outcomes during the early development of treatments (Fleming and Powers 2012).Hence, selecting the best treatment considering multiple outcome measures becomes a relevant issue for most patient populations.
Few articles have attempted to address the treatment selection in the face of multiple responses.The creation of a composite outcome via a linear combination of outcomes was advocated by Lizotte, Bowling, and Murphy (2012) In a contrasting approach, Laber, Lizotte, and Ferguson (2014) and Ertefaie et al. (2016) provide continuous/dynamic learning methods for selecting set-valued treatment regimes when there are multiple responses.Rather than a single optimal treatment, the general recommendation in set-valued selections is a set of possible treatments that are "no worse" than any other in that set.When there are multiple treatments, methods proposed in both articles above conduct the selection at each of several stages comparing two treatments at a time.Another article that addresses multi-response treatment selection is by Lizotte and Laber (2016) where these authors focus on a multi-objective sequential optimization method that again gives a non-dominated treatment set that is commonly known as a Pareto optimal collection among possible treatments.On the other hand, Butler et al. (2017) provide a selection method for only two treatments where a single treatment is recommended using patient survey data in addition to clinical data.
To address the multiple response issue, Siriwardhana, Datta, and Kulasekera (2020) used a rank aggregation method to obtain an ordered list of treatments based on multiple ranked lists of treatment performance measures.However, a limitation of that technique is that the method does not fully address relationships among multiple outcomes.In treatment selection problems dealing with multiple responses, outcome measures often tend to be correlated, especially when a collection of clinical or behavioral outcomes are used.For example, scores relating to cognitive improvement in multiple cognitive domains are typically positively associated (Li, Lindenberger, and Sikstr€ om 2001;Cheng et al. 2012).Similarly, when multiple surrogate markers are used on a specific clinical outcome, naturally, such markers are assumed to be correlated.As such, there are many cases where the natural dependency among responses can not be ignored, demanding that treatment selection methods be capable of addressing potential dependencies among responses, in the multiple response setting.
In this paper, we consider the selection of the optimal treatment among K possible treatments for a patient using his or her baseline characteristics when multivariate outcomes (responses) are to be considered.First, to handle statistical issues arising due to high dimensional covariates, each patient is assigned a score based on his/her covariate values.Next, the set of conditional means of multivariate responses are estimated at a given patient score, using a smooth mean estimation step for K treatment options.Finally, a rank aggregation concept that is also capable of handling potential dependencies among ranked lists is applied to find an overall ranking for the K options.The "optimal" treatment for the given score is defined as the treatment estimated to be the best-ranked option in the overall ranked list.The proposed method allows one to use apriori opinions on the importance of each response in determining the best treatment procedure.Empirical studies show that the proposed method has very desirable properties in terms of the selection frequency of the best treatment and an aggregated average gain.The article also demonstrates an application of the proposed technique to a real dataset resulted from an HIV clinical trial.
The remainder of the article is organized as follows.In Section 2, we discuss the proposed methodology.Section 3 includes simulation results followed by a real data illustration in Section 4. The main body of the paper ends with a discussion in Section 5.

Treatment selection
In this section, we describe the proposed procedure on selecting the optimal personalized treatment based on multiple outcome measures.Suppose we observe J continuous response variables for each patient undergoing a treatment selected from K possible treatments and, without loss of generality, suppose larger values of each component of the J dimensional response vector are indicative of better outcomes.Let Y Ã k ¼ ðY Ã 1k , :::, Y Ã Jk Þ 0 indicate the vector of responses for the kth treatment with an associated r dimensional covariate vector X.We assume, there exists a natural correlation structure R k among Y Ã k responses.In this work we assume that we have data from a randomized clinical trial (RCT) study that provides responses and covariate information of a set of patients randomized into K arms.It should be noted that in practice, using the data resulted from a RCT experiment one cannot observe the full J Â K matrix of counterfactuals ðY Ã 1 , :::, Y Ã K Þ for a single patient; hence, one cannot obtain a sample from the joint distribution of ðY Ã 1 , :::, Y Ã K , XÞ: Rather, one observes K independent pairs of observations ðY k , X k Þ from marginal distributions of ðY Ã k , XÞ for k ¼ 1, :::, K, where Y k ¼ ðY 1k , :::, Y Jk Þ 0 (see Siriwardhana et al. (2019, Siriwardhana, Datta, andKulasekera 2020) for a detailed discussion).Further, assume that we can use a patient's covariate value X is to obtain a lower dimensional composite patient score UðXÞ, as described in Siriwardhana, Datta, and Kulasekera (2020) say to summarize each patient's characteristics.
Here, we consider pairs of independent observations ðY k , X k Þ from the marginal distribution of ðY Ã k , XÞ, k ¼ 1, :::, K to select the optimal treatment for K treatments using vectors of smoothed conditional means for each treatment.We define and vectors l k ðuÞ ¼ ðl 1k ðu 1 Þ, :::, l Jk ðu J ÞÞ 0 for U ¼ ðU 1 , :::, U j Þ 0 and u ¼ ðu 1 , :::, u J Þ 0 where components of these vectors correspond to each 1, :::, J response.Although we suppress the dependence of U j s on the covariate vector for brevity, in all developments below, quantities related to U j s are functions of X.
In our proposed approach, we rank the K values for each component of l k ðuÞ vectors (k ¼ 1, :::, K) to get size K vectors v j ðuÞ ¼ ðv j1 ðuÞ, :::, v jK ðuÞÞ 0 , where v jk ðuÞ is the rank of l jk k ¼ 1, :::, K among l j ðuÞ for each j (here j ¼ 1, ::, J) with the largest l jk ðu j Þ value given the rank 1.We next produce an overall rank by following an aggregation technique as a basis to find the optimal treatment.
In their previous approach Siriwardhana, Datta, and Kulasekera (2020) used a aggregation method by Pihur, Datta, andDatta (2007, 2009) to combine these rank vectors v j ðuÞ; j ¼ 1, :::, J to get an overall ranking of treatments v Ã ðuÞ ¼ ðv Ã 1 , :::v Ã K Þ for a given patient score u: They defined the optimal treatment as k Ã ðuÞ ¼ arg min However, the rank aggregation method by Pihur, Datta, andDatta (2007, 2009) does not account for dependencies among the J rank lists, which could potentially lead to unwanted estimation errors when such dependencies exist.In many real life problems dealing with multiple responses, natural associations among responses can not be ignored.In finite sample cases, unignorable dependencies among estimated response means due to correlated responses, subsequently translate to correlated rank lists when one uses rankings of those estimated means.In the current work we propose to use an aggregation technique that can be utilized under correlated rank lists for the proposed treatment selection concept.We detail the proposed aggregation step in the sequel.
We link the jth component Y jk of the response vector Y k for the kth treatment and covariates X k via a Single Index Model (SIM).The SIM formulation provides great flexibility and reasonable efficiency in modeling many types of data.This model is expressed as, for j ¼ 1, :::, J and k ¼ 1, :::, K, where each b jk is a r-vector of parameters, g jk is an unknown smooth link function and jk s are error terms with E½ jk jX ¼ 0: We assume independence of jk s across k ¼ 1, :::, K for a fixed j where these terms are correlated across js for any given k.
RCT data used in our approach are of the form ðY ki , X ki Þ where Y ki ¼ ðY 1ki , :::, Y Jki Þ 0 and Y jki indicates the jth component of the response for the ith individual under treatment k with associated covariate values X ki , i ¼ 1, :::, n k : The relationship (3) between response and covariates for such a sample can be written as Following Siriwardhana et al. (2019) and Siriwardhana, Datta, and Kulasekera (2020) we define a score vector UðXÞ for a patient with covariate X as follows.First define, Next, define the jth component of the combined overall score vector as The overall score is given as UðXÞ ¼ ðU 1 ðXÞ, :::, U J ðXÞÞ 0 where U j ðXÞ ¼ ðS j ðXÞ, d j ðXÞÞ 0 for j ¼ 1, :::, J: Since model functions defined in (3) contain unknown parameters, components of these score vectors should be estimated using a standard function estimation method.In the literature, there are many different estimation techniques available for estimating the link function and the index vector of a SIM, allowing us to use one out of several available reasonable estimation methodd to estimate the gs and the b s (Ichimura, Hall, and Hardle 1993;Hristache et al. 2001;Yu and David Ruppert 2002).We adopt the Hristache et al. (2001) procedure in our simulations and data analysis in the sequel.
In particular, for any given vector X ¼ x, let and Û j ðxÞ ¼ ð Ŝj ðxÞ, dj ðxÞÞ 0 ; j ¼ 1, :::, J: As suggested in Siriwardhana, Datta, and Kulasekera (2020), a suitable estimator for l jk ðu j Þ, k ¼ 1, :::, K at a given u j ¼ ðs j , d j Þ 0 can be obtained using the smooth mean estimator given by, where w is a kernel function with x !0 and Ð xðtÞdt ¼ 1, and h jk , k ¼ 1, :::, K are a set of smoothing parameters.Here I(A) is the indicator of A. The bandwidth selection for estimating l ki s is a challenging issue.However, as Siriwardhana, Datta, and Kulasekera (2020) suggested, methods given in Wand and Jones (1995) for kernel smoothing provides a reasonable solution for this estimation problem.
First, find the distance c j ðv, û0 Þ ¼ gðv, v j ðû 0 ÞÞ, for a given v 2 V K , a rank list of length K, where V K is all permutations of integers f1, ::, Kg: Here we propose to use Spearman's rank distance (Pihur, Datta, and Datta 2007) for gð:Þ, which is defined as where M 1 , :::, M K is a list of real values, r 1 , :::, r K are the ranks of M k , k ¼ 1, :::, K and v 1 , :::, v K is a permutation of integers 1, :::, K and q is a positive number.
Let Cðv, û0 Þ ¼ ðc 1 ðv, û0 Þ, :::, c J ðv, û0 ÞÞ 0 for a fixed rank vector v: Now, we define an overall rank distance Lðv, û0 Þ by where Dðv, û0 Þ is a suitable dispersion matrix for Cðv, û0 Þ: Here K is a diagonal matrix with diagonal elements ðs 1 , :::, s J Þ which signify the practical importance of 1, ::, J responses including the views of patients and clinicians.As such, Lðv, û0 Þ can be considered as a Mahalanobis type weighted distance from Cðv, û0 Þ to the origin IR J þ : We propose to minimize Lðv, û0 Þ, with respect to v 2 V K , for an estimated score û0 corresponding to a new patient.Suppose this minimum occurs at a vector We then define the optimal treatment as We approximate Dðv, û0 Þ by a suitable weighted dispersion matrix Dðv, û0 Þ, calculated from Cðv, ûki Þ s using lj ðû jki Þ s, i ¼ 1, :::, n k ; k ¼ 1, :::, K, corresponding to training samples from each treatment group.The weights W ki s, say, for i ¼ 1, :::, n k ; k ¼ 1, :::, K, for this calculation are developed using a localization approach to be described below.For a given v 2 V K we find where mð:Þ is the weighted average given by mðCðv, ûÞÞ ¼ In developing Dðv, û0 Þ, we use a localizing weighting scheme to localize Cðv, ûki Þ vectors at û0 ¼ ðû 10 , :::, ûJ0 Þ 0 using K training samples of size n k each in order to achieve a reasonable approximation for Dðv, û0 Þ: Our approach is the following.For the ith training observation in the kth group with a score value ûjki ¼ ðŝ jki , djki Þ, a weight is defined as hjk using the bandwidth hjk and kernel xð:Þ for j ¼ 1, :::, J and k ¼ 1, :::, K: Next, define an overall weight as the product of J individual weights For this purpose, we propose to use the same bandwidth h jk and kernel xð:Þ, that were used in the calculation of ljk ðû jki Þ s, corresponding to the ith training patient in the kth set; i ¼ 1, :::, n k , k ¼ 1, :::, K: Even though such a selection of weights may not be the optimal, as observed from empirical results, these weights seem to be reasonable and flexible to implement the proposed procedure.
In summary, the proposed approach can be implemented as follows: Using the historical RCT data (training set), estimate the set of J Â K single index models that are relating covariates and J response variables corresponding to each treatment k, k ¼ 1, :::, K (i.e., model ( 4)).
Estimate scores for all patients in the training data set and the new patient.Use the criteria given in formula (6) to estimate the jth sub-component of the score with respect to jth response ûj : Find the estimated overall score û ¼ ðû 1 , :::, ûJ Þ in each case.

Empirical studies
In this section we present a simulation study that investigates the properties of the proposed procedure in finite samples.
We performed a series of simulations with the proposed procedure under various settings to evaluate its performance.Primarily, we focused on the accuracy of treatment assignment of a new (test) observation using estimated values of l jk functions from a set of training data.This simulation study was performed for K ¼ 2, 3 treatments with response dimension J ¼ 2, 3, and 4. To illustrate the personalized medicine treatment concept, in our simulations we selected our model sets such that each model in a set dominates other competing models for some combination of covariate values.In particular, none of the considered models fully dominate other models within the whole covariate space.Hence, subjects with distinct covariates vectors could experience the highest response from different treatments.
In our study, we first simulated K independent multivariate (dimension J) samples of size n (n ¼ 50 or n ¼ 100) per group.The components of the r dimensional covariate vectors X were generated independently from the UðÀ1, 1Þ distribution, where r was fixed at 10. Using various link functions and index vectors, we obtained the treatment responses from model (3) for each k.We examined the performance of the proposed methodology under a set of highly nonlinear regression models given by sine/cosine function that followed the SIM structure.Such trigonometric functions pose most difficulties in smooth estimation and therefore we believe these models somewhat present worst case scenarios with respect to estimation of the optimal treatment.A set of model functions used for the study is provided in Table 1.For a given k, k ¼ 1, :::, K, the errors were generated from either a J dimensional multivariate normal distribution or a J dimensional multivariate double exponential distribution with zero mean and a compound symmetric correlation matrix where the off-diagonal values were chosen from the set f0:1, 0:5, 0:9g: The R package mvtnorm (Genz et al. 2017) was used for the generation of these random vectors in the multivariate normal case, where the dispersion parameter r N was chosen from the set f0:1, 0:3, 0:5g: We used the R package LaplacesDemon (Statisticat, LLC 2017) for generating Double Exponential random variables with dispersion parameter r D chosen from f0:1, 0:3, 0:5g: Once the K samples were generated, we estimated the corresponding SIMs followed by an estimation of scores at each covariate value.SIMs were estimated by the procedure given in Hristache et al. (2001) using Epanechnikov kernels.Then, a new covariate value X 0 was generated in the same manner as above, and for its corresponding estimated score û0 , we calculated ljk ðû j0 Þ for k ¼ 1, :::, K; j ¼ 1, :::, J: Similarly, we calculated ljk ðû jki Þ s for the ith patient in the kth training set; i ¼ 1, :::, n k ; k ¼ 1, :::, K, for dimension j, j ¼ 1, :::, J: Next we obtained rank distances c j ðv, ûki Þ ¼ gðv, v j ðû jki ÞÞ for a given a rank list v 2 V K , where V K is the set of all permutations of the integers f1, :::, Kg and gð:Þ is the Spearman's footrule distance function with q ¼ 1.This produced a vector of rank distances CðvÞ ¼ ðc 1 ðv, ûki Þ, :::, c J ðv, ûki ÞÞ 0 : Next, following the procedure in (11), a localized dispersion matrix Dðv, û0 Þ, at the neighborhood of û0 ¼ ðû 10 , :::, ûJ0 Þ was obtained as an approximation for Dðv, û0 Þ, based on the K training sets of size n k each, k ¼ 1, :::, K: This was followed by the calculation of Lðv, û0 Þ given in (8) for a chosen K matrices and we then use the proposed procedure given in (10) to estimate corresponding kÃ ðû 0 Þ: The kernel function in this estimation was taken to be a Normal (0, 1) probability density function.We chose all bandwidths by the algorithm given by Wand and Jones Wand and Jones (1995) for each h jk , k ¼ 1, :::, K; j ¼ 1, :::, J: Table 1.Sets of smooth mean functions used for generating treatment responses.We choose the common vector C to be a unit vector; C ¼ ð 1 ffiffiffi ffi 10 p , ::: 1 ffiffiffi ffi 10 p Þ 0 , for all combinations of j and k.
Next, to define the "correct selection" for the above covariate value X 0 , we follow the following approach.First, we generated K new response vectors, Ỹk0 ¼ ð Ỹ 1k0 , :::, Ỹ Jk0 Þ 0 , k ¼ 1, :::, K, each with mean vector ðg 1k ðb 0 1k X 0 Þ, :::, g Jk ðb 0 Jk X 0 ÞÞ 0 for k ¼ 1, :::, K, corresponding to this X 0 using model (3) where the errors were generated independently from the same error distribution that was used to generate the K original samples.Our approach will be to define quantities similar to those given in ( 8), ( 9) and ( 10) using these new observation vectors.
To do that, we ranked rows of the J Â K matrix Ỹ 0 created from response vectors Ỹk0 ; among Ỹ j10 , :::, Ỹ jK0 for each j, j ¼ 1, ::, J: Then we proceeded to the calculation of rank distance c j Ỹ0 ðvÞ ¼ gðv, v j ð Ỹ0 ÞÞ and the corresponding distance vector C Ỹ0 ðvÞ ¼ ðc 1 Ỹ0 ðvÞ, :::, c J Ỹ0 ðvÞÞ 0 for a fixed v: Next we obtained a suitable dispersion matrix, say D Ỹ0 ðvÞ) for C Ỹ0 ðvÞ as follows.

Ỹ0
: The treatment assignment for the new patient was considered to be correct if kÃ ¼ k: We repeated this procedure 1000 times for each model and error distribution combination.Frequencies of correct treatment assignments for a representative set of cases are given in the Tables 2 and 3, and Supplementary Tables S1-S6.
Simulations results demonstrate reasonable selection accuracies in the each scenario considered.More importantly, the selection frequency remained consistent at large values of the response correlation, indicating the potential of the proposed technique for such cases.As to be responses.Now, letting C be the theoretical maximum for any k 1K ðuÞ within the whole covariate domain, we average k 12 ðû 0 Þ=C and k 1K ðû 0 Þ=C for the 1000 new test cases and denote them by X 12 and X 1K , respectively.
Note that, measures given by Xs quantify the potential gains/losses in terms of the weighted aggregation of mean outcomes for borderline cases (i.e., mis-classifications close to the decision boundary) as well as cases that are furthest from the optimal case according to the proposed procedure.We report a few results in Tables 4 and 5. Positive values of X 12 and X 1K indicate aggregated relative average gains in expected treatment outcomes by our treatment selection technique.As noted in tables, all X values are positive and they are higher for X 1K cases than those for corresponding X 12 s indicating that our proposed procedure results gains with respect to average responses.

ACTG-175 HIV clinical trial
In this section, we provide another illustration of our proposed method using, the data resulted from the ACTG 175 Clinical Trial (Hammer et al. 1996).This clinical trial was a randomized, double-blinded, placebo-controlled clinical trial that was conducted to compare single nucleoside or two nucleosides antiviral medications in adults infected with human immunodeficiency (HIV-1) whose T-cell CD4 counts were in the range of 200 to 500 per cubic millimeter.The study randomized HIV-1-infected patients to one of four daily regimens: 600 mg of zidovudine (arm-0), 600 mg of zidovudine plus 400 mg of didanosine (arm-1), 600 mg of zidovudine plus 2.25 mg of zalcitabine (arm-2), or 400 mg of didanosine (arm-3).The data set by Juraska et al. (2017) contains information on 2,136 HIV-1-infected subjects.Arms 0, 1, 2, and 3 contain 532, 519, 524, and 561 patients, respectively.
Our intention in this data analysis is merely to demonstrate what would be the optimal treatment for a new subject if one were to use the proposed treatment selection based on individual patient characteristics when training samples with corresponding covariate values for multiple treatments are available in advance.In this illustration, rather than assuming there is a standard treatment (i.e., control) that is being compared with experimental treatment(s) as in a typical clinical trial, we take the stance that given multiple treatments can be used to treat a patient, what would be optimal for the individual based on his/her characteristics.In this study, subjects were examined periodically, capturing their T-cells counts (i.e., CD4 T helper cells and CD8 cytotoxic T cells), that are critical components in the human immune system.The scientific literature on HIV/AIDS often declares CD4 cells as the primary T-cell type that is suppressing the HIV cell Table 4. X 12 and X 13 values calculated using 1000 test cases by the proposed method.Three treatments with four responses, using equal weights.
Error correlation ðqÞ q ¼ 0 q ¼ 0:5 q ¼ 0:9 Error dist.replication, by signaling various other cells for immune response.For an HIV patient, the severity of the disease progression is directly measured by the decline in CD4 counts.The important role of the CD8 cell is typically referred to as the antibody reaction against cancers and various types of other viruses.However, some studies have illustrated the important role of the CD8 during the early stages of HIV progression (e.g., see, Streeck and Nixon 2010).
In previous a study, Siriwardhana, Datta, and Kulasekera (2020) used this dataset to demonstrate personalized treatment strategies with respect to log-transformed dual outcomes given by CD4 and CD8 counts of a patient after 20 weeks.However, their analysis did not account for correlation among the two responses, whereas the current work has the capability to handle such dependencies.We estimated the Pearson product-moment correlation between CD4 and CD8 responses as 0.251 using the observed data.In this analysis also we used log transformed CD4 and CD8 counts after 20 weeks as our responses.As covariates, we used log-CD4 and log-CD8 counts at baseline, age, weight, and the number of months a patient received the pre-antiviral therapy.
As a training set, we selected 200 patients at random from each group and considered their data as outcomes of a RCT trial in which those patients were randomly assigned to one of the four treatments.We used the training data to estimate the corresponding SIMs for two outcomes: log -CD4, and log -CD8, together with above covariates.Considering the remaining 1336 cases as "new" patients, we applied the proposed treatment selection for those cases.In this illustration, we ranged the "importance" weights from 0 to 1 for each response.This procedure was repeated for 1000 random training (i.e., size ¼ 200, per group) and testing (i.e., size ¼ 1336) sets to obtain average assignments for those 1000 partitions.
Error correlation ðqÞ q ¼ 0 q ¼ 0:5 q ¼ 0:9 Error dist.64.50%, 15.48%, and 16.76% test patients are proposed to be assigned to Arms 0, 1, 2, and 3, respectively; whereas the corresponding assignment for weights 0.5:0.5 (equal weights) are 5.28%, 61.41%, 14.31%, and 19.00%.The overall pattern of these assignments indicates that only a few patients are proposed to be assigned to zidovudine alone, which was used as the control arm in the study (arm-0).Our analysis using individual characteristics shows that Treatment Arm 1 appears to be a choice for most patients while the control Arm is the least chosen.It is interesting that the marginal analysis of ACTG-175 data (Hammer et al. 1996) concluded, that treatment Arms 1, 2 and, 3 all slowed the progression of HIV disease at a higher level compared to Arm 0 although no Arm among 1, 2 and 3 was selected over others.Proposed assignments in our approach appear to have classified many patients for Arm 1 compared with others.
Extending our analysis we conducted the following study to examine the effectiveness of proposed treatment selection for test patients compared to their original random assignment.
Suppose CD4 and CD8 responses indicated by j ¼ 1 and j ¼ 2, respectively.For a given testing set, consider the tth test subject with a covariate value X t who was randomly assigned to a particular arm in the original study.Suppose the individual's estimated score is ût ¼ ðû 1t , û2t Þ, where û1t ¼ ðs 1 , d 1 Þ and û2t ¼ ðs 2 , d 2 Þ are marginal estimated scores based on CD4 and CD8 responses respectively.In our approach, we attempted estimating the gain in conditional means of log À CD4, and log-CD8 values for an individual, comparing the proposed assignment versus his/her original assignment.Let v, v ¼ 0, ::, 3, be the treatment option suggested by the proposed technique for the test patient based on his/her estimated score ût and let v 0 , v 0 ¼ 0, ::, 3 be the treatment group patient was assigned in the original trial by the random assignment.We define the average conditional gain D j for a given j, j ¼ 1, 2, as where Y jv 0 is the observed value by the original assignment.Following, Siriwardhana, Datta, and Kulasekera (2020), we estimated EðY jv jû t Þ using a smooth mean estimator given by, where n v is the number of patients originally assigned to treatment arm v, v ¼ 0, ::, 3, w is a kernel function with x !0 and Ð xðtÞdt ¼ 1, and h jv 's, j ¼ 1, 2 are corresponding smoothing parameters.We proceeded the estimation using Gaussian kernels and the corresponding smoothing parameters were determined by the method suggested in Wand and Jones (1995) for kernel smoothing.Next, we estimated the overall gains Dj , j ¼ 1, 2 values by averaging rescaled Djt values, seemingly suggesting that if one were able to use prior data to estimate the personal score for a new patient, then the proposed assignment is more beneficial for him/her than an assignment based on a clinical trial which may only find the best treatment for the general population.

Discussion
In this article, we proposed a novel personalized treatment plan to select the optimal treatment from a set of multiple treatments when the outcome measures are multivariate.This is an extension of the method by Siriwardhana, Datta, and Kulasekera (2020), which uses the rank aggregation.Our method assumes that there are historical RCT data that contains patient characteristics and response variables with respect to each treatment, and we use such information to select the best treatment option for a new patient.The proposed method incorporates potential dependencies among responses as opposed to other existing methods in the literature.Our empirical studies show that the new method performs very satisfactorily in selecting the optimal treatment in a multiple treatment setting.Our analysis of a real clinical trials dataset which has multiple treatment options reveals possible changes if one were to use multiple outcome measures as opposed to a single measure.
As the working model, we use the semiparametric Single Index Models to relate responses and subject covariates.This model offers great flexibility with modeling smooth complex relationships, allowing us to handle many real-life situations.The proposed method can also be applied using quantile regression single index models to handle problems with many different error types, deviating from the Gaussian error structure.
In our approach, we defined W ki s similar to the geometric means of w jki s.In a few empirical studies performed under arithmetic means, we observed comparable performances similar to current selection for W ki s.In our empirical study and real data analysis, we used s ¼ 1 and Spearman's distance function.However, the choice of ss and the choice of the distance function gð:Þ can be influential in the optimal treatment selection.Although the specification of response weight is an important step, we consider this as more of a subjective/qualitative issue than a quantitative issue, where one can guide the choice of optimal weights based on patient qualities and preferences and advice of physicians.When dealing with responses that are related to factors that govern the cure from the disease or degradation of the patient due to various aspects of the disease, then a team of clinicians may be better suited to determine what weights should be applied.On the other hand, for outcome measures that correspond to the quality of life, behavior, financial impact, etc., a physician with consultation from relevant support/advisory groups may be more appropriate.For situations where sufficient data is available from various studies, weight selection using more of a data-based method might be appropriate.However, we have not investigated those aspects in the current work.
There can be several extensions of this work.The proposed method requires a collection of complete data records from RCT studies, but when dealing with life-threatening or terminal conditions, priority is typically given to survival type measures that could be subject to censoring (e.g., overall survival, disease-free survival, etc.).We have not explored that avenue in this work.Also, although we focused only on continuous outcomes, binary and ordinary outcomes are also common (e.g., relapse, response level, cure, etc.) and therefore selection methods for such responses would be valuable.In addition, another exploration is possible modifications needed in order to use the proposed method with data from other sources such as observational studies.

Table 2 .
Accuracies of treatment selection in 1000 test cases using the proposed technique for the case of three treatments and four responses (equally weighted).

Table 6 .
Treatment assignment summary for ACTG-175 Clinical Trial data, by the proposed method selecting both CD4 and CD8 counts as clinical response, using weights s CD4 and s CD8 for CD4 and CD8 counts, respectively.

Table 7 .
DCD4 proposed method compared to the original assignment for test patients, using weights s CD4 and s CD8 for CD4 and CD8 counts, respectively.The 10-fold cross-validation technique was used for estimating DCD4 ¼ 1336 is the number of test patients and S Dj is the standard deviation of Djt s.It is rea- c and DCD8 c by the