High dimensional precision medicine from patient-derived xenografts

12/13/2019 ∙ by Naim U. Rashid, et al. ∙ 0

The complexity of human cancer often results in significant heterogeneity in response to treatment. Precision medicine offers potential to improve patient outcomes by leveraging this heterogeneity. Individualized treatment rules (ITRs) formalize precision medicine as maps from the patient covariate space into the space of allowable treatments. The optimal ITR is that which maximizes the mean of a clinical outcome in a population of interest. Patient-derived xenograft (PDX) studies permit the evaluation of multiple treatments within a single tumor and thus are ideally suited for estimating optimal ITRs. PDX data are characterized by correlated outcomes, a high-dimensional feature space, and a large number of treatments. Existing methods for estimating optimal ITRs do not take advantage of the unique structure of PDX data or handle the associated challenges well. In this paper, we explore machine learning methods for estimating optimal ITRs from PDX data. We analyze data from a large PDX study to identify biomarkers that are informative for developing personalized treatment recommendations in multiple cancers. We estimate optimal ITRs using regression-based approaches such as Q-learning and direct search methods such as outcome weighted learning. Finally, we implement a superlearner approach to combine a set of estimated ITRs and show that the resulting ITR performs better than any of the input ITRs, mitigating uncertainty regarding user choice of any particular ITR estimation methodology. Our results indicate that PDX data are a valuable resource for developing individualized treatment strategies in oncology.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 39

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The complexity of human cancer is reflected in the molecular and phenotypic diversity of patient tumors (Polyak, 2011). This diversity results in heterogeneity in response to treatment, which complicates clinical decision making. The complexity of human cancer is also reflected in the high failure rate of new therapies entering oncology clinical trials, highlighting limitations in the ability of standard preclinical models to evaluate new therapies (Tentler et al., 2012). A recent study utilized patient-derived xenografts (PDXs) to perform a large-scale screening in mice to evaluate a large number of FDA-approved and preclinical cancer therapies (Gao et al., 2015). Genomic information and observed treatment responses were used to identify efficacious therapies that standard cell line model systems had missed, and also validate known associations between genomic biomarkers and differential response to treatment. The results from this PDX study mirrored those seen in human patients. Thus, PDX models can be used to evaluate in vivo therapeutic response and discover novel biomarkers to inform individualized treatment decisions.

PDX models are based on the transfer of primary human tumors directly from the patient into an immunodeficient mouse (Siolas and Hannon, 2013). Briefly, pieces of primary solid tumors are collected from patients by surgery or biopsy (Hidalgo et al., 2014). The collected tumor pieces from an individual patient are then implanted into mice subcutaneously to create a PDX line, whereby tumor size and rate of tumor growth after implantation may be measured over time. After the tumor reaches sufficient size, the line may be expanded by passaging directly from the implanted tumor into additional genetically identical mice. Through this expansion, multiple treatments may be applied to mice originating from the same PDX line, allowing for the application of multiple treatments to the same patient tumor. High throughput genomic assays such as RNA sequencing (RNA-seq) and DNA sequencing may be performed on the original patient tumor. Features of the original tumor will be retained throughout line expansions (Hidalgo et al., 2014), making PDX models ideal for learning how to personalize cancer treatment, given observed feature-response associations.

Personalized treatment recommendations in oncology have traditionally centered around the classification of patients into subgroups (Sargent et al., 2005). In some cases, patient subgroups may be derived from predictive models based upon genomic biomarkers (Parker et al., 2009). Yet, significant heterogeneity in treatment response may still be observed within such subgroups (Metzger-Filho et al., 2012; Chen et al., 2014), and assignment of optimal treatment is predicated upon accurate subgroup prediction. An alternative approach to precision medicine is the estimation of an individualized treatment rule (ITR), a map directly from the patient covariate space to the space of allowable treatments that can be used to make treatment decisions. The optimal ITR is defined as the one that maximizes the mean of a clinical outcome, such as treatment response, when applied in a population. Examples of such covariates may include patient clinical information, such as laboratory assay results, or high dimensional genomic data, such as gene expression or mutation data from a patient’s tumor. As such, treatment recommendations based upon an optimal ITR may result in improved clinical outcomes by harnessing individual-specific molecular and clinical features not captured by subgroup-based approaches.

A number of methods have been proposed to estimate an optimal ITR. One approach is to fit a regression model for treatment response given a set of applied treatments and patient covariate information. The optimal treatment is then the one providing the maximum predicted response in a new patient, given the fitted regression model and the covariate information for that patient (Qian and Murphy, 2011). An example of this approach is Q-learning (Murphy, 2005; Zhao et al., 2009; Qian and Murphy, 2011; Schulte et al., 2014). Direct search methods, including outcome weighted learning (OWL) (Zhao et al., 2012; Zhou et al., 2017; Liu et al., 2016; Chen et al., 2017), doubly robust ITR estimators (Zhang et al., 2012, 2013), and marginal structural models (Robins et al., 2008; Orellana et al., 2010)

estimate the optimal ITR using inverse probability weighting rather than regression. In direct search methods, the class of ITRs is prespecified, while in other approaches, the class of ITRs is implied by the modeling process. Recent advances in machine learning for causal inference have produced a number of estimators for the conditional average treatment effect that could be used to estimate an ITR for a binary treatment decision

(Imai et al., 2013; Athey and Imbens, 2016; Wager and Athey, 2018). Other methods for estimating optimal ITRs include marginal mean models (Murphy et al., 2001) and Bayesian predictive modeling (Ma et al., 2016). Regardless of the approach, many existing methods may directly utilize high-dimensional genomic biomarker data. However, such methods were not designed for PDX studies, where the application of multiple treatments within a subject (PDX line) results in correlated outcomes, and the number of available treatments is large.

In this manuscript, we utilize the wealth of clinical and biomarker data generated by the Novartis PDX study (Gao et al., 2015) to estimate optimal ITRs for treating several human cancers. Given the correlated outcomes within PDX lines, the large number of treatments, and the high dimension of the covariate space, it is difficult to fit a nonparametric model to the conditional mean response without large amounts of data. Thus, we explore a number of ways of imposing structure on the conditional mean response, including reducing the dimension of the covariate space, grouping treatments that result in a similar mean response, and constructing a treatment tree to convert the problem of selecting from a large set of treatments to a sequence of binary comparisons. The result is a multi-step procedure, where each step can be thought of as imposing structure on the model for the conditional mean response in a way that alleviates the challenges posed by PDX data and takes advantage of the unique structure of PDX data. Such multi-step procedures are well-studied and have shown good performance in many applications (Love et al., 2014; Bourgon et al., 2010; Soneson and Delorenzi, 2013; Rashid et al., 2014). We show that the proposed multi-step procedure achieves improved performance over standard methods in certain settings. We examine the various modeling decisions that are made at each step of the multi-step procedure and demonstrate the use of super-learning (Luedtke and van der Laan, 2016) to improve prediction performance by aggregating the proposed models.

In Section 2, we describe the large-scale PDX data set that motivated this work and use it to highlight the challenges associated with estimating optimal ITRs using PDX data. We present our methodological approaches in Section 3. We present results from our data analyses in Section 4. In Section 5, we conclude with a discussion to compare and contrast the various modeling decisions we made and discuss the clinical implications of our findings. Additional details, including software for reproducing our work, are given in the Supplemental Materials.

2 Large-Scale PDX Drug Screen for Treatment Response

2.1 Data Overview

The Novartis PDX study (Gao et al., 2015) established a total of 1075 PDX lines corresponding to a variety of human cancers. A subset of these lines were genomically profiled prior to treatment for gene expression (399 lines), copy number analysis (375), and mutations (399). In addition, 281 lines were enrolled in the drug response study. Our study utilizes 190 PDX lines with complete genomic and response data (Supplementary Figure 1). Five types of cancer are represented among these 190 lines (Supplementary Figure 2): breast cancer (BRCA), melanoma (CM), colorectal cancer (CRC), non-small cell lung cancer (NSCLC), and pancreatic cancer (PDAC). A median of 21 treatments were applied per PDX line across different cancers (Supplementary Figure 2). Certain PDX lines had fewer mice and therefore fewer treatments than the total number of treatments available for a particular cancer type. One mouse per PDX line was set as a control line and did not receive treatment.

In total, 3487 mice, expanded from the 190 PDX lines with complete data, were used for our study. Details regarding the biomarker data utilized in this study are given in Section 6 of the Supplemental Materials. Briefly, each of the 190 PDX lines had 22,665 genes measured for gene expression via RNA-seq, 23,854 genes measured for gene-level copy number estimates via copy number array, and between 159 and 293 mutations (25 and 75 percentiles) identified via DNA sequencing. In total, a union set of approximately 60,000 features are available for ITR estimation. Because each genomic assay was performed on the patient tumor prior to implantation, all mice expanded from the same PDX line share the same set of genomic biomarkers.

2.2 Study Design and Response Variables

Mice from each PDX line were treated with either single agents or combinations. A total of 38 unique therapies were applied, administered either as a single agent (36 total administered) or in combination with other agents (26 total combinations administered). Certain treatments were limited to particular cancers, whereas others were applied across cancers. Each mouse was treated for a minimum of 21 days. Tumor size was evaluated twice weekly by caliper measurements, and the approximate volume of the tumor was calculated using the formula (, where is the major tumor axis and is the minor tumor axis.

Two measures were utilized to summarize response to treatment: best average response (BAR) and time to tumor doubling (TTD). BAR is defined as , where is the day in which the th measurement was taken, is tumor volume at day 0, and is tumor volume at measurement taken on day (Gao et al., 2015). BAR is a measure of the maximum observed tumor shrinkage from baseline over measurements taken at least 10 days after start of treatment, scaled by time since baseline. More negative values of BAR indicate a better response. We used BAR in the analysis so that larger values indicate a better response. BAR mirrors similar criteria for assessing response in human clinical trials (RECIST, Therasse et al., 2000)

. TTD is defined as the number of days from baseline that the tumor doubled in size from its baseline measurement. Due to skewness in the distribution of TTD, we used the natural log for analysis.

2.3 Implications of PDX Data for ITR Estimation

The unique structure of PDX data—the application of multiple treatments to mice implanted with the same tumor—makes PDX data ideally suited for estimating optimal ITRs. Comparing responses between mice within the same line does not amount to observing true treatment effects due to the inherent variability that exists across mice; however, the improved precision of observing responses to multiple treatments applied to the same tumor may substantially improve the performance of estimated ITRs. The method we propose involves a sequence of decisions between two groups of treatments. At each step, we aim to choose the group that contains the best treatment; thus, the proposed method directly uses multiple responses within each line that would not be available without PDX data.

Existing methods for estimating optimal ITRs have typically been designed for a small number of treatments (two, for example). In this case, the set of treatments is large (20 or greater). Modeling the conditional response is difficult in the presence of a large set of treatments due to the large number of terms that would be included in the model. In a preclinical study like this, it is more important to identify groups of promising treatments than to fully assess all pairwise comparisons between treatments. Thus, we reduce the size of the treatment set by grouping treatments with similar mean responses using hierarchical clustering. This allows us to adaptively create treatment groups which have approximately similar treatment effects.

Finally, the set of genomic biomarkers is high-dimensional. Approximately 60,000 genomic features are available, and many of the available features may exhibit either low variability or low expression across samples. While some methods for estimating an optimal ITR can handle such high dimension, we implement some common preprocessing steps similar to existing statistical methods in bioinformatics to reduce the dimension prior to ITR estimation. We also utilize some additional steps to mitigate the ultra-high dimension of the predictor space, which we motivate and justify in the next section.

3 Methods

3.1 Overview

Let denote response and let

denote a vector of covariates. Let

denote treatment, where exactly one of is equal to 1 and the rest equal to 0, with indicating that treatment is given. The conditional mean response can be expressed as

(1)

If we were to obtain estimators , , an estimator for the optimal ITR would be . However, fitting model (1) nonparametrically in the PDX setting would be difficult without large amounts of data, due to the large number of treatments and high dimension of the covariate space. Therefore, we propose to impose structure on model (1) to ameliorate these difficulties. The result is a multi-step procedure to reduce the dimension of the feature space, reduce the size of the treatment set, and estimate optimal ITRs for the reduced treatment set using the reduced feature space. The steps of our procedure are described in Sections 3.1.1-3.1.3 below. Various modeling decisions must be made at each step, and some alternative choices for these decisions are discussed in Section 3.3. Because it is not immediately obvious what the optimal set of modeling decisions is, we apply super-learning (Luedtke and van der Laan, 2016) to combine variants of the proposed method with different embedded models. For complete details on our methodology, we refer the reader to Section 7 of the Supplemental Materials.

3.1.1 Preprocessing (step 1).

In the first step, we perform a preprocessing of genomic features to remove features without sufficient variance (Supplemental Materials Section 7.1). Genomic features are filtered out at this step due to low expression or low variability. This is a well-studied technique for dimension reduction prior to analysis

(Love et al., 2014; Bourgon et al., 2010; Soneson and Delorenzi, 2013; Rashid et al., 2014). This screening is performed separately for each cancer type. A summary of the number of genes and features remaining after this step is given in Supplementary Table 2. In addition, treatments applied in less than 90% of PDX lines within a cancer type were filtered out for that cancer type (see Section 7.1 of the Supplemental Materials). We summarize the number of treatments applied per PDX line per cancer after treatment filtering in Supplementary Table 1.

3.1.2 Supervised screening (step 2).

In the second step, we further reduce the dimension of the feature space by selecting likely prognostic and predictive genomic features (Supplemental Materials Section 7.2). Genes are ranked using Brownian distance covariance (Székely et al., 2009), evaluating the dependence between a vector of gene-level predictors for each PDX line and the bivariate response (BAR, TTD) for a given treatment (prognostic) or difference in response between a pair of treatments (predictive). After ranking genes, we select the top genes (, Supplementary Table 2) and use all available platforms for the selected genes, giving us features corresponding to each value of . Screening in this way imposes structure on model (1) by forcing the to be constant in some of the covariates for certain . This strategy helps to alleviate the difficulties caused by the large number of genomic features and is similar in nature to sure independence screening (Fan and Lv, 2008), used in ultra-high-dimensional regression problems to reduce the feature dimension to a more moderate size prior to application of variable selection methods. All of our analyses are repeated for each of top genes to yield insights into the performance of estimated ITRs based on differing numbers of features. Screening is performed separately for each cancer type.

Because the set of genomic features consists of multiple platforms per gene, the feature set resulting from selecting the top genes may still contain a large number of features. In addition, predictors generated from different platforms (e.g., RNA-seq and copy number) on the same gene may be correlated, and gene expression may be correlated across different genes. Thus, there may be a lower-dimensional feature space that would be sufficient for ITR estimation. To address this, we applied a further dimension reduction step using deep learning autoencoders (DAE, Wang and Laber, 2017)

, a variant of deep neural networks

(Vincent et al., 2010). See Section 7.3 of the Supplemental Materials for details. This allows us to evaluate whether a lower dimensional representation of the feature set improves ITR estimation at the cost of additional computational burden.

Deep learning autoencoders build a nonlinear prediction model to predict all feature variables from a low-dimensional representation of the original features, with dimension chosen by cross-validation. The reconstruction error of this approach is several orders of magnitude lower than principal components analysis across cancer types, particularly for

(Supplementary Tables 3 and 4). This indicates that linear dimension reduction techniques are not sufficient to capture the information contained in the data. Dimension reduction using DAEs imposes structure on model (1) by forcing the , to depend on only through the low-dimensional summary of the original features. The dimensions for each feature set following application of deep learning autoencoders within each cancer are given in Supplementary Table 5. After applying steps 1 and 2 to obtain low-dimensional sets of covariates, we use the low-dimensional sets of covariates to estimate optimal ITRs.

3.2 Estimation of Treatment Tree-based ITRs

Estimating an ITR to select from a large number of treatments is challenging for two reasons. First, fitting model (1) in the presence of a large number of treatments would be difficult due to the large number of treatment feature interaction terms. Second, the resulting ITR would be difficult to interpret and implement. Instead of an ITR which selects one treatment from the full set of treatments, it would be more useful to have an ITR which selects a group of treatments from a partition of the full set of treatments, where treatments in the same group can be expected to lead to similar responses. Ideally, we would create groups of treatments such that the conditional mean response function is the same for all treatments in the same group. However, this treatment grouping is unknown. In Section 3.2.1

, we describe a heuristic technique for approximating this treatment grouping using hierarchical clustering. This allows us to create different groupings of treatments by cutting the tree (the dendogram resulting from the clustering) in different places. We then define a class of ITRs that can be expressed as a sequence of binary decisions, one for each step of the tree. This tree structure is distinct from tree-based regimes which use a tree structure to represent the final estimated ITR

(Laber and Zhao, 2015; Zhang et al., 2015). A variety of methods could be used to construct decision rules for each step of the tree. We introduce two: a regression-based approach (Section 3.2.3) and a direct search approach (Section 3.2.4).

3.2.1 Treatment tree construction.

Denote the mouse corresponding to the PDX line within the cancer type with subscript , where , , and , with representing the number of treatments applied to PDX line in cancer type , and representing the number of PDX lines for cancer type (Supplementary Figure 2). We let , i.e., for each PDX line, up to mice were expanded to receive a maximum of treatments in cancer type , one treatment applied per mouse (Supplementary Figure 3). For certain lines, treatments may have been applied, as the number of mice per line varied. For ease of notation, we assume that and that the th mouse in each PDX line received the same treatment within cancer type.

For treatment in cancer , we define the treatment response vector , where

is scaled to have standard deviation 1 (rows in Supplementary Figure 6, left panel). Because there are inherent baseline differences in response between PDX lines, we first center the response values of each PDX line (columns in Supplementary Figure 6, left panel) with respect to the “null” response within each PDX line. Using PDX lines with

mice, we calculate the Euclidean distances between and , , to get a distance between each pair of treatments. For a fixed constant, , we group the nearest neighbors of the “untreated” treatment response vector to form a “null” set of treatments, denoted by , containing those treatments producing low or no response. Then, for each , we compute the vector of centered observed outcomes, , where is the vector of averaged treatment responses for each of the treatments in . Thus, is the difference between the response to treatment and the average response in the null treatment group.

Next, we perform hierarchical clustering on the centered treatment response vectors and take the resulting dendogram as a treatment tree (Supplementary Figure 6, right panel). For a fixed constant, , we can construct treatment groups by cutting the tree steps from the root node. For cancer type , we label the treatment groups determined by and as and denote the set of all possible treatment groups in cancer given by . Although depends on and , we will write in place of to simplify notation. The hierarchical clustering groups treatments together when they result in similar mean responses, in contrast to grouping treatments by predefined characteristics such as molecular target or mechanism of action. Grouping treatments in this way forces certain , in model (1) to be equal, similar to a fusion penalty (Tibshirani et al., 2005). This strategy helps to alleviate the difficulties caused by the large number of treatments and allows the estimated ITRs to select a group of treatments likely to produce similar outcomes (Wu, 2016). We average responses within the resulting treatment groups, yielding one response for each PDX line and treatment group in . An ITR can be estimated for each value of and , and the optimal values of and can be selected using cross-validation.

3.2.2 Identification of the optimal ITR.

Let be the vector of genomic features, where we write to make it clear that the domain of the feature space is data-dependent (see Sections 3.1.1 and 3.1.2). An ITR is a mapping . For each treatment , define the potential outcome to be the outcome that would be observed under treatment (Rubin, 1978). Within cancer type , the set of potential outcomes consists of , , . We make the assumption that for all , all , and all . That is, we assume that the expected value of the potential outcomes for mouse and from the same PDX line would be equal if they both received the same treatment. The standard assumption of positivity is not needed because each PDX line is assigned to every treatment used for that cancer type. On the other hand, while the standard assumption of no unmeasured confounders is still needed in our setting, the primary source of confouding is due to the assignment of mice to treatment within PDX line and not due to the genetic features of the PDX line since each line receives all treatments. Hence we assume that the process of assigning treatment to mice is exchangeable and homogeneous, and thus the no unmeasured confounding assumption obtains. Define the value of an ITR, , by . Let be the set of ITRs which can be expressed as a sequence of decision rules, one for each step of the treatment tree, starting from the root node and proceeding until a leaf node is reached. The optimal ITR associated with the treatment tree is , i.e., the ITR that maximizes value within the class.

3.2.3 Treatment tree-based Q-learning.

In this section, we propose an extension of Q-learning (Murphy, 2005; Zhao et al., 2009, 2011; Schulte et al., 2014) to estimate the optimal ITR associated with the treatment tree. Let , denote the set of treatment groups downstream to the left () and right () arms of a node at step , where . Starting at the step corresponding to the lowest node, we compute for and . Here is the mean observed reward (centered response) in PDX line across the treatments belonging to treatment group . We then fit a regression model for based on , separately for and , to obtain , the estimated conditional mean reward in treatment group at step given genomic features. If either or are missing for a given line, then that line did not receive that particular treatment; these observations are removed before fitting the regression model. The estimated optimal decision rule at step of the tree for an individual with genomic features is given by

(2)

i.e., the treatment group with the highest predicted clinical response given is the estimated optimal treatment group at step . We repeat the above process for step , except only utilizes the observations for each PDX line from the optimal treatment group selected at the previous step. Thus, the proposed method directly utilizes the multiple responses observed for each PDX line.

Step corresponds to the highest node in the treatment tree consisting of the non-null treatments, and step corresponds to the split between the null group (centered by their own means) and the non-null treatments from the treatment tree. The decision rule is determined in a similar manner at step , evaluating whether any non-null treatment should be applied. At , the response vector for the null group is a vector of zeros, and the decision rule at this step simplifies to determining whether the expected value of the response under the optimal non-null treatment is greater than 0. A number of techniques could be used to fit the regression models at each stage, representing different ways of imposing structure on the , in model (1). We discuss various embedded models that could be used in more detail in Section 3.3.

Given the sequence of estimated decision rules , the optimal treatment for a new individual given their set of predictors is obtained by following the decision rules sequentially from step downward until one arrives at a terminal node on the tree. Rather than recommending a sequence of treatments, as in standard Q-learning, the estimated optimal decision rule at each step of the tree directs the user through either the left or right arm at each node, creating a path through the tree to arrive at a terminal node representing the optimal treatment group for that individual.

3.2.4 Treatment tree-based outcome weighted learning.

Outcome weighted learning (OWL) estimates the optimal ITR for selecting between two treatments by maximizing an inverse probability weighted estimator of the value function over a fixed class of decision rules. Thus, unlike Q-learning, OWL does not rely a fitted regression model. Our extension of OWL adopts the same tree structure as in the previous section. However, instead of fitting a separate regression model for each arm at split , the optimal decision rule at split is estimated directly using weighted classification methods (Zhou et al., 2017; Liu et al., 2016; Chen et al., 2017). When treatment is binary, any decision rule can be written as for some decision function, . We will assume that the decision boundary at each split is smooth. Thus, at each split, we use a class of decision rules defined by some class of smooth functions , e.g., a reproducing kernel Hilbert space. At the th split of the treatment tree, let , where

(3)

Here, denotes the empirical measure of the data used in the split, is a penalty term for the decision function , is a tuning parameter which is selected using cross-validation, and is a space of functions. As in Q-learning, the observations included when computing the minimizer in equation (3) are all of those corresponding to the estimated optimal treatment group at the previous step. Thus, our extension of OWL also directly utilizes the multiple outcomes observed per PDX line. We discuss options for selecting in Section 3.3; different options for impose different specific forms that the , must take in model (1).

3.2.5 Theoretical justification of tree-based ITRs.

In this section, we show that the value of the optimal ITR associated with a treatment tree increases strictly with each step down the tree, until the step where there is no heterogeneity in conditional mean response for treatments within the same group. Let , , be a grouping of treatments, that is, is a partition of the set of all treatments. Conditioning hereafter on the selected set of features , let be an ITR, i.e., implies that treatment group is selected for a subject with features . We assume that if , the decision maker will select one of the treatments from with equal probability. Define the value of with respect to to be

where is the number of treatments in . The optimal ITR is given by

and the value of the optimal ITR (also called the optimal value) is

We call a grouping (partition) finer than if any can be written as the union of sets in . We obtain the following result:

Theorem 3.1.

If is finer than , then . Furthermore, equality holds if and only if, for any and such that , with probability one.

Proof.

For with , we have for any that

The above holds with equality if and only if is the same for all , or equivalently, for all , with probability one. ∎

Theorem 3.1 states that a finer partition of treatments never leads to a decrease in value of the optimal ITR and that the change in value of the optimal ITR is zero only when the conditional mean response is constant across treatments within groups in the finer partition.

Theorem 3.2.

For a grouping and any set , further partition into two nonempty groups to obtain a finer partition, . If for any such finer partition, , , it holds that

with probability one for all .

Proof.

From Theorem 3.1, we obtain that for any and such that ,

with probability one. Suppose that the proposition does not hold. Then, there exists some such that for all for some with positive probability. For any , order the such that are nondecreasing with at least two adjacent values different, say positions and in the ordered sequence. Then, if we let and , we obtain that

for all . Because has positive probability, this creates a contradiction. ∎

Theorem 3.2 implies that if no binary splits at one branch will lead to an increase in the optimal value, then all treatments within that branch are homogeneous with respect to conditional mean response. Thus, each further partition created by cutting the tree at a lower step will lead to an increase in the optimal value function until there is no heterogeneity in response across treatments in the same group.

3.3 Modeling Decisions

The framework proposed in this paper involves creating a tree of treatments using hierarchical clustering and estimating an ITR as a sequence of binary decisions, one for each step of the tree. Within this general framework, there are a number of specific modeling decisions that need to be made. In this section, we describe different variants of the proposed method that result from making different modeling decisions.

The tree-based ITR estimation method involves an embedded method to estimate the optimal decision rule at each step. Fitting a regression model at each step results in an extension of Q-learning (see Section 3.2.3

). A number of regression models could be used as the embedded regression model; in our analyses, we used linear models with a LASSO penalty and random forests. In Section 

4, these are referred to as QL and , respectively. At each split, our extension of Q-learning uses the maximum outcome across treatments downstream to the left and right of the split as the observed outcomes when fitting the regression model. Alternatively, we could obtain predicted maximum outcomes using the regression model and use the predicted maximum outcomes to fit the model at the next step. These are analogous to the pseudo-values used in standard versions of Q-learning (Zhao et al., 2009, 2011), and would allow application of the proposed tree-based Q-learning method to data that do not come from a PDX study. We examine both strategies in Section 4, and we refer to Q-learning with observed outcomes and with pseudo-values as and , respectively.

Constructing an inverse probability weighted estimator (IPWE) at each step, rather than fitting a regression model (as in Q-learning), results in an extension of OWL (see Section 3.2.4). The IPWE is maximized over a class of functions. We use both the class of linear functions and the reproducing kernel Hilbert space associated with the Gaussian kernel function (Zhao et al., 2012; Chen et al., 2017). In Section 4, these are referred to as and , respectively.

The second modeling decision that must be made is the selection of the covariate set. We applied the proposed method using the dimensional set of genomic features resulting from screening using Brownian distance covariance, for each

(Supplementary Table 1). We also applied the proposed method to the lower dimensional set of features extracted from

genes using the DAE, for each (Supplementary Table 5). Analyses using the features extracted from the DAE are labeled with the subscript “dl" in Section 4.

A final variant of the method that we explored involves replacing the observed outcomes with model-predicted outcomes before estimation. We fit a random forest to predict outcomes based on covariates alone (features and treatments) and replaced the observed outcomes with predictions based on the fitted model for all later stages of the analysis. This approach acts to denoise the observed outcomes. Analyses utilizing this approach are labeled with the subscript “smoothed". In Section 4, we report analyses using different combinations of the modeling decisions described above to capture synergistic effects of the various modeling decisions. In addition to estimating tree-based ITRs as proposed in Section 3.2.1, we also estimated ITRs by fitting model 1 using linear models with the LASSO penalty and random forests. These “off-the-shelf" methods were included to compare to the proposed method.

3.3.1 Super-learning.

The various modeling decisions outlined in Section 3.3 above result in many variants of the proposed method. While certain variants may work better than others in certain settings, the optimal choice of modeling decisions may not always be clear. A natural analysis to try in this case is to combine a set of input ITRs estimated using various embedded models in the hopes that the resulting ITR performs better than any of the input ITRs. To accomplish this, we apply the super-learning algorithm of Luedtke and van der Laan (2016). Briefly, Luedtke and van der Laan (2016) propose combining a number of existing methods using cross-validation to calculate a linear combination of latent functions to maximize the value function. In our implementation, there is not an explicit latent function due to the fact that we are using a sequential treatment tree as our model. To approximate the value of the latent function with respect to a single treatment, we utilize the predicted reward from one of our sub-models at a given node in the tree whose direct children include our goal treatment. To optimize the superlearner, we use simulated annealing to estimate the coefficients. Multiple chains are used whose starting points are selected from a set of randomly generated coefficients. In Section 4, variants of super-learning with different sets of input ITRs are referred to with the subscript sl, followed by the number of ITRs that are combined to produce the superlearner.

3.4 Performance Measures for Individualized Treatment Rules

We used five fold cross-validation to evaluate the performance of ITRs estimated using the proposed method with the various modeling decisions as described in Section 3.3. Within each cancer type, we divided PDX lines into five folds. An optimal ITR was estimated using the training data set that results from holding out each fold, and the estimated value was calculated on the held-out fold as , where denotes the empirical measure taken over mice in the held-out fold. Tuning parameters (such as and for estimating tree-based ITRs and the OWL penalty parameter) were selected using cross-validation within each training data set. The estimated value of the estimated ITR, denoted by is calculated by averaging the value estimates resulting from each fold. We also calculated the standard deviation of value estimates across folds.

Different cancer types may result in different marginal mean outcomes. To facilitate comparisons across cancer types, we also computed the observed value, , defined as the sample average of centered responses for all mice that received non-null treatments, and the optimal value, , defined as the sample average across PDX lines of the maximum centered response across treatments, i.e., . The observed value and optimal value can be used to define two metrics for evaluating estimated ITRs: proportion of optimal value, defined as , and ratio to observed value, defined as .

4 Results

All analyses were performed separately for each cancer type. We applied various combinations of the modeling decisions outlined in Section 3.3. Each modeling variant was applied to the feature set resulting from screening with different values of (see Section 3.1.2), utilizing the original features or the features extracted from the DAE. We present results for BAR and defer results for TTD to Section 5 of the Supplemental Materials.

Figure 1 illustrates the strong linear relationship between the mean value under each estimated ITR and the optimal value for the associated cancer type and treatment grouping. Note that in the majority of cases the optimal values do not vary within a cancer type. However, in some cases there is variability in the optimal value due to the fact that we select and separately for each analysis. Each point in Figure 1 represents a particular estimated ITR within a particular cancer. We also observe a similar relationship between the mean and observed values for each estimated ITR. This suggests that the marginal mean outcome differs across cancer types. This is not unexpected, as some cancers are known to be more sensitive to available treatment options (e.g., CM), and others less so (e.g., PDAC). This observation motivates our use of and to evaluate performance of estimated ITRs. We present results for here and defer results for to Section 4 of the Supplemental Materials.

Figure 1: Original (left) and scaled (right) values (for best average response) from all analyses, across methods, cancer types, and . Optimal values for each method vary significantly by cancer type. The value of each estimated ITR is correlated with the optimal value across cancer types. We normalize the estimated value for each method by the optimal value to allow for comparisons between cancers. The metric, called “proportion of optimal," provides a measure of how close the value of an estimated ITR is to the optimal value.

4.1 Relative Performance of Methods Pooling Across Conditions

We summarize for each variant of the proposed methods and the “off-the-shelf" methods in Figure 2, pooling results across different cancer types and values of . From this figure, several striking trends are apparent. For example, the application of data smoothing prior to ITR estimation boosts the overall performance for many methods, such as Q-Learning with embedded linear models (light and dark green), and OWL methods (light and dark grey). This pre-smoothing is less beneficial for Q-Learning methods utilizing embedded random forests (salmon and red), as the pre-smoothing itself is performed using random forests. Q-learning methods using pseudo-values (red, ) performed similarly to their counterparts (pink) across conditions. In general, Q-learning performed better than OWL overall across variable conditions. In addition, using the lower dimensional features extracted from the DAE did not show significant benefit for most approaches with the exception of the OWL methods using the linear kernel, which we will show later to be sensitive to the dimension of the feature space. Q-Learning methods with non-linear embedded models (salmon, red) showed much better robustness to various modeling choices than linear ones, and showed to be especially helpful in OWL (dark grey). Lastly, simpler off-the-shelf methods showed lower performance and much higher variability across conditions compared to methods utilizing the treatment tree approach. Relative to the LASSO, almost all methods utilizing the treatment tree performed better than the simpler off the shelf methods.

We also find that utilizing a weighted combination of ITRs using the superlearner approach (dark blue) resulted in the best overall performance across conditions. Increasing the number of ITRs included in the superlearner had the effect of boosting performance while also reducing variability in performance across conditions. Here, the SL4 combined all four Q-Learning methods from Fig 2 utilizing pre-smoothing, SL6 includes the addition of methods and , SL8 adds and , and SL16 add all Q-learning methods utilizing DAE features. OWL methods were excluded from the superlearner due to computational time, and we expect additional performance increases through their inclusion. This suggests that the superlearner approach can boost performance and also mitigate user uncertainty regarding the selection of various modeling approaches in cases where the optimal approach may not be clear beforehand.

The examination of pooled overall and relative performance in Figure 2 is informative when the analyst is unsure what approach is best for their data and wants to use the “safest” approach. For example, the results in Figure 2 can be used to determine a sequence of modeling decisions that exhibited good overall performance across a variety of conditions and achieved low variability across conditions. For example, , or Q-learning utilizing smoothed outcomes and observed-values, achieved good performance on average and low variability in performance across conditions. Alternatively, one may elect to combine ITRs from multiple methods using the superlearner approach to avoid choosing an individual method, at the cost of additional computational burden. In addition, our results suggest the following general conclusions: linear methods may benefit from the use of data smoothing prior to ITR estimation, further dimension reduction of the predictor space via DAEs is beneficial for OWL methods using the linear kernel, and nonlinear embedded models are more robust to various modeling choices than linear embedded models.

Figure 2: Overall performance across methods, pooled over cancer types and number of features (top). for each method is also normalized to the LASSO in each condition to highlight the relative performance of each approach (bottom). This relative measure was constructed by subtracting the pertaining to the LASSO from that of the other methods within each combination of cancer type and value.

4.2 Impact of on Performance

We now examine the relative performance of our estimated ITRs by , the number of genes used in estimation. Figure 3 contains boxplots of pooled over cancer types and stratified by estimation method and (Supplementary Table 1). Prior to pooling we adjust in each cancer type and method by the performance at to more clearly delineate changes with respect to dimension. Most methods did not show strong trends in performance across , with the exception of the class of methods (light grey). In this class, decreases with increasing numbers of genes. However, the same trend did not appear when utilizing the Gaussian kernel. Slight downward trends were also observed for Q-learning with an embedded linear model (light and dark green). These observations together suggest that methods in which a linear decision rule is estimated at each step of the treatment tree may be sensitive to the dimension of the feature space.

The optimal set of features for ITR estimation (indexed by ) is the set of all those features for which at least one of the , in model (1) is not constant. While the optimal set of features is unknown, the results in Figure 3

indicate that the proposed treatment tree-based ITR estimation method performs well regardless of the number of features selected, with the exception of linear OWL and Q-learning with an embedded linear model. These results lead to the conclusion that a smaller set of features post-supervised screening may be optimal for methods with an embedded linear model, and an embedded nonlinear model can be used to provide robustness against selecting the “wrong” set of features.

Figure 3: Overall trends in performance for each method across , aggregated over cancer types.

4.3 Overall Performance by Cancer

The performance resulting from specific modeling decisions varies across cancer types. Figure 4 contains boxplots of for the best performing set of modeling decisions within the classes defined by the colors in Figure 2. To select the best performing variant in each class, we chose the one with the largest averaged across values of within that particular class. The boxplots in Figure 4 contain over values of . We show the full set of results corresponding to all methods and cancers in the Section 3 of the Supplemental Materials. Within certain cancers, such as breast cancer (BRCA), specific modeling decisions do not result in large differences in performance. This may be due in part to the fact that a small number of treatments appear to work well uniformly across samples in BRCA (Supplementary Figure 8). In contrast, for pancreatic cancer (PDAC) and non-small cell lung cancer (NSCLC), greater heterogeneity in response exists across treatments (Supplementary Figures 9 and 10). In addition, we find that in almost all cancer types, the superlearner tended to perform better than or similar to all other classes of methods, suggesting its use when it is unclear which individual method may be optimal for a particular dataset.

Table 1 lists the best overall set of modeling decisions for each cancer type, along with the genomic features that were most important for selecting treatments using the best performing estimated ITR. For cancers where the best performing ITR resulted from the DAE predictors or OWL methods using the Gaussian kernel, it is difficult to determine which genes were the most important. For these cases, we selected important features using the second-best ITR (Reference Method in Table 1). For all cancer types, the best performing ITR resulted from the tree-based approach rather than an “off-the-shelf" method. In addition, despite their sensitivity to dimension, the class of methods were represented as the best ITR in two out of five cancers. This suggests that as long as the correct feature dimension is selected beforehand, OWL method can perform well relative to other methods. In practice however, the optimal dimension is difficult to ascertain unless one evaluates multiple candidate feature sets, as we have in this manuscript.

The treatments most frequently recommended by the best performing ITR for each cancer, along with the corresponding values of and , are given in Supplementary Table 6. The treatment tree for the optimal ITR varied in the values of across cancer types, suggesting variability in the amount of response heterogeneity across cancer types. For example, in CM, where response tended to be strong overall (Figure 1), , suggesting that relatively more treatments had similar response profiles across PDX lines. The selected value of for the best performing ITR was low for each cancer type, indicating that only a small number of treatments were found to be effectively the same as “untreated" based on the hierarchical clustering (Supplementary Table 6). For “off-the-shelf" methods, and by design since no grouping of treatments was performed.

We list the average (unconditional) response for each of the treatments within cancer type in Supplementary Table 8, calculated as the sample average response across mice treated with each of the treatments across PDX lines. In BRCA, the treatment with the larger mean response was also the most recommended treatment (LEE011 + everolimus). In PDAC, however, there is less variability in average response across treatments, and the best performing ITR recommends BKM120 + binimetinib to 18 PDX lines and abraxane + gemcitabine to 18 PDX lines (see Supplementary Table 7).

Figure 4: Performance of the best estimated ITR in each class across cancer types.
Cancer Method Lsup Reference Method Top 5 Predictors
BRCA owllinearsmoothdl 50 qlrf2 COL1A1.rna,CTCFL.rna,FMNL3.rna,SRPX.rna,HMCN1.cn
CM owllineardl 100 ql1smooth ETV7.rna,DPYSL3.rna,LEPREL1.rna,GABRE.cn,ATP2B1.rna
CRC ql1smooth 100 ql1smooth FGG.rna,ALPK1.rna,WDR27.cn,DIDO1.mut,C10orf26.rna
NSCLC owlkernelsmooth 100 ql2smooth POFUT2.rna,TNNI3.cn,TUBG1.cn,NAT8L.cn,PTPRE.cn
PDAC ql1smooth 100 ql1smooth CTH.rna,DUSP4.cn,TPP2.rna,ACVR1B.rna,ZNF264.cn
Table 1: Best performing method and number of predictors for each cancer. Top 5 most important predictors are listed. Predictors pertaining to the next best performing method were provided if the top performing method utilized deep learning (Reference Method). Predictors ending in .rna are from the gene expression data, .cn from the copy number data, and .mut from the mutation dataset.

4.4 ITR Performance When Limiting Features to a Single Genomic Platform

The three genomic platforms utilized in this study provide a wealth of information for ITR estimation and biomarker discovery. However, the high dimension of the feature space provides practical and computational challenges. Predictors from the same gene may be correlated (for example, gene expression and gene copy number) and may therefore be redundant. Evaluating three genomic platforms for each PDX line’s original tumor increases cost and amount of tumor tissue required. For these reasons, we also explored the performance of ITRs estimated using only the RNA-seq gene expression platform, a common genomic assay performed by biomedical researchers. We repeated the same process described in Section 3, but using only features resulting from RNA-seq.

The overall conclusions are largely similar to those in Figure 2 (see Supplementary Figure 11). When we compare the difference in between the ITR estimated using the full feature set and the ITR estimated using RNA-seq only, we find that most methods perform similarly using only RNA-seq data (Supplementary Figure 12). Q-learning and OWL methods with embedded linear models did slightly better than corresponding ITRs trained on all three platforms, which is likely related to the relative sensitivity of these methods to the dimensions of the feature space (smaller in the RNA-seq only analysis). Otherwise, most methods evaluated in RNA-seq only performed similarly to ITRs trained on all three platforms. Due to the smaller feature space in the RNA-seq only analysis we did not elect to run the deep learning variants of each method. Overall, these results suggest that utilizing a single genomic platform may be a more cost-effective option that will result in estimated ITRs with comparable performance.

4.5 Top Genomic Features in NSCLC

Lung cancer is one of the leading causes of cancer related deaths in the United States, and NSCLC accounts for the majority of clinical cases of lung cancer (Ettinger et al., 2010). Therapeutic agents used to treat NSCLC include paclitaxel, which interferes with cellular microtubular dynamics and cellular division through the targeting of tubulin (Wise et al., 2000), and cetuximab, which targets the epidermal growth factor receptor (EGFR) (Pirker et al., 2009). Since the mechanism of action differs between these two treatments, it is not surprising that we observed significant heterogeneity in response between these two treatments in this study (Supplementary Figure 10). Here, we examine the relationship between response and the genomic features found to be the most important for making decisions using the best performing estimated ITR.

The best performing ITR for NSCLC resulted from (see Table 1). Given that that gaussian kernel for OWL does not allow direct interpretation of its predictors, we utilize the reference method in this cancer to examine the role of each selected predictor with respect to response. We determined the most important genomic features for this ITR using the following approach. For each of the splits in the associated treatment tree, we computed the absolute value of the regression coefficient for each feature in the model fit at a given node and retained the feature with the largest magnitude value at each node. Cross-validation selected ; therefore, we retained a set of 13 top features, where 10 of these were unique (Figure 5). We then calculated Spearman’s rank correlation coefficient between each selected feature and response, separately within each non-null treatment. This resulted in the matrix of feature-treatment pairwise correlation coefficients seen in Figure 5. We then clustered the columns of this matrix, pertaining to the seven selected genomic features, using Euclidean distance between the vector of correlation coefficients in each column. We also clustered the rows of this matrix, pertaining to the non-null treatments, using the tree structure that was previously constructed for , rather than the correlations.

The treatment with one of the strongest correlation to Tubulin Gamma 1 (TUBG1) copy number is the therapy paclitaxel, suggesting that a higher copy number of TUBG1 may potentiate response in patients being targeted with agents such as paclitaxel. This is notable because paclitaxel directly targets tubulin (Kumar, 1981), of which TUBG1 plays a major role. Furthermore, cetuximab is the only treatment that is anticorrelated with TUBG1 copy number. This unique relationship with TUBG1 is reflected in the treatment tree, as cetuximab is the only member of a branch furthest away from all other non-null treatments. Cetuximab also exhibits a unique mechanism of action, as the only monoclonal antibody EGFR inhibitor in our data set. The relationship between TUBG1 copy number and response to cetuximab and paclitaxel is displayed in Figure 5. These results indicate that that the observed correlations between treatment response and the top features determined from reflect the role these features play as important variables in the best performing ITR.

Figure 5: Top genomic features selected by in NSCLC. Cells are colored by the magnitude of their Spearman’s correlation between response to treatment (rows) and expression of top features (columns). Genomic features were clustered by their vector of Spearman’s correlation to each treatment (using Euclidean distance), whereas treatments were grouped based on the treatment tree constructed for . The patterns of treatment-feature correlations tended to coincide with the predetermined treatment groupings.

5 Discussion

In this paper, we introduced several approaches to ITR estimation using PDX data. The unique structure of a PDX study, where multiple treatments are applied to samples from the same human tumor implanted into mice, naturally lends itself to precision medicine. The substantially improved precision that results from the PDX structure may result in better performing ITRs. However, PDX data also pose a number of challenges, including a large number of unordered treatments and a high-dimensional feature space. These factors make it difficult to nonparametrically model the conditional mean of the response. The method we propose involves of sequence of steps that alleviate these difficulties. Our method involves screening the covariates to find a lower dimensional feature space, constructing a treatment tree that allows for recasting ITR estimation as a sequence of binary decisions, and estimating a decision rule at each split to select the arm that contains the optimal treatment. Because we aim to select the arm that contains the optimal treatment at each split rather than the arm with the largest average response, our estimation technique utilizes the maximum response downstream of each arm for each line. Thus, our estimation technique incorporates the unique structure of the PDX data by directly using the multiple responses observed per PDX line. We’ve shown that the method we propose not only produces high-quality ITRs, but also identifies genes that are known to be associated with response to treatment (e.g., the TUBG1 gene shown in Figure 5).

The methods we propose requires making a number of modeling decisions at various stages of the pipeline, including selecting the dimension of the feature space, choosing embedded models, and selecting tuning parameters, among others. We demonstrated various combinations of these modeling decisions in our analyses. While certain variants performed better than others for certain cancer types, the treatment tree-based approach outperformed “off-the-shelf" methods overall. Reducing noise by using random forest-predicted outcomes (smoothing) improved performance of the estimated ITRs in most settings, and basing ITRs on DAE-extracted features improved performance in the presence of a linear embedded model. In particular, Q-learning with smoothed outcomes and pseudo-values performed well across all settings and is a good first choice for estimating ITRs from PDX data. We recognize that the models, tuning parameters, and implementation discussed here, despite our best efforts, may not be optimal. The method studied here could potentially be improved through careful tuning. Studying the proposed method further, including through extensive simulation experiments, could yield more concrete recommendations as to which embedded models perform best. Our implementation of a superlearner consisting of multiple estimated ITRs was shown to improve performance above individual ITRs. This approach helps mitigate user uncertainty regarding the best choice of ITR estimation approach or modeling choices for a given problem, where one may combine the results from multiple ITRs using the superlearner for improved performance.

Many of the modeling decisions in this paper were made with computational intensity in mind. In our analyses, we utilized two large computing clusters at the University of North Carolina at Chapel Hill, the Killdevil cluster with 9500 computing cores and the Longleaf cluster with 3600 computing cores. Future improvements on the proposed method must be implemented in a computationally efficient way. We also recognize that there are many methods we could have used that may have performed better, and that different methodological approaches entirely may be able to improve upon these results. Nevertheless, the analyses in this paper demonstrate the potential for using PDX studies to inform precision medicine.

The assumption that ITRs estimated from PDX data are applicable to humans is crucial to this work. This assumption is based on decades of biomedical research on generalizing PDX results to humans. A major conclusion of Gao et al. (2015) was that the responses observed in their PDX lines correlated with the responses observed in the human patients from which the tumors were taken. Many prior studies have shown that PDX models show stronger correlation with human response compared to traditional cell line models (Whittle et al., 2015; Scott et al., 2014; Rosfjord et al., 2014). Prior work has also shown that the tumor microenvironment, consisting of stromal cells and other tissue, may impact tumor activity and response to treatment. In PDX models, the microenvironment surrounding the implanted human tumor is non-human. However, several recent studies have suggested that tumor recruitment of mouse stroma in PDX models mirrors that seen in humans (Roife et al., 2016; Wang et al., 2017) and may show similarity in treatment response when specifically targeted.

A key next step for this research will be to plan and conduct validation studies in humans to determine if the biomarkers and ITRs discovered here can be used to improve outcomes in human patients. We note that ITRs based on only one genomic platform (RNA-seq) resulted in comparable performance to ITRs based on three genomic platforms. Since independent validation studies will be more expensive if more genomic platforms are needed, validating the ITRs based only on RNA-seq data using a study of human cancer patients would be a low-cost initial step toward validating the results in this paper. One advantage of the proposed method is that, in some settings, we could apply the estimated ITRs to external studies that included only a subset of the treatments applied in this study. Given a sequence of decision rules (pertaining to each step of the treatment tree), one could start at the lowest node that contains all of the treatments of interest and follow the tree to arrive at a recommended treatment, rather than starting from the top of the tree. We also note that, while our results allow for comparing the performance of estimated ITRs across different cancer types and modeling strategies, the data utilized for this paper do not allow for comparing the performance of ITRs estimated from PDX studies to ITRs estimated from human trials. Another important step forward for validating these results will be to compare the treatment rules discovered here to those estimated from human trials to determine the value of PDX studies for precision medicine.

The design of a PDX study plays a key role in ITR estimation, and designing high-quality PDX studies is another key next step for this research. Our results suggest that the observed responses in PDX studies are noisy. Having replicates within PDX line, i.e., multiple mice per line assigned to the same treatment, may improve the performance of the estimated ITRs. The design could also be improved by having a larger number of distinct tumor lines to ensure sufficient representation of cancer heterogeneity across a diverse spectrum of cancer patients.

While our research has, in some ways, raised more questions than it has answered, we feel that the treatment rules, biomarkers, and other results we have discovered here are interesting in their own right and merit further research, including validation studies. Future research in this area will allow us to make measurable advances in applying PDX studies in precision medicine and translating the results into clinical practice. Plans for such future research is underway.

References

  • S. Athey and G. Imbens (2016) Recursive partitioning for heterogeneous causal effects. Proceedings of the National Academy of Sciences 113 (27), pp. 7353–7360. Cited by: §1.
  • R. Bourgon, R. Gentleman, and W. Huber (2010) Independent filtering increases detection power for high-throughput experiments. Proceedings of the National Academy of Sciences 107 (21), pp. 9546–9551. Cited by: §1, §3.1.1.
  • J. Chen, H. Fu, X. He, M. R. Kosorok, and Y. Liu (2017) Estimating individualized treatment rules for ordinal treatments. arXiv:1702.04755 [stat]. Note: arXiv: 1702.04755 External Links: Link Cited by: §1, §3.2.4, §3.3.
  • Z. Chen, C. M. Fillmore, P. S. Hammerman, C. F. Kim, and K. Wong (2014) Non-small-cell lung cancers: A heterogeneous set of diseases. Nature Reviews Cancer 14 (8), pp. 535–546. Cited by: §1.
  • D. S. Ettinger, W. Akerley, G. Bepler, M. G. Blum, A. Chang, R. T. Cheney, L. R. Chirieac, T. A. D’Amico, T. L. Demmy, A. K. P. Ganti, et al. (2010) Non-small cell lung cancer. Journal of the National Comprehensive Cancer Network 8 (7), pp. 740–801. Cited by: §4.5.
  • J. Fan and J. Lv (2008) Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B 70 (5), pp. 849–911. Cited by: §3.1.2.
  • H. Gao, J. M. Korn, S. Ferretti, J. E. Monahan, Y. Wang, M. Singh, C. Zhang, C. Schnell, G. Yang, Y. Zhang, et al. (2015) High-throughput screening using patient-derived tumor xenografts to predict clinical trial drug response. Nature Medicine 21 (11), pp. 1318–1325. Cited by: §1, §1, §2.1, §2.2, §5.
  • M. Hidalgo, F. Amant, A. V. Biankin, E. Budinská, A. T. Byrne, C. Caldas, R. B. Clarke, S. de Jong, J. Jonkers, G. M. Mælandsmo, et al. (2014) Patient-derived xenograft models: An emerging platform for translational cancer research. Cancer Discovery 4 (9), pp. 998–1013. Cited by: §1.
  • K. Imai, M. Ratkovic, et al. (2013) Estimating treatment effect heterogeneity in randomized program evaluation. The Annals of Applied Statistics 7 (1), pp. 443–470. Cited by: §1.
  • N. Kumar (1981) Taxol-induced polymerization of purified tubulin. Mechanism of action.. Journal of Biological Chemistry 256 (20), pp. 10435–10441. Cited by: §4.5.
  • E. Laber and Y. Zhao (2015) Tree-based methods for individualized treatment regimes. Biometrika 102 (3), pp. 501–514. Cited by: §3.2.
  • Y. Liu, Y. Wang, M. R. Kosorok, Y. Zhao, and D. Zeng (2016) Robust hybrid learning for estimating personalized dynamic treatment regimens. arXiv:1611.02314 [stat]. Note: arXiv: 1611.02314 External Links: Link Cited by: §1, §3.2.4.
  • M. I. Love, W. Huber, and S. Anders (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15 (12), pp. 550. Cited by: §1, §3.1.1.
  • A. R. Luedtke and M. J. van der Laan (2016) Super-learning of an optimal dynamic treatment rule. The International Journal of Biostatistics 12 (1), pp. 305–332. Cited by: §1, §3.1, §3.3.1.
  • J. Ma, F. C. Stingo, and B. P. Hobbs (2016) Bayesian predictive modeling for genomic based personalized treatment selection. Biometrics 72 (2), pp. 575–583. Cited by: §1.
  • O. Metzger-Filho, A. Tutt, E. de Azambuja, K. S. Saini, G. Viale, S. Loi, I. Bradbury, J. M. Bliss, H. A. Azim Jr, P. Ellis, et al. (2012) Dissecting the heterogeneity of triple-negative breast cancer. Journal of Clinical Oncology 30 (15), pp. 1879–1887. Cited by: §1.
  • S. A. Murphy, M. J. van der Laan, J. M. Robins, and C. P. P. R. Group (2001) Marginal mean models for dynamic regimes. Journal of the American Statistical Association 96 (456), pp. 1410–1423. Cited by: §1.
  • S. A. Murphy (2005) A generalization error for Q-learning. Journal of Machine Learning Research 6 (Jul), pp. 1073–1097. Cited by: §1, §3.2.3.
  • L. Orellana, A. Rotnitzky, and J. M. Robins (2010) Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, part i: main content. The International Journal of Biostatistics 6 (2). Cited by: §1.
  • J. S. Parker, M. Mullins, M. C. Cheang, S. Leung, D. Voduc, T. Vickery, S. Davies, C. Fauron, X. He, Z. Hu, et al. (2009) Supervised risk predictor of breast cancer based on intrinsic subtypes. Journal of Clinical Oncology 27 (8), pp. 1160–1167. Cited by: §1.
  • R. Pirker, J. R. Pereira, A. Szczesna, J. Von Pawel, M. Krzakowski, R. Ramlau, I. Vynnychenko, K. Park, C. Yu, V. Ganul, et al. (2009) Cetuximab plus chemotherapy in patients with advanced non-small-cell lung cancer (FLEX): An open-label randomised phase III trial. The Lancet 373 (9674), pp. 1525–1531. Cited by: §4.5.
  • K. Polyak (2011) Heterogeneity in breast cancer. The Journal of Clinical Investigation 121 (10), pp. 3786–3788. Cited by: §1.
  • M. Qian and S. A. Murphy (2011) Performance guarantees for individualized treatment rules. The Annals of Statistics 39 (2), pp. 1180. Cited by: §1.
  • N. U. Rashid, A. S. Sperling, N. Bolli, D. C. Wedge, P. Van Loo, Y. Tai, M. A. Shammas, M. Fulciniti, M. K. Samur, P. G. Richardson, et al. (2014) Differential and limited expression of mutant alleles in multiple myeloma. Blood 124 (20), pp. 3110–3117. Cited by: §1, §3.1.1.
  • J. Robins, L. Orellana, and A. Rotnitzky (2008) Estimation and extrapolation of optimal treatment and testing strategies. Statistics in Medicine 27 (23), pp. 4678–4721. Cited by: §1.
  • D. Roife, B. Dai, Y. Kang, M. V. R. Perez, M. Pratt, X. Li, and J. B. Fleming (2016) Ex vivo testing of patient-derived xenografts mirrors the clinical outcome of patients with pancreatic ductal adenocarcinoma. Clinical Cancer Research 22 (24), pp. 6021–6030. Cited by: §5.
  • E. Rosfjord, J. Lucas, G. Li, and H. Gerber (2014) Advances in patient-derived tumor xenografts: From target identification to predicting clinical response rates in oncology. Biochemical Pharmacology 91 (2), pp. 135–143. Cited by: §5.
  • D.B. Rubin (1978) Bayesian inference for causal effects: the role of randomization. The Annals of Statistics 6 (1), pp. 34–58. Cited by: §3.2.2.
  • D. J. Sargent, B. A. Conley, C. Allegra, and L. Collette (2005) Clinical trial designs for predictive marker validation in cancer treatment trials. Journal of Clinical Oncology 23 (9), pp. 2020–2027. Cited by: §1.
  • P. J. Schulte, A. A. Tsiatis, E. B. Laber, and M. Davidian (2014) Q- and A-learning methods for estimating optimal dynamic treatment regimes. Statistical Science: A Review Journal of the Institute of Mathematical Statistics 29 (4), pp. 640–661 (eng). External Links: ISSN 0883-4237, Document Cited by: §1, §3.2.3.
  • C. L. Scott, H. J. Mackay, and P. Haluska Jr (2014) Patient-derived xenograft models in gynecological malignancies. In American Society of Clinical Oncology educational book/ASCO. American Society of Clinical Oncology. Meeting, pp. e258. Cited by: §5.
  • D. Siolas and G. J. Hannon (2013) Patient-derived tumor xenografts: Transforming clinical samples into mouse models. Cancer Research 73 (17), pp. 5315–5319. Cited by: §1.
  • C. Soneson and M. Delorenzi (2013) A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics 14 (1), pp. 91. Cited by: §1, §3.1.1.
  • G. J. Székely, M. L. Rizzo, et al. (2009) Brownian distance covariance. The Annals of Applied Statistics 3 (4), pp. 1236–1265. Cited by: §3.1.2.
  • J. J. Tentler, A. C. Tan, C. D. Weekes, A. Jimeno, S. Leong, T. M. Pitts, J. J. Arcaroli, W. A. Messersmith, and S. G. Eckhardt (2012) Patient-derived tumour xenografts as models for oncology drug development. Nature Reviews Clinical Oncology 9 (6), pp. 338–350. Cited by: §1.
  • P. Therasse, S. G. Arbuck, E. A. Eisenhauer, J. Wanders, R. S. Kaplan, L. Rubinstein, J. Verweij, M. Van Glabbeke, A. T. van Oosterom, M. C. Christian, et al. (2000) New guidelines to evaluate the response to treatment in solid tumors. Journal of the National Cancer Institute 92 (3), pp. 205–216. Cited by: §2.2.
  • R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight (2005) Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society: Series B 67 (1), pp. 91–108. Cited by: §3.2.1.
  • P. Vincent, H. Larochelle, I. Lajoie, Y. Bengio, and P. Manzagol (2010)

    Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion

    .
    Journal of Machine Learning Research 11 (Dec), pp. 3371–3408. External Links: ISSN ISSN 1533-7928 Cited by: §3.1.2.
  • S. Wager and S. Athey (2018) Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association 113 (523). Cited by: §1.
  • L. Wang and E. Laber (2017)

    Sufficient markov decision processes

    .
    Submitted. Cited by: §3.1.2.
  • X. Wang, A. D. Mooradian, P. Erdmann-Gilmore, Q. Zhang, R. Viner, S. R. Davies, K. Huang, R. Bomgarden, B. A. Van Tine, J. Shao, et al. (2017) Breast tumors educate the proteome of stromal tissue in an andividualized but coordinated manner. Sci. Signal. 10 (491), pp. eaam8065. Cited by: §5.
  • J. R. Whittle, M. T. Lewis, G. J. Lindeman, and J. E. Visvader (2015) Patient-derived xenograft models of breast cancer and their predictive power. Breast Cancer Research 17 (1), pp. 17. Cited by: §5.
  • D. O. Wise, R. Krahe, and B. R. Oakley (2000) The -tubulin gene family in humans. Genomics 67 (2), pp. 164–170. Cited by: §4.5.
  • T. Wu (2016) Set valued dynamic treatment regimes. Ph.D. Thesis, The University of Michigan. Cited by: §3.2.1.
  • B. Zhang, A. A. Tsiatis, M. Davidian, M. Zhang, and E. Laber (2012) Estimating optimal treatment regimes from a classification perspective. Stat 1 (1), pp. 103–114. Cited by: §1.
  • B. Zhang, A. A. Tsiatis, E. B. Laber, and M. Davidian (2013) Robust estimation of optimal dynamic treatment regimes for sequential treatment decisions. Biometrika 100 (3), pp. 681–694. Cited by: §1.
  • Y. Zhang, E. B. Laber, A. Tsiatis, and M. Davidian (2015) Using decision lists to construct interpretable and parsimonious treatment regimes. Biometrics 71 (4), pp. 895–904. Cited by: §3.2.
  • Y. Zhao, D. Zeng, A. J. Rush, and M. R. Kosorok (2012) Estimating individualized treatment rules using outcome weighted learning. Journal of the American Statistical Association 107 (499), pp. 1106–1118. Cited by: §1, §3.3.
  • Y. Zhao, M. R. Kosorok, and D. Zeng (2009) Reinforcement learning design for cancer clinical trials. Statistics in Medicine 28 (26), pp. 3294–3315. Cited by: §1, §3.2.3, §3.3.
  • Y. Zhao, D. Zeng, M. A. Socinski, and M. R. Kosorok (2011) Reinforcement learning strategies for clinical trials in nonsmall cell lung cancer. Biometrics 67 (4), pp. 1422–1433 (en). External Links: ISSN 1541-0420, Document Cited by: §3.2.3, §3.3.
  • X. Zhou, N. Mayer-Hamblett, U. Khan, and M. R. Kosorok (2017) Residual weighted learning for estimating individualized treatment rules. Journal of the American Statistical Association 112 (517), pp. 169–187. Cited by: §1, §3.2.4.