INTRODUCTION

A dynamic treatment regime is a set of sequential decision rules that provides a formal process for making treatment decisions at a series of key decision points throughout disease progression. At each decision point, the decision rule accepts patient history as input and returns a recommended treatment from among the feasible options. An optimal treatment regime leads to the most favorable outcome on average when it is used to select treatments for a patient population. Identifying an optimal regime is an integral component of precision medicine, an approach to healthcare that focuses on tailoring treatment decisions to patient characteristics.

For patients with potentially life-threatening conditions, identification of an optimal treatment regime that maximizes remaining life is of particular interest. Formally, the remaining life of a patient at time t, given the patient has survived up to time t, is referred to as residual life. Sepsis is a life-threatening medical condition that involves severe infections with organ dysfunction and is a leading cause of death worldwide (Singer et al., 2016). Although international guidelines for sepsis treatment have been established, treating septic patients remains highly challenging, as the heterogeneity of the septic patient population results in differing responses to medical intervention (László et al., 2015). Consequently, identifying optimal decision rules that maximize residual life and select treatment based on patients' individualized characteristics is highly desirable for septic patient care.

We are especially interested in identifying optimal treatment regimes for septic patients in the intensive care unit (ICU) using their electronic medical record (EMR) data. In particular, we study a data set of septic patients constructed from the Multiparameter Intelligent Monitoring Intensive Care database (MIMIC-III), a freely available database comprised of de-identified medical records for over 40,000 critical care patients (Johnson et al., 2016). Given a number of patient characteristics, including admission records, longitudinal vital signs, and longitudinal laboratory measurements, we seek to estimate an optimal treatment regime that maximizes residual life for the septic patients in MIMIC-III.

An overview of methods for estimating optimal treatment regimes based on censored time-toevent outcomes can be found in Chapter 8 of Tsiatis et al. (2020). One class of methods seeks to directly optimize an estimator for the value of a regime, where the value of a regime is defined in terms of the expectation of a given function of the event time, and optimization is conducted over all treatment regimes in a restricted class. Zhao et al. (2015), Bai et al. (2017), Cui et al. (2017), Zhou et al. (2021), and Wang et al. (2022) introduced direct optimization methods for regimes limited to only a single decision point, while Jiang et al. (2017a), Jiang et al. (2017b), Hager et al. (2018), Zhao et al. (2020), Xue et al. (2022), and Choi et al. (2023) presented methods for estimating optimal multistage regimes. In practice, implementing direct optimization methods with observational EMR data can be computationally prohibitive. First, EMR data often include a large number of covariates. Thus, the class of possible regimes derived from EMR data is likely to be large, rendering direct optimization computationally intensive. Second, estimators for the value of a regime often involve weighting by patients' propensities for being consistent with the given regime. In EMR data with a large number of possible regimes, the number of patients observed to be consistent with any given regime may be low, resulting in unstable value estimates. These complications become especially severe in settings with a large number of decision points.

Given the limitations of direct optimization, Q-learning has become a popular method for estimating optimal treatment regimes from observational data. Q-learning is a backward induction algorithm introduced in the reinforcement learning literature to solve multistage decision problems (Bellman, 1957;Watkins and Dayan, 1992). The Q-learning algorithm considers three key elements at each stage: a state, an action, and a reward. Q-learning estimates a sequence of decision rules, or a policy, that maximizes the expected sum of rewards. At each stage, a decision rule accepts the accumulated states and actions as input and returns a recommended action. Like direct optimization, Q-learning estimates become increasingly unstable as the number of stages increases. However, the Q-learning algorithm can more readily estimate optimal policies in settings with numerous covariates and low levels of observed consistency with studied policies.

In the clinical setting, a stage corresponds to a decision point, a state corresponds to patient history, an action corresponds to treatment, and a policy corresponds to a dynamic treatment regime. Moreover, a reward is defined to be some quantitative measure of patient health. The Qlearning algorithm estimates an optimal treatment regime that maximizes the expected sum of rewards using a backward recursive approach. First, the algorithm estimates an optimal decision rule at the final decision point, then the algorithm sequentially estimates an optimal decision rule at each previous decision point. At each decision point, a Q-model is posited for the expectation of the sum of current and future rewards, given patient history and treatment. To estimate optimal decision rules, each Q-model is optimized with respect to treatment, with all future treatment decisions taken to be optimal (Murphy, 2005).

The Q-learning framework has been extended to accommodate censored time-to-event outcomes. Zhao et al. (2011), Goldberg and Kosorok (2012), Huang et al. (2014), andZhao et al. (2020) all introduced Q-learning methods to estimate optimal treatment regimes that maximize expected survival time beyond the first decision point, where survival time is restricted to a constant. For regimes comprised of a fixed number of decision points, Zhao et al. (2011) presented a Qlearning method that used support vector regression to fit the Q-models. Goldberg and Kosorok (2012) extended the Q-learning algorithm to accommodate a flexible number of decision points, positing Q-models of an unspecified functional form adjusted by inverse probability of censoring weights. Huang et al. (2014) implemented a related approach that relied on fitting accelerated failure time models. Zhao et al. (2020) extended the Q-learning algorithm presented by Goldberg and Kosorok (2012) to allow shared decision rule parameters across decision points. Recent innovations in the field include an imputation-based Q-learning approach presented by Lyu et al. (2023), which uses non-or semiparametric models to estimate optimal treatment rules, followed by multiple imputation to predict optimal potential survival times. Additionally, Illenberger et al. (2023) introduced a combined Q-learning and policy-search method to estimate optimal list-based treatment regimes with a constraint imposed on expected treatment costs. Notably, all of the presented Q-learning methods defined the reward at each decision point to be a measure of the incremental amount of time between the given decision point and the next decision point or failure, whichever comes first.

Alternative reinforcement learning algorithms have been proposed for estimating optimal treatment regimes with censored time-to-event outcomes. One approach, called A-learning, employs a similar backward recursive strategy to Q-learning, but estimates the optimal decision rule at each decision point by optimizing only the portion of the regression outcome that involves differences between treatments (Murphy, 2003;Robins, 2004). While this makes A-learning robust to misspecification of the outcome models, A-learning also requires modeling the probability of the observed treatments given patient history at each decision point. Simoneau et al. (2020) and Zhang et al. (2022) introduced A-learning methods to estimate optimal treatment regimes that maximize expected restricted survival time beyond the first decision point. Alternatively, Cho et al. (2022) presented a backward recursive algorithm to estimate an optimal treatment regime that maximizes either mean survival time or survival probability at a given timepoint. As opposed to Qlearning and A-learning, the method presented by Cho et al. (2022) does not express the value of a regime in terms of a sum of rewards. At each decision point, Cho et al. (2022) appends conditional survival probability information, rather than a reward, in a backward recursive fashion.

To estimate optimal treatment regimes with censored time-to-event outcomes, we present a new Q-learning method based on an intuitive reward definition. While the aforementioned Qlearning methods define the reward to be the incremental time between decision points or failure, we argue that in clinical application, a more interpretable reward is a patient's remaining survival time, or residual life. Using this reward specification, the Q-learning algorithm estimates an optimal treatment regime that maximizes the expected sum of residual life, where the sum is taken across all decision points reached by the patient. We refer to this quantity as "expected cumulative residual life." To formulate the proposed Q-learning algorithm, we first introduce the residual life value estimator (ReLiVE), which estimates the expected cumulative residual life of a patient under a fixed treatment regime. We then build on ReLiVE to present the residual life value estimator Qlearning method (ReLiVE-Q). By defining the reward at each decision point to be residual life, ReLiVE-Q estimates an optimal treatment regime that maximizes expected cumulative residual life. To guarantee a finite expectation, we study restricted residual life, which we compute based on event times restricted to a fixed constant.

In addition to reward interpretability, ReLiVE-Q offers two advantages that make it especially suitable for estimating optimal treatment regimes from EMR data with censored time-to-event outcomes. First, the relationship between patient history, treatment, and event times is inherently complex and unlikely to be well-represented by a parametric model. At each decision point, ReLiVE-Q fits flexible, nonparametric Q-models using the random forest algorithm (Breiman,4 • G. Rhodes et al. 2001). Second, ReLiVE-Q incorporates longitudinal covariates into the Q-models using lowdimensional summaries of the trajectories, called context vectors, that are constructed using modern machine learning techniques (Rhodes et al., 2023). This is important because EMR data often contain a large number of longitudinal covariates measured over a potentially dense time grid, and the number and timing of measurements often differs among patients. Thus, including the entire patient history directly in the Q-models is often computationally prohibitive. Moreover, while summary statistics of the longitudinal covariates could be included in the Q-models, simple summaries are unlikely to provide an adequate synthesis of the complex information contained in the trajectories.

In Section 2, we present the statistical framework for ReLiVE and ReLiVE-Q, and in Section 3, we detail the ReLiVE and ReLiVE-Q methodology. In Section 4, we demonstrate the utility of ReLiVE-Q in a simulation study, and in Section 5, we apply ReLiVE-Q to estimate optimal treatment regimes for septic patients in MIMIC-III. In both the simulation study and MIMIC-III data application, we demonstrate that ReLiVE-Q estimates personalized treatment regimes that optimize a clinically meaningful function of residual life. We conclude with a discussion of implications and open problems in Section 6.

STATISTICAL FR A ME WOR K 2.1. Potential outcomes and treatment regimes

We focus our study on settings where treatment decisions are made at fixed intervals, as is the scenario for the studied MIMIC-III data set. Consider a series of K decision points at which treatment decisions must be made, where the decision points occur at fixed times τ 1 , . . . , τ K with τ 1 = 0. Let x 1 denote the vector of covariates measured at τ 1 . Then the patient history at decision point 1 is h 1 = (τ 1 , x 1 ). For k = 1, . . . , K, define A k to be the finite set of all available treatment options at decision point k, and let a k ∈ A k be the treatment administered at decision point k. We assume all a k ∈ A k are feasible for all subjects at all decision points k = 1, . . . , K, as this assumption holds for the treatments studied in MIMIC-III. For k = 2, . . . , K, let x k denote the vector of covariate information accumulated between decision points k -1 and k. Moreover, let τk = (τ 1 , . . . , τ k ), xk = (x 1 , . . . , x k ), and āk = (a 1 , . . . , a k ), where k = 1, . . . , K, and let ā = āK . Then for k = 2, . . . , K, a patient who does not experience the event of interest prior to decision point k will have patient history h k = ( τk , xk , āk-1 ) at decision point k. In contrast, a patient who experiences the event of interest at time t between decision points k -1 and k will have patient history

For each decision point, k = 1, . . . , K, let d k (h k ) be a decision rule. If h k indicates that the event of interest has not occurred prior to decision point k, then d k maps h k to A k . Else, if h k indicates that the event of interest has occurred prior to decision point k, then d k (h k ) returns null. For notational simplicity, suppress the dependence of the decision rule on patient history, and let To formalize the definition of the residual life that would occur if an individual were to follow a given regime d ∈ D, we define the potential outcomes that would be achieved if a randomly selected patient were treated according to d. To this end, we first define the potential outcomes associated with a given sequence of treatments ā. Let κ * (ā) be the potential number of decision points a randomly selected patient would reach if administered treatments ā, where 1

) be the vector of potential covariate information that would occur between decision points k -1 and k if a patient was administered treatments āk-1 . Further, let T * (ā κ * (ā) ) be the potential time-to-event that would be observed if a patient was administered treatments āκ * (ā) . Define the set of all possible potential outcomes to be

As discussed in Sections 6.2.3 and 8.3.2 of Tsiatis et al. (2020), it can be shown that the potential outcomes associated with a given regime d ∈ D can be defined in terms of W * . Let T * (d) be the Biostatistics, 2024, 00, 0 • 5 potential time-to-event under regime d. To avoid infinite event times, also specify a restricted event time

, L}] is the potential number of decision points reached under regime d. For k = 2, . . . , κ * L (d), let X * k ( dk-1 ) be the vector of potential covariate information that would occur between decision points k -1 and k if a randomly selected patient was treated according to regime dk-1 . Denote the set of potential outcomes associated with regime

Because sepsis is a life-threatening medical condition, we specify T * (d) to be the potential time of death under regime d. Then, restricted residual life under regime d at a given time τ is min{T * (d), L}τ . We specify the reward at each decision point k = 1, . . . , κ * L (d) to be the patient's restricted residual life under regime d at τ k . Let the cumulative restricted residual life of a patient under regime d be the sum of the patient's restricted residual life under d, taken across all decision points reached by the patient. We define the value of regime d ∈ D, V(d), to be the expected value of cumulative restricted residual life under d. That is, we define

An optimal treatment regime d opt ∈ D satisfies the condition that

Thus, an optimal regime d opt maximizes expected cumulative restricted residual life across regimes d ∈ D. Based on observed EMR data, we first aim to estimate V(d) for any fixed regime d ∈ D. We then aim to estimate an optimal regime d opt and its value V(d opt ) from the observed data.

Data and assumptions

We now describe the observed data. As is conventional in the survival analysis literature, let T > 0 and C > 0 denote the potential times to death and censoring, respectively. We observe only U = min(T, C) and = I(T < C), an indicator of whether death is observed ( = 1) or censored ( = 0). Accounting for the restricted lifetime L > 0, we observe the restricted outcome U L = min(U, L) and associated indicator L = I{min(T, L)

Then the observed history at decision points k = 1, . . . , κ L is H k = ( τk , Xk , Āk-1 ), and the observed history at decision points

We aim to estimate value (2.1) for a given regime d ∈ D, and to estimate an optimal regime d opt and its value, based on the observed data. Because T * (d) is defined in terms of W * , we must be able to express relevant functionals of the distribution of the potential outcomes in terms of the observed data. Such expression is possible under three key, standard assumptions: the Stable Unit Treatment Value Assumption (SUTVA), the Sequential Randomization Assumption (SRA), and the positivity assumption. SUTVA, commonly referred to as the consistency assumption, states that the observed data are the same as those that would potentially be achieved under the treatment decisions observed to be administered. Formally, we assume κ L = κ * L ( Āκ L ) and T = T * ( Āκ L ). We also assume X k = X * k ( Āk-1 ) for k = 2, . . . , κ L . SRA states that the treatment administered at each decision point is independent of the set of potential outcomes conditional on the history. For k = 1, . . . , K, we assume W * ⊥ A k |H k , κ L ≥ k, where ⊥ denotes independence. The positivity assumption states that all treatment options at each decision point are represented in the data. This is necessary to specify the distributions of the potential outcomes for all d ∈ D in terms of the observed data. Formally, we assume

For simplicity, we also assume censoring is independent of treatment assignment, patient characteristics, restricted lifetime, and potential outcomes. See Section 6 for a discussion on relaxing this assumption.

METHODS

Methodological foundation of ReLiVE

We first introduce the methodological foundation for ReLiVE, our proposed estimator for value (2.1) of a fixed regime. For a given regime d ∈ D, we construct an estimator for V(d) by defining K backward recursive Q-functions. We incorporate inverse-probability weights in the Q-functions to account for censoring. For

where k = 1, . . . , K. Under SUTVA and the independent censoring assumption, K k (u|h k , a k ) = P(C > u)/P(C > τ k ). See Section 1 of the supplementary material for details.

For a fixed regime d ∈ D, we define the Q-function at decision point K to be

Similarly, we define the Q-function at decision points k = K -1, . . . , 1, to be

Under SUTVA, SRA, and the positivity and independent censoring assumptions, V d 1 (h 1 ) is an unbiased estimator of V(d) if the true Q-functions are known. In Section 2 of the supplementary material, we provide a proof demonstrating

In practice, the Q-functions for regime d are unknown and must be estimated from the data. In Section 3.5, we present the value estimator ReLiVE, which uses a flexible approach to model and estimate the Q-functions.

Methodological foundation of ReLiVE-Q

Next, we introduce the methodological foundation for ReLiVE-Q, our proposed Q-learning method for estimating an optimal treatment regime and its value. Building on the Q-functions presented in Section 3.1, we characterize an optimal treatment regime d opt ∈ D such that V(d opt ) ≥ V(d) for all d ∈ D. At decision point K, define the Q-function for d opt to be

and let

SRA, and the positivity and independent censoring assumptions, it can be shown that

) is an optimal treatment regime satisfying V(d opt ) ≥ V(d) for all d ∈ D, and the value of an optimal regime is V(d opt ) = E{V d opt 1 (h 1 )}. See Sections 7.2.3-7.2.4 of Tsiatis et al. (2020) for details. In practice, the Q-functions for d opt are unknown and must be estimated from the data. In Section 3.6, we present ReLiVE-Q, a Q-learning method that uses a flexible approach to estimate the Q-functions for d opt .

Downloaded from https://academic.oup.com/biostatistics/advance-article/doi/10.1093/biostatistics/kxae002/7604037 by North Carolina State University Libraries user on 18 July 2024 Biostatistics, 2024, 00, 0 • 7

Representation of longitudinal covariates

In practice, certain covariates, such as sex and birth date, are collected on patients only at baseline. Denote the vector of covariates measured only at baseline as S. Other covariates, such as vital signs and laboratory values, may be measured repeatedly during a patient's follow up. Let M be the set of patient-specific measurement times, and denote the vector of longitudinal covariates measured at time t ∈ M as Z(t) = {Z 1 (t), . . . , Z p (t)}. For longitudinal covariates measured at time t = τ k , k = 1, . . . , κ L , we follow the convention that longitudinal covariates are measured immediately prior to making treatment decisions. Then X 1 = {S, Z(0)}, and

When estimating the Q-functions for ReLiVE and ReLiVE-Q, it is desirable to leverage all available patient information. Recall that EMR data often contain a large number of longitudinal covariates potentially measured over a dense time grid, and the number and timing of measurements often differs among patients. To incorporate the trajectories of the longitudinal covariates into the Q-models specified in Sections 3.5 and 3.6, it is expedient to construct lower-dimensional summary representations of the longitudinal trajectories. Let f (•) be a vector-valued function, and define vector

to be a function of the longitudinal covariate measurements accumulated by decision point k. Then, X1 = {S, Z(0)}, and for k = 2, . . . , κ L , we define Xk = {S, ζ (k)}. The observed patient history is then given by H 1 = (τ 1 , X 1 ) = {τ 1 , S, Z(0)} at decision point 1, and H k = ( τk , Xk , Āk-1 ) = { τk , S, ζ (k), Āk-1 } at decision points k = 2, . . . , κ L . At decision points k = 1, . . . , K, we train the Q-models for ReLiVE and ReLiVE-Q over the predictor space Xk , as described in Section 3.4.

Intuitively, it is desirable to select a function f (•) that summarizes the longitudinal trajectories from baseline to decision point time. Several simple summary functions are commonly used in practice. The baseline vector ζ B (k) contains the longitudinal measurements recorded at baseline, the last-value carried forward vector ζ L (k) contains the most recently observed measurements, and the average vector ζ A (k) contains the average of the measurements observed through decision point k. Because these simple summary functions are unlikely to capture the complex nature of the longitudinal trajectories, we propose synthesizing the longitudinal covariates using window-specific context vectors constructed via long short-term memory (LSTM) autoencoders, as described in Rhodes et al. (2023). The window-specific context vector ψ l (τ ) is an encoded representation of the trajectory of longitudinal covariate Z l (•) from baseline to time τ > 0, where l = 1, . . . , p. Specifically, we propose defining ζ C (k) = {ψ 1 (τ k ), . . . , ψ p (τ k )}. A detailed description of the studied forms of ζ (•) can be found in Section 3 of the supplementary material.

Q-function estimation strategy

We introduce methods to estimate the Q-functions for ReLiVE and ReLiVE-Q in Sections 3.5 and 3.6, respectively. Although extension to more general treatment settings is possible, we focus on settings with two treatment options such that A k = {0, 1} for k = 1, . . . , K. We model the Qfunctions using a non-parametric approach to account for the inherently complex relationship between patient history, treatment, and residual life. Specifically, we estimate

) at each decision point k = 1, . . . , K using distinct random forests, where the outcomes used to train the forests are presented in Sections 3.5 and 3.6. A random forest is an ensemble learning algorithm that conducts estimation by combining the output of numerous decision trees. For a detailed overview of the random forest algorithm, see Breiman (2001). At decision point 1, we train the random forests over the predictor space X1 = {S, Z(0)}. Given a function ζ (•), at decision points k = 2, . . . , K, we train the random forests over the predictor space Xk = {S, ζ (k)}. We train the random forest for Q d k (h k , 0) on patients having U L > τ k and A k = 0, and we train the random forest for Q d k (h k , 1) on patients having U L > τ k and A k = 1.

ReLiVE: value estimation for a fixed regime

We now introduce ReLiVE, which estimates value (2.1) for a fixed regime

is the Kaplan-Meier estimator for the survival function of censoring time C at time t (Kaplan and Meier, 1958). Then Kk (u|h k , a k ) = Kk (u). Beginning with decision point K, given patient history H K = h K , we train the random forests for

Given m patients, we define ReLiVE as

where Vd

3.6. ReLiVE-Q: estimation of an optimal regime and its value

Next, we introduce ReLiVE-Q, which estimates an optimal treatment regime d opt ∈ D and its value V(d opt ). Beginning with decision point K, given patient history H K = h K , we train the random forests for

where

, where a k = 0 is the standard treatment. We estimate an optimal treatment regime to be dopt = ( dopt 1 , . . . , dopt K ). Given m patients, we estimate the value of an optimal treatment regime to be

Testing and validation procedures

To evaluate the estimation performance of ReLiVE-Q, we implement testing and validation procedures that compare value estimates for an optimal treatment regime to value estimates for two fixed treatment regimes, the observed treatment regime d obs and the no treatment regime d no . At each decision point k = 1, . . . , K, d obs specifies d k (h k ) = a k , where a k is the treatment observed in the data at decision point k, and d no specifies d k (h k ) = 0. Thus, d obs represents the treatment regime that was actually administered by clinicians in practice, and d no represents the treatment regime that always administers the standard treatment. In the testing and validation procedures, we compare N value estimates for each regime d ∈ {d opt , d obs , d no }. A detailed description of the testing and validation procedures is provided in Section 4 of the supplementary material.

In summary, the testing procedure obtains N value estimates Vtest (d) for each regime d ∈ {d opt , d obs , d no } using N unique training/testing data sets. The procedure computes Vtest (d) for the Biostatistics, 2024, 00, 0 • 9 fixed regimes d ∈ {d obs , d no } by conducting ReLiVE separately on the training and testing data sets. The testing procedure also estimates dopt and its value Vtest (d opt ) via a cross-validated adaptation of ReLiVE-Q. Moreover, the validation procedure obtains N value estimates Vval (d) for each regime d ∈ {d opt , d obs , d no } by exploiting the known data generation process in a simulation study. Given a simulated data set of patient covariates and a regime d, the validation procedure randomly generates N potential survival times T * (d) for each patient. Using each of the N randomly generated sets of T * (d), the procedure computes Vval (d) d), L}τ j ] across patients.

We implement the testing and validation procedures using simulated data in Section 4, and we implement the testing procedure using the MIMIC-III data set in Section 5. In both sections, we repeat the testing and validation procedures four times, defining ζ (•) in terms of the baseline vector ζ B (•), the last-value carried forward vector ζ L (•), the average vector ζ A (•), and the context vectors ζ C (•), as described in Section 3.3. Note, Vtest (d opt ), Vtest (d obs ), and Vtest (d no ) are dependent on the fitted Q-models, as is Vval (d opt ). Thus, these value estimates are dependent on the functional form of ζ (•). In contrast, Vval (d obs ) and Vval (d no ) are independent of the fitted Q-models, so these estimates do not depend on ζ (•). We evaluate the estimation performance of ReLiVE-Q by studying the results of the testing and validation procedures in Sections 4 and 5.

SIMUL ATION STUDY 4.1. Simulation strategy

We design a simulation study to evaluate the estimation performance of ReLiVE-Q. We provide a detailed description of the simulation strategy in Section 5 of the supplementary material. In summary, we conduct the testing and validation procedures described in Section 3.7 to obtain N = 500 value estimates for each studied regime using simulated data. We consider K = 4 decision points that occur at times (τ 1 , . . . , τ 4 ) = (0, 3, 6, 9), and we generate a data set of m = 10, 000 patients. For each patient, we generate a single treatment variable A k ∈ A k = {0, 1} at k = 1, . . . , 4, and we generate a single covariate measured only at baseline, S. We consider ten measurement times M = (0, 1, . . . , 9), and we generate three longitudinal covariates Z l (t) = B l (t) + l (t) at each t ∈ M, where l = 1, 2, 3. We generate censoring times C and impose a restricted lifetime of L = 50. We conduct two separate analyses, each using a distinct survival time generation process. First, we generate T according to an accelerated failure time (AFT) model. Second, we generate T according to a Cox proportional hazards model. We provide the code used to conduct the simulation study in the supplementary material.

Simulation results

The 500 testing value estimates Vtest (d) and validation value estimates Vval (d) are plotted in Figures 1 and 2 for d ∈ {d opt , d obs , d no }. In both the AFT and Cox analyses, the distributions of Vtest (d) and Vval (d) are consistently higher for the optimal treatment regime than for the observed treatment regime or the no treatment regime. Thus, the simulation study supports that V( dopt ) > V(d obs ) and V( dopt ) > V(d no ), which allows us to conclude that ReLiVE-Q successfully produces reasonable estimates of an optimal treatment regime. Moreover, in both the AFT and Cox analyses, Q-models fit with ζ B (•) result in lower value estimates for the optimal treatment regime than those fit with ζ A (•), ζ L (•), or ζ C (•). In the AFT analysis, Q-models fit with ζ L (•) or ζ C (•) result in the highest value estimates for the optimal treatment regime, while in the Cox analysis, Q-models fit with ζ C (•) result in the highest. Thus, the simulation study supports that sophisticated functions are necessary to synthesize the trajectories of the longitudinal covariates, and that representing the longitudinal covariates with context vectors can lead to improved optimal treatment regime estimation via ReLiVE-Q.

Figure 1. AFT simulation study: boxplots of the value estimates from the validation procedure (left) and the testing procedure (right) for an optimal treatment regime (Opt), the observed treatment regime (Obs), and the no treatment regime (No). For scenarios dependent on the Q-models, value estimates are presented using the baseline vector (B), the average vector (A), the last-value carried forward vector (L), and the context vector (C).

Figure 2. Cox simulation study: boxplots of the value estimates from the validation procedure (left) and the testing procedure (right) for an optimal treatment regime (Opt), the observed treatment regime (Obs), and the no treatment regime (No). For scenarios dependent on the Q-models, value estimates are presented using the baseline vector (B), the average vector (A), the last-value carried forward vector (L), and the context vector (C).

MIMIC-III DATA APPLICATION 5.1. MIMIC-III data description

Next, we evaluate the estimation performance of ReLiVE-Q by conducting the testing procedure described in Section 3.7 on the EMR data of septic patients in the MIMIC-III database. MIMIC-III is a freely available database comprised of deidentified medical records for over 40,000 patients who stayed in the critical care units at Beth Israel Deaconess Medical Center between 2001 and 2012 (Johnson et al., 2016). MIMIC-III contains data on patients' demographics, vital signs, laboratory measurements, medications, imaging reports, chart notes, procedure codes, diagnostic codes, hospital stay, and survival. For a complete description of the MIMIC-III database, refer to Johnson et al. (2016).

In 2016, the definitions and clinical criteria for sepsis and septic shock were updated in the Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3) (Singer et al., 2016). Sepsis-3 defines sepsis as a "life-threatening organ dysfunction caused by a dysregulated host response to infection" and provides clinical criteria for diagnosing septic patients (Singer et al., 2016). Komorowski (2019) developed code to identify patients in MIMIC-III fulfilling the Sepsis-3 criteria. Komorowski's code pulls relevant physiological parameters for septic patients from up to 24 h preceding their sepsis diagnosis until 48 h after. The code aggregates the data into 4-h time windows, recording an appropriate summary statistic when several measurements were taken in the same time window. We use Komorowski's code to construct our studied data set of m = 20, 952 patients, only 15% of whom have observed restricted survival times, i.e. L = 1.

We estimate optimal treatment regimes for septic patients in the ICU using the MIMIC-III data set. We conduct the testing procedure described in Section 3.7 to obtain N = 100 value estimates for the studied regimes. Technical details are provided in Section 7 of the supplementary material. We consider K = 10 decision points that occur every 4 h at (τ 1 , τ 2 , . . . , τ 10 ) = (0, 4, . . . , 36). We study a single treatment A k ∈ {0, 1} at each decision point, k = 1, . . . , 10, where A k = 1 if the patient is provided vasopressors in the 4-h time window [τ k , τ k + 4) and A k = 0 otherwise. We specify a restricted lifetime of L = 40, and we study a predictor space containing a single baseline covariate S and 13 longitudinal covariates Z(•). The baseline covariate of interest is an indicator of whether the patient was previously admitted to the ICU during the given hospital stay. The longitudinal covariates of interest include a longitudinal indicator of mechanical ventilator dependence, as well as 12 longitudinal vital signs and laboratory values: albumin, calcium, magnesium, hemoglobin, arterial lactate, arterial pH, fraction of inspired oxygen (FiO 2 ), peripheral oxygen saturation (SpO 2 ), Sequential Organ Failure Assessment (SOFA) score, respiratory rate, heart rate, and systolic blood pressure.

MIMIC-III results

The 100 testing value estimates Vtest (d) are plotted in Fig. 3 for d ∈ {d opt , d obs , d no }. The distribution of Vtest (d) is consistently higher for the optimal treatment regime than for the observed treatment regime or the no treatment regime. These results support that V( dopt ) > V(d obs ) and V( dopt ) > V(d no ). Thus, we can conclude that ReLiVE-Q successfully uses EMR data to produce reasonable optimal treatment regime estimates for septic patients in the ICU.

Again, Q-models fit with ζ B (•) result in lower value estimates for the optimal treatment regime than those fit with

result in the highest value estimates for the optimal treatment regime. This application suggests that representing longitudinal covariates with context vectors can lead to improved optimal treatment regime estimation via ReLiVE-Q. In Section 8 of the supplementary material, we evaluate the importance of each studied covariate for accurately estimating the Q-functions.

DISCUSSION

In the simulation study and application to MIMIC-III, we demonstrate that the optimal treatment regime estimated via ReLiVE-Q results in higher estimates of value (2.1) than the no treatment and observed treatment regimes. Thus, we can expect cumulative restricted residual life to be higher on average for patients who follow the estimated optimal treatment regime, as compared to patients who always receive the standard treatment, and compared to patients who received the observed Figure 3. MIMIC-III application: boxplots of the value estimates from the testing procedure for an optimal treatment regime (Opt), the observed treatment regime (Obs), and the no treatment regime (No). Value estimates are presented from Q-models using the baseline vector (B), the average vector (A), the last-value carried forward vector (L), and the context vector (C). treatments administered by clinicians. Thus, ReLiVE-Q leverages patient history to estimate personalized treatment regimes that maximize a clinically meaningful function of residual life. This finding is especially important for patients diagnosed with sepsis, as septic patients experience varying responses to treatment for the life-threatening condition. Moreover, we demonstrate that synthesizing longitudinal covariates with context vectors in ReLiVE-Q leads to improved optimal treatment regime estimation, as compared to using simpler summary statistics.

We limit our study to clinical settings with two treatment options that are feasible for all patients at all decision points, where decision points occur at fixed intervals. Extending our methods to accommodate a finite number of treatment options is straightforward. Our methods can be extended to settings where not all treatments are feasible for all patients under the feasible sets framework outlined in Tsiatis et al. (2020). By re-framing the decision point times as potential outcomes dependent on patient history, our methods can also be extended to accommodate subjectspecific decision point times. See Sections 8.3.1-8.3.2 of Tsiatis et al. (2020) for an example of this approach. Further research is required to determine how ReLiVE-Q performs in these more complex clinical settings.

We assume censoring is independent of treatment assignment, patient characteristics, restricted lifetime, and potential outcomes, and we estimate the censoring weights via the Kaplan-Meier method. However, it may be possible to relax this assumption and take censoring to be noninformative in the sense that the cause-specific hazard of censoring as a function of patient history and potential outcomes does not depend on the potential outcomes. In this case, a time-dependent Cox model could be used to incorporate covariate information into the censoring weight estimates. See Section 8.3.2 of Tsiatis et al. (2020) for details.

In this article, we focus on accurately estimating optimal treatment regimes. Accordingly, we model the Q-functions using flexible, nonparametric random forests. A drawback to modeling the Biostatistics, 2024, 00, 0 • 13 Q-functions using random forests is that the estimated optimal treatment rules are not readily interpretable. ReLiVE-Q could be implemented with simpler Q-models to improve decision rule interpretability, though further research is required to determine how this would affect estimation.

ReLiVE-Q defines the reward at each decision point to be residual life. This quantity is more interpretable than the rewards used in the alternative Q-learning methods described in Section 1, defined as the incremental amount of time between the given decision point and the next decision point or failure. However, ReLiVE-Q's improvement in reward interpretability comes with the trade-off of a less interpretable value function. ReLiVE-Q defines the value of regime d ∈ D to be the expected value of cumulative restricted residual life under d summed across all decision points reached by the patient, while the Q-learning methods discussed in Section 1 more simply define the value of regime d to be the expected value of restricted residual life under d.