Dynamic treatment regimes are a set of time-adaptive decision rules that can be used to personalize treatment across multiple stages of care. Grounded in causal inference methods, dynamic treatment regimes identify variables that differentiate the treatment effect and may be used to tailor treatments across individuals based on the patient’s own characteristics – thereby representing an important step toward personalized medicine. In this manuscript we introduce Penalized Spline-Involved Tree-based Learning, which seeks to improve upon existing tree-based approaches to estimating an optimal dynamic treatment regime. Instead of using weights determined from the estimated propensity scores, which may result in unstable estimates when weights are highly variable, we predict missing counterfactual outcomes using regression models that incorporate a penalized spline of the propensity score and other covariates predictive of the outcome. We further develop a novel purity measure applied within a decision tree framework to produce a flexible yet interpretable method for estimating an optimal multi-stage multi-treatment dynamic treatment regime. In simulation experiments we demonstrate good performance of Penalized Spline-Involved Tree-based Learning relative to competing methods and, in particular, we show that Penalized Spline-Involved Tree-based Learning may be advantageous when the sample size is small and/or when the level of confounding of the outcome is high. We apply Penalized Spline-Involved Tree-based Learning to the retrospectively-collected Medical Information Mart for Intensive Care dataset to identify variables that may be used to tailor early fluid resuscitation strategies in septic patients.
Personalized medicine is built upon the understanding that patients are uniquely heterogeneous in their existing and emergent comorbidities, as well as their tolerance of, response to, and even preference for different treatments. Given the increasing prevalence of chronic health conditions, as well as the rapid increase in healthcare expenditures overall, large scale initiatives to deliver personalized medicine are underway. One such avenue to advance personalized medicine is through dynamic treatment regimes (DTRs),1,2 the statistical methods of which are grounded in causal inference. DTRs, also known as adaptive interventions, are a series of stage-specific decision rules that map a patient’s measured baseline and time-varying characteristics to a treatment assignment at each successive stage. One particular objective within the field is to estimate an optimal DTR such that, if the population of interest were to receive treatment consistent with this regime, overall patient-level outcomes would be optimized.
DTR estimation methods are often classified based on their degree of dependence on parametric assumptions. With the abundance of observational data available to us, flexible estimation methods, which are able to account for what is expected to be a complex relationship among variables of interest, are often desired. Additionally, because optimal DTR estimation is largely an exploratory process and collaboration with clinician-scientists is critical, the need for interpretability in an estimated optimal DTR is paramount. As a result, flexible and robust methods that yield interpretable results—for example, those with a decision tree-type structure—have been enjoying much popularity. Over the past decade tree-based methods have evolved from the ability to handle a single stage and/or binary treatment setting3–5 to a multi-stage setting with multiple treatment options per stage,6–8 which better reflects how care for chronic health conditions is delivered in practice. Zhang et al.8 estimate an optimal multi-stage DTR using a decision list; however, computational demands restrict each statement to a maximum of two covariates and, additionally, the unidirectional growth of decision lists precludes correction of estimation error(s) that may have occurred at previous steps. Sun and Wang6 propose a stochastic decision tree search process that estimates counterfactual outcomes using bayesian additive regression trees, which can be quite complicated and computationally intensive. Tao et al.7 develop tree-based reinforcement learning (T-RL), which cleverly embeds into the decision tree framework a purity measure devised from the augmented inverse probability weighted (AIPW) estimate of the counterfactual mean outcome. Although the AIPW estimator is consistent and doubly robust for the counterfactual mean outcome, it has been well established that IPW-style estimators are unstable when weights are highly variable, which will often be the case with low propensity of treatment assignment and/or as the number of stages increases.
In this article, we propose Penalized Spline-Involved Tree-based (PenSIT) Learning, which seeks to improve upon existing tree-based approaches for estimating an optimal multi-stage multi-treatment DTR. While conceptually similar to the implementation of T-RL, PenSIT Learning makes use of a different purity measure—one that uses a Penalized Spline-Involved (PenSI) estimator of counterfactual outcomes developed from the penalized spline of propensity prediction method used for missing data.9–11 Specifically, we predict missing counterfactual outcomes for the treatments not assigned to patients using regression models that incorporate a penalized spline function of the propensity assigning that treatment and other covariates predictive of the outcome. The PenSI estimator of the counterfactual mean outcome, like the AIPW estimator, is consistent and retains the property of double robustness against model misspecifications, which may lend more stability and provide improved performance (e.g., a higher percentage of observations correctly classified to their optimal multi-stage DTR) under certain data generating mechanisms. PenSIT Learning estimates stage-specific optimal decision rules using backward induction, beginning with the final stage, to remove bias arising from confounding by indication in the multi-stage treatment situation. PenSIT Learning is a viable alternative to T-RL for estimation of an optimal multi-stage DTR and may be advantageous in correctly identifying the optimal multi-stage treatment sequence when the underlying DTR is tree-based, sample sizes are small, and/or when the level of confounding is high. Additionally, due to the modeling flexibility afforded by PenSIT Learning, it may be preferred to T-RL in many situations.
In this article, we introduce relevant notation and formulation, followed by PenSIT Learning methodology. We then present simulation results and, to follow, describe the application of PenSIT Learning to data obtained from the Medical Information Mart for Intensive Care (MIMIC-III) Clinical database to estimate an optimal, two-stage DTR reflecting restrictive or liberal fluid resuscitation strategies designed to minimize a measure of multi-organ dysfunction.
Notation and formulation
Notation
Suppose we are estimating a -stage treatment regime () in which one of treatments (; is administered to every subject . Treatment received by the th individual at the th stage is denoted , with assumed to be categorical. As is customary, a capital letter denotes the random variable with a lower case letter denoting the realized value. We suppress the patient-level indicator when it can be safely omitted without confusion. Variables collected and available when making the th treatment decision are denoted . Following the th stage treatment , measurements are collected on a set of covariates , which may also include an intermediate reward outcome, . We denote the full covariate and treatment history prior to the decision at stage as . A final outcome of interest is a clinically relevant, prespecified function of observed, stage-specific intermediate reward outcomes, the higher the better by convention. Common functions for , for example, include the sum or last value. Our fully observed data, then, represent the collection of independent and identically distributed multivariate observations from subjects in our population of interest and is summarized as follows: . Now, further to our estimation goal, we let denote a -stage DTR. Each stage-specific decision rule maps history to the -th treatment decision, that is, . Thus, we can more specifically express as .
Following Rubin’s potential outcomes framework,12 we use , or simply , to denote the counterfactual outcome, known interchangeably as a “potential outcome,” for a patient treated with conditional on prior treatment history, . Similarly, identifies the counterfactual outcome under a treatment sequence and denotes the counterfactual outcome under regime . Using the mean counterfactual outcome, , to evaluate performance of a DTR , the optimal DTR, , is the one that satisfies
for all , where is the class of all potential regimes. Our statistical goal, therefore, is to estimate an interpretable, optimal, -stage DTR, , using observational data such that, if all patients were to be assigned to multi-stage treatment using this regime, the expected counterfactual outcome of our population of interest would be maximized: .
Link to observed data
As mentioned previously, only one of the counterfactual outcomes is observed, making estimation of impossible without a series of assumptions. Therefore, we make the following three foundational assumptions: consistency, positivity, and ignorability.13
Consistency: The potential outcome under the observed treatment agrees with that of the observed outcome. For example, in the final stage , we can express this as: , where is an indicator function that returns a value of if the argument is true and a value of otherwise. Consistency further assumes that there is no interference between units, which means that one patient’s observed and counterfactual outcomes are independent of the treatments of all other patients.
Positivity: An assumption of positivity is fulfilled if there is a positive probability for each subject of being assigned at the th stage treatment decision conditional on history . Expressed mathematically, under the positivity assumption for all and for all , where represents the propensity score for subject and represents a positive constant.
Ignorability: Also known as the assumption of no unmeasured confounders (NUCA), ignorability implies that data on all variables that are associated both with the assignment of and the outcome have been observed. Furthermore, ignorability implies that the counterfactual outcomes are independent of the treatment given the propensity score. For example, .10,11,14 Under the assumptions of consistency, positivity, and ignorability, it can be shown that . Following the derivation for a single stage in Tao and Wang,15 we can express the optimal decision rule for the final stage as:
Further, following the theory of Rosenbaum and Rubin,14 conditioning on all the full history, , is mathematically the same as conditioning on the corresponding propensity score. Therefore, we also have:
Note that equations (1) and (2) are solvable using observed data. Based on the above, we formally introduce our proposed PenSIT Learning in the following section.
As is well known in the classification and regression tree (CART) literature16 and also as discussed in tree-based optimal DTR estimation literature,3,7 a decision tree is a widely-used machine learning technique constructed by recursive partitioning of the covariate space that is used to identify features and their associated cutpoints that are best able to describe the relative homogeneity of another variable or outcome. The result of a decision tree analysis can be represented in a tree-like structure with nodes representing distinct features and leaves capturing those observations identifed as most similar or “pure.” A purity metric, , where refers to a parent node indexed by and (with its complement ) refers to a specific partition applied to , is a criterion used to determine which of these binary splits will be applied at . Purity measures used previously in the CART literature, for example, include misclassification rates for binary outcomes or residual sum of squares for continuous outcomes.17 In addition to a purity measure, estimation in the tree-based learning context relies upon a set of criteria needed to establish whether a parent node will or will not be partitioned. These criteria, often termed “stopping criteria,” are determined by pre-specified, user-defined inputs and are discussed in the section “PenSIT Learning: Implementation of PenSIT Learning Node Splitting and Stopping Rules.”
Here it is important to note the differences between estimation using CART versus estimation of an optimal DTR using tree-based methods. First, CART is a supervised learner whereas estimation of an optimal DTR is unsupervised. Supervised learning, a term used in the machine learning sphere, refers to estimation of a model in the presence of the outcome/label/flag; unsupervised learning, on the other hand, performs its estimation in the absence of knowledge of the outcome/label/flag. In our case, the true label, that is, the true optimal treatment, is unknown for each subject, and thus our method of estimation is unsupervised. In CART, the object of estimation, which is often an outcome of interest, is directly observed whereas, for optimal DTR estimation, the object of estimation is an optimal treatment sequence that is not directly observed and is therefore considered unsupervised. Secondly, while CART is used for prediction of a target variable based on observed covariates, the goal of optimal DTR estimation represents the crux of personalized medicine: to estimate stage-specific and dynamic decision rules for assigning treatment to individuals based on their unique demographic or disease-specific characteristics such that overall (counterfactual) outcomes for the population will be maximized. Given that optimal DTR estimation has a causal goal, purity measures used previously in this context are derived from estimators of the counterfactual mean outcome; these include an IPW-based estimator,3 an AIPW-based estimator,7 and others. As will be shown in the following sections, we replace the expression identified in IPW-based estimators with the expression introduced in equation (2).
PenSIT estimation for final stage
For simplicity, denote . We propose to model as:
where is a pre-specified transformation of the propensity score, for example, logit or identity; denotes a penalized spline with fixed knots18,19 indexed using the parameters ; and refers to a parametric function of other covariates in , indexed by the parameters . A truncated linear spline basis with fixed knots was used, with equal knot spacing and ,19,20 where refers to the sample size with observations . A penalty term was introduced through the use of linear mixed models,11,19 with random effects for the truncated linear spline basis matrix. Note in equation (3) that the conditional mean outcome is modeled semiparametrically using covariates predictive of the outcome and a penalized spline of a function of the propensity score, which makes the model more flexible and easier to achieve correct model specification. We then propose to estimate the counterfactual mean outcome under treatment , that is, , as , where we introduce the PenSI estimator as , and refers to the empirical mean operator. Following Zhou et al.,11 is consistent and doubly robust for (proof found in Web Appendix A).
Assuming observations are independent and identically distributed across all individuals and following some unspecified multivariate distribution , and also assuming consistency, positivity, and ignorability, is consistent and doubly robust for if either of the following conditions are met:
The conditional mean model consisting of and is correctly specified.
The parametric component of the conditional mean model, , is misspecified but the propensity models and the relationship between the outcome and the function of the propensity score are correctly specified.
The consistency of offers a valid large sample estimate of the counterfactual mean outcome under treatment conditional on prior treatment ; additionally, the double robustness of the PenSI estimator provides two opportunities for estimation consistently. Furthermore, Zhang and Little10 maintain that, by regressing the outcome on a spline of the logit of the propensity score, condition (2) in Proposition 3.1 is met under relatively weak conditions due to the modeling flexibility of a spline, which requires minimal assumptions about the relationship between the outcome and the propensity for treatment assignment, .
Given our goal to estimate a treatment regime that maximizes the counterfactual mean outcome following regime , we estimate using our proposed PenSI estimator as follows:
Based on this formulation we then propose a purity measure, , suitable for constructing a tree when estimating an optimal treatment rule at the final, th stage:
Specifically, refers to the maximum empirical version of the expected counterfactual outcome under decision rule when node is split according to partition such that patients in subset are assigned treatment while patients in the complementary set are assigned to , for .
PenSIT estimation for stages
Similar to other tree-based optimal DTR estimation methods (e.g. Sun and Wang,6 Tao et al.7), estimation proceeds in a backward recursive manner, beginning with estimation of the th stage decision rule. This is important to account for time-varying confounding by indication and effects related to treatments received at earlier stages but mediated through later treatment(s), both of which can result in biased estimation. Therefore, in the context of estimation of an optimal multi-stage DTR, the optimal -stage decision rule relies upon the patient receiving the optimal treatment at all future stages.
Following our exposition in the previous sections, is estimated within the tree-based construct using the PenSI purity measure for the -th stage, , introduced in equation (4). In order to generalize for estimation of the th stage decision rule, for , we now introduce additional notation. Let refer to the predicted counterfactual pseudo-outcome (defined below) at stage under treatment while all future treatments are optimal, which is potentially never actually observed. The assumption under an optimal, multi-stage treatment assignment regime is that the long-term outcome is maximized. Therefore, when estimating the decision rule for the th treatment stage, we must account for the fact that the patient was treated with the optimal treatment at all future stages. To this end we construct a stage-specific pseudo-outcome for any stage prior to the last, which represents the predicted counterfactual outcome at the -th stage contingent upon the patient receiving the optimal treatment at all future stages. Mathematically this can be expressed as: . Under an assumption of consistency, . Positivity in the multi-stage setting was introduced in the “Notation and Formulation: Link to Observed Data” section and the assumption of ignorability can be expressed for the th stage estimation as . Therefore, similar to that introduced in equation (2), the optimal decision rule at the th stage can be expressed as a function of the predicted pseudo-outcome, with :
We define and propose to estimate the mean pseudo-outcome under treatment as , where, using the notation and modeling choices introduced in equation (3):
Assuming consistent estimation at all future stages through backward induction and following Proposition 3.1 and Zhou et al.,11 is a consistent and doubly robust estimator for .
The associated PenSI purity measure used at the th treatment stage, that is, , can then be defined as follows, where node is split by the partition identified by based on rule , which assigns treatment to patients in the set defined by and assigns to those in , for :
Details related to calculation of , a key component of the PenSI purity measure at each stage, are provided in Web Appendix B.
Selection of tuning parameters for tree-based estimation
Several user-defined inputs are needed to implement PenSIT Learning. First, a positive value, , must be specified in order to determine whether a potential split of node by partition identifies a meaningful difference in purity, that is, . We recommend that be selected to represent a level of clinical or practical significance determined based on clinical knowledge or practical rationale, although may also be chosen adaptively from the data.7,17,16,21,22
Two other user-specified tuning parameters are also necessary at each stage to perform PenSIT Learning: the minimum number of observations that can fall into each of the terminal nodes, , and a maximum depth to which the tree is allowed to grow, . Generally, the smaller the minimum node size and the larger the depth, the more complex the estimated optimal decision rule will be, leading to the potential of overfitting.6 An optimal range for the minimum node size in a CART-type analysis between 1 and 20 has been suggested.23 Similarly, a depth of and if often considered a good starting point.21
Implementation of PenSIT Learning node splitting and stopping rules
As discussed above, tree-based estimation is performed using backward induction, commencing with the th stage and ending with Stage . Input for the tree-based partitioning at each stage include estimated counterfactual outcomes for all and user pre-specified , , , as mentioned above. Using the framework of Tao et al.,7 the following terminal criteria determine when a node becomes a terminal node in the th stage estimation:
If the size of a node is less than twice the minimum node size , that is, , becomes a terminal node.
If all possible splits of result in child nodes with fewer than observations, that is, , for all possible partitions , then becomes a terminal node.
If the tree depth reaches the pre-specified depth , all nodes at depth become terminal nodes.
The process of recursively splitting the tree into successively smaller partitions is conducted as follows. Begin with root node , for .
At node , evaluate the three terminal criteria above.
If at least one termination criterion is satisfied, no splits of the node are carried out. Assign a single best treatment to all subjects in : , where refers to the purity in the absence of a split.
If no termination criteria are met, determine the best split as follows: .
– If , no split of the node is carried out. Assign a single best treatment to all subjects in :
– If , split into child nodes and as determined by .
If all nodes are terminal nodes, stop. If not, set and repeat Step (1).
Simulation studies
We consider an observational study for a two-stage DTR with treatment options per stage. For each independent individual we generate a set of or baseline covariates, , from a multivariate normal distribution with a mean of and an exchangeable correlation structure defined by correlation coefficient , where refers to the cardinality, or size, of the vector . First stage treatment is generated from a binary distribution, with the data generating mechanisms reflecting low, moderate, or high degree of confounding (as described in Web Appendix C).
The stage 1 optimal decision rule assuming an underlying tree-type DTR is and when an underlying non-tree-type DTR structure is assumed. The intermediate reward outcome following Stage 1 is generated as , where . This reflects an unequal penalty dependent on one of the observed covariate values if the patient was not treated according to their optimal therapy; this is intended to add an additional degree of complexity into the data generating scenario and be more reflective of data that may be encountered in a real world setting. Second stage treatment is generated from a binomial distribution and also reflects either a low, moderate, or high degrees of confounding; refer also to Web Appendix C. The stage 2 optimal decision rule assuming an underlying tree-type DTR is and, when an underlying nontree-type DTR structure is assumed, . The final outcome where , with . Under optimal treatment allocation and assuming independence across observations, .
We compared the performance of PenSIT Learning to tree-based reinforcement learning (T-RL), Q-Learning using linear modeling (Q-Linear), and Q-Learning using nonparametric modeling (Q-NP). T-RL is a tree-based DTR estimation method that uses a purity measure constructed using the AIPW estimator of the counterfactual mean outcome.7 Q-Learning using linear regression modeling assumes a linear and additive relationship between the covariates and the expected outcome (using the lm function in R). Q-Learning with nonparametric modeling allows a more flexible relationship for the Q-functions, which are estimated using random forest prediction (using randomForest in R). Both T-RL and PenSIT Learning use random forest prediction to generate stage 1 pseudo-outcomes. Across all simulation studies we assume that the parametric component of the conditional mean model in PenSIT Learning, , is incorrectly specified for all in order to evaluate performance under a scenario more reflective of the real world, noting that we would expect improved performance over the results presented when is correctly specified. Performance of PenSIT Learning and T-RL are evaluated under both a correctly- and incorrectly-specified propensity model. Although consistent estimation of the counterfactual mean outcome is not ensured when both the parametric component of the prediction model and the propensity model are misspecified, we present results for an incorrectly-specified propensity model in order to demonstrate performance in this scenario, which may be likely to occur in practice. With a test set of size (), performance is evaluated using two metrics: (1) the percentage of observations correctly classified to their optimal treatments, , and (2) the estimated counterfactual mean outcome , which reflects the expected counterfactual outcome had everyone in the patient population of interest been treated “optimally” based on the regime estimated using each respective method. For each simulation design setting we tabulate the median and interquartile range (IQR) of and across all Monte Carlo iterations. For all stages within each simulation experiment, is set as a 5% improvement although we find in supplemental simulations (results not shown) that PenSIT Learning appears to be more sensitive to perturbations of than does T-RL, particularly with a low degree of confounding. For all simulations we set the minimum node size and tree-depth at 20 and 5, respectively. We tabulate results based on the logit transformation of the estimated propensity score; we also explored the use of the identity transformation and found performance differences of the identity transformation compared to the logit transformation to be negligible (results not shown). We also investigated performance using different knot sizes, from 5 to 35 knots, and also using stepwise model selection based on improvement of AIC; no meaningful differences based on these modifications were found (results not shown), suggesting robustness of this method to the choices of knot size and knot placement.
As can be seen generally across all our simulation results (Tables 1 to 3; Supplemental Tables S1 to S4 in Web Appendix D), PenSIT Learning performance improves with sample size when the number of baseline covariates are held fixed, as expected. When the sample size is fixed, performance generally worsens as the number of baseline covariates increases, although this is most apparent with smaller sample sizes. When the underlying DTR structure is tree-type, PenSIT Learning performs well across all data generating settings (Tables 1 and 3), although performance is best when the level of confounding is moderate or high. When and , for example, PenSIT Learning is able to correctly identify the optimal treatments more than 95% of the time when confounding is moderate or high, but this falls to around 85% correct assignment of the optimal treatments with a low degree of confounding. Additionally, the percentage of correctly-classified treatments for PenSIT Learning is similar—typically within 1% across all covariate cardinality (Table 1). When confounding is low (Supplemental Tables S2 and S4 in Web Appendix D), we observe a high degree of variability for PenSIT Learning in both the estimated percentage of correctly-classified observations and in the estimated counterfactual mean outcome. When and , for example, the interquartile range of the correctly-classified optimal treatments is more than 13%, compared with less than half that value when confounding is moderate or high. With a sample size of (top panels, Table 3, Supplemental Table S4 in Web Appendix D), PenSIT Learning under an assumed tree-type DTR reveals similar performance to that discussed above, exceeding 90% correct classification with moderate and high levels of confounding across all settings. When the underlying DTR is nontree-type, PenSIT Learning performance is modest across all sample sizes (Tables 2 and 3, Supplemental Tables S2 to S4 in Web Appendix D), correctly identifying the optimal treatments between 75% and 80% of the time, with performance improving both as the level of confounding decreases and as the number of covariates decreases.
Performance summary % [median, (IQR)] and [median, (IQR)] for estimation of an optimal two-stage dynamic treatment regime (DTR) with two possible treatments per stage, assuming an underlying tree-type DTR structure with varying degrees of confounding for larger sample sizes of and . Generated with specified training dataset sample size () and test dataset size; No.Var.H number of variables in covariate history ; Propensity model is generated using either “correct” or “incorrect” specification; generated using multivariate normal distribution with using exchangeable correlation structure (); PenSIT Learning Penalized Spline-Involved Tree-based learning; T-RL Tree-based Reinforcement Learning; Q-Linear Linear Q-Learning; Q-NP Nonparametric Q-Learning; % opt percent of test set classified to its optimal treatment using a treatment rule estimated using the applicable method; IQR interquartile range; refers to the estimated counterfactual mean outcome under the estimated optimal DTR. Under optimal treatment allocation .
No.Var.H = 20
No.Var.H = 50
Method
% opt
% opt
Moderate degree of confounding
Q-Linear
70.0 (2.6)
7.0 (0.14)
64.5 (3.2)
6.8 (0.14)
Q-NP
89.1 (6.3)
7.7 (0.18)
84.3 (6.8)
7.6 (0.19)
Correct
T-RL
96.6 (6.3)
7.9 (0.22)
96.0 (5.8)
7.9 (0.22)
PenSIT
96.9 (5.5)
7.9 (0.21)
96.6 (5.1)
7.9 (0.20)
Incorrect
T-RL
96.1 (6.3)
7.9 (0.23)
94.4 (12.7)
7.8 (0.33)
PenSIT
96.9 (5.1)
7.9 (0.20)
96.6 (5.1)
7.9 (0.20)
Q-Linear
72.0 (2.2)
7.1 (0.12)
69.2 (2.1)
7.0 (0.13)
Q-NP
97.0 (2.0)
7.9 (0.10)
96.2 (2.9)
7.9 (0.12)
Correct
T-RL
98.0 (2.5)
7.9 (0.12)
97.5 (2.8)
7.9 (0.14)
PenSIT
98.4 (1.9)
8.0 (0.14)
98.3 (1.9)
7.9 (0.14)
Incorrect
T-RL
97.8 (2.6)
7.9 (0.13)
97.1 (3.0)
7.9 (0.14)
PenSIT
98.4 (1.9)
8.0 (0.12)
98.3 (1.8)
7.9 (0.14)
Higher degree of confounding
Q-Linear
69.6 (2.7)
7.0 (0.14)
64.7 (3.1)
6.8 (0.15)
Q-NP
79.7 (4.7)
7.5 (0.18)
77.6 (4.4)
7.4 (0.18)
Correct
T-RL
92.1 (15.4)
7.8 (0.42)
92.9 (13.9)
7.8 (0.38)
PenSIT
96.9 (4.5)
7.9 (0.18)
96.6 (5.4)
7.9 (0.21)
Incorrect
T-RL
91.0 (16.5)
7.7 (0.42)
88.0 (15.1)
7.7 (0.38)
PenSIT
96.9 (4.7)
7.9 (0.18)
97.0 (4.7)
7.9 (0.18)
Q-Linear
71.6 (2.3)
7.1 (0.12)
69.1 (2.3)
7.0 (0.13)
Q-NP
89.4 (9.2)
7.8 (0.25)
85.1 (8.8)
7.6 (0.26)
Correct
T-RL
95.4 (13.0)
7.8 (0.28)
95.5 (12.9)
7.9 (0.32)
PenSIT
98.5 (1.7)
7.9 (0.12)
98.2 (2.2)
7.9 (0.12)
Incorrect
T-RL
95.9 (11.2)
7.9 (0.28)
93.9 (16.6)
7.8 (0.41)
PenSIT
98.5 (1.7)
7.9 (0.12)
98.3 (2.0)
7.9 (0.12)
Performance summary % [median, (IQR)] and [median, (IQR)] for estimation of an optimal two-stage dynamic treatment regime (DTR) with two possible treatments per stage, assuming an underlying nontree-type DTR structure with varying degrees of confounding for larger sample sizes of and . Generated with specified training dataset sample size () and test dataset size; No.Var.H number of variables in covariate history ; Propensity model is generated using either “correct” or “incorrect” specification; generated using multivariate normal distribution with using exchangeable correlation structure (); PenSIT Learning Penalized Spline-Involved Tree-based learning; T-RL Tree-based Reinforcement Learning; Q-Linear Linear Q-Learning; Q-NP Nonparametric Q-Learning; % opt percent of test set classified to its optimal treatment using a treatment rule estimated using the applicable method; IQR interquartile range; refers to the estimated counterfactual mean outcome under the estimated optimal DTR. Under optimal treatment allocation .
Method
No.Var.H = 20
No.Var.H = 50
% opt
% opt
Moderate degree of confounding
Q-Linear
78.0 (4.0)
7.3 (0.17)
69.2 (3.9)
7.0 (0.17)
Q-NP
80.0 (4.6)
7.4 (0.16)
75.8 (4.9)
7.3 (0.15)
Correct
T-RL
76.4 (6.1)
7.2 (0.20)
76.0 (6.8)
7.2 (0.21)
PenSIT
78.0 (4.6)
7.3 (0.16)
78.2 (4.9)
7.3 (0.16)
Incorrect
T-RL
75.8 (6.5)
7.2 (0.21)
75.0 (6.7)
7.2 (0.22)
PenSIT
78.0 (4.6)
7.3 (0.16)
78.1 (4.8)
7.3 (0.16)
Q-Linear
82.4 (3.2)
7.4 (0.13)
76.5 (3.0)
7.3 (0.16)
Q-NP
83.8 (2.9)
7.5 (0.12)
81.1 (3.2)
7.5 (0.14)
Correct
T-RL
77.1 (5.3)
7.3 (0.15)
76.9 (5.0)
7.3 (0.16)
PenSIT
78.7 (4.0)
7.3 (0.14)
78.2 (3.7)
7.3 (0.15)
Incorrect
T-RL
77.2 (5.1)
7.3 (0.15)
77.0 (5.3)
7.3 (0.17)
PenSIT
78.7 (4.1)
7.3 (0.14)
78.3 (3.8)
7.3 (0.15)
Higher degree of confounding
Q-Linear
73.6 (5.5)
7.2 (0.21)
65.3 (5.2)
6.9 (0.21)
Q-NP
75.2 (5.4)
7.3 (0.19)
70.9 (5.5)
7.2 (0.20)
Correct
T-RL
73.2 (6.7)
7.2 (0.23)
72.8 (7.8)
7.1 (0.27)
PenSIT
75.7 (5.5)
7.2 (0.16)
75.4 (5.2)
7.2 (0.16)
Incorrect
T-RL
72.6 (7.4)
7.1 (0.27)
72.2 (8.5)
7.1 (0.26)
PenSIT
75.4 (5.4)
7.2 (0.15)
75.6 (5.4)
7.2 (0.17)
Q-Linear
78.7 (4.0)
7.3 (0.16)
73.1 (3.8)
7.1 (0.16)
Q-NP
80.2 (3.5)
7.5 (0.13)
77.4 (3.9)
7.4 (0.14)
Correct
T-RL
73.7 (6.3)
7.2 (0.20)
74.2 (6.2)
7.2 (0.19)
PenSIT
76.0 (4.1)
7.3 (0.15)
76.4 (4.4)
7.3 (0.14)
Incorrect
T-RL
74.6 (6.3)
7.2 (0.19)
74.3 (6.2)
7.2 (0.18)
PenSIT
76.0 (4.2)
7.3 (0.15)
76.2 (4.4)
7.3 (0.14) 6
Performance summary % [median, (IQR)] and [median, (IQR)] for estimation of an optimal two-stage dynamic treatment regime (DTR) with two possible treatments per stage to evaluate performance in smaller samples () for both tree- and nontree-type DTRs with varying degrees of confounding. Generated with training dataset sample of size with test dataset size; No.Var.H number of variables in covariate history ; Propensity model is generated using either “correct” or “incorrect” specification; generated using multivariate normal distribution with using exchangeable correlation structure (); PenSIT Learning Penalized Spline-Involved Tree-based learning; T-RL Tree-based Reinforcement Learning; Q-Linear Linear Q-Learning; Q-NP Nonparametric Q-Learning; % opt percent of test set classified to its optimal treatment using a treatment rule estimated using the applicable method; IQR interquartile range; refers to the estimated counterfactual mean outcome under the estimated optimal DTR. Under optimal treatment allocation .
Method
No.Var.H = 20
No.Var.H = 50
% opt
% opt
Tree-type DTR
Moderate degree of confounding
Q-Linear
67.5 (3.3)
6.9 (0.15)
58.9 (4.0)
6.6 (0.16)
Q-NP
78.0 (7.4)
7.4 (0.23)
72.6 (7.0)
7.2 (0.26)
Correct
T-RL
94.6 (12.0)
7.8 (0.34)
93.1 (13.0)
7.8 (0.34)
PenSIT
94.7 (11.2)
7.8 (0.27)
93.2 (12.7)
7.8 (0.31)
Incorrect
T-RL
93.4 (14.8)
7.8 (0.38)
86.8 (17.2)
7.6 (0.46)
PenSIT
94.7 (11.4)
7.8 (0.28)
93.2 (12.9)
7.8 (0.31)
Higher degree of confounding
Q-Linear
67.1 (3.7)
6.9 (0.17)
58.9 (4.4)
6.6 (0.16)
Q-NP
73.8 (5.8)
7.3 (0.20)
70.0 (6.0)
7.2 (0.22)
Correct
T-RL
89.4 (15.5)
7.7 (0.43)
86.6 (18.2)
7.6 (0.47)
PenSIT
94.6 (8.5)
7.8 (0.25)
91.7 (11.7)
7.7 (0.32)
Incorrect
T-RL
86.1 (19.8)
7.7 (0.50)
86.4 (11.3)
7.6 (0.25)
PenSIT
94.4 (8.6)
7.8 (0.27)
92.4 (11.4)
7.8 (0.30)
Nontree-type DTR
Moderate degree of confounding
Q-Linear
73.0 (4.8)
7.1 (0.20)
61.7 (5.0)
6.7 (0.22)
Q-NP
75.4 (5.6)
7.3 (0.18)
69.4 (7.3)
7.1 (0.24)
Correct
T-RL
74.7 (7.8)
7.2 (0.23)
74.5 (8.2)
7.2 (0.25)
PenSIT
77.2 (6.1)
7.3 (0.17)
77.0 (6.3)
7.3 (0.19)
Incorrect
T-RL
74.3 (8.0)
7.2 (0.25)
72.0 (9.0)
7.1 (0.31)
PenSIT
77.0 (6.3)
7.3 (0.17)
76.5 (6.3)
7.2 (0.19)
Higher degree of confounding
Q-Linear
68.6 (6.1)
7.0 (0.22)
57.6 (6.0)
6.6 (0.25)
Q-NP
69.9 (7.8)
7.1 (0.27)
62.8 (8.9)
6.9 (0.36)
Correct
T-RL
71.1 (9.3)
7.1 (0.33)
70.3 (10.1)
7.0 (0.34)
PenSIT
75.0 (6.4)
7.2 (0.18)
74.3 (7.1)
7.2 (0.21)
Incorrect
T-RL
71.7 (11.3)
7.1 (0.36)
61.6 (22.8)
6.8 (0.71)
PenSIT
75.1 (6.4)
7.2 (0.19)
74.4 (7.9)
7.2 (0.24)
For a fixed and a true, underlying tree-type DTR, PenSIT Learning is preferred to competing methods when the level of confounding is high. Furthermore, the improvement of PenSIT Learning over other methods is most pronounced when the covariate cardinality is large or when the sample size is small. With moderate confounding, PenSIT Learning performance is comparable to T-RL across all settings and is comparable to Q-NP when , although PenSIT Learning achieves a clear advantage over T-RL when the sample size is small and the propensity model is incorrectly specfied. When the level of confounding is low, T-RL is preferred over PenSIT Learning for a fixed across all sample sizes. With a nontree-type DTR structure, PenSIT Learning is generally preferred across all data generating scenarios, including different levels of confounding settings, although the improvement over T-RL is modest, with both methods correctly classifying the optimal treatment regime between 75% and 80% of the time. Under low confounding and low covariate cardinality with a nontree-type DTR, Q-NP would be preferred to both PenSIT Learning and T-RL. Relative to PenSIT Learning, Q-Linear exhibits inferior performance across all sample sizes and levels of confounding when the true DTR structure is tree-type. When the sample size is small and the DTR is nontree-type, PenSIT Learning achieves slightly better optimal treatment classification rates than does Q-NP and Q-Linear, particularly when the number of covariates increases.
Application of PenSIT Learning to MIMIC-III data
Sepsis is a clinical syndrome characterized by systemic inflammation and infection and is associated with one of the highest rates of mortality among conditions commonly treated in emergency departments (EDs) and intensive care units (ICUs).24 Sepsis is routinely treated using fluid resuscitation, antibiotics, and may also include treatment with vasopressors, mechanical ventilation, and others. The established clinical guidelines for treating sepsis, known as the “Surviving Sepsis Campaign,”25 strongly recommend that resuscitation of at least 30 mL/kg of intravenous (IV) fluid be given within the first 3 h. However, this recommendation is given with a stated “low quality of evidence” due to the fact that results across studies have been inconsistent with indirect evidence, imprecise results, and a likelihood of bias. Therefore, due to the paucity of strong evidence as to the most beneficial fluid resuscitation strategy in the early hours of treatment, using electronic medical record data from the Medical Information Mart for Intensive Care III (MIMIC-III),26 we estimate an optimal two-stage DTR in adult septic patients admitted to the medical ICU (MICU) after presenting to the ED with the goal of minimizing a measure of multi-organ failure. Refer to additional information about cohort eligibility and data analysis in Web Figure 1 and Web Appendix D. As MIMIC-III is a publicly available dataset, the construction and de-identification of which was approved by the institutional review boards (IRBs) for the owner and host of these data,26 there was a need for neither further local IRB approval nor subject-level informed consent for this research.
The estimated two-stage treatment strategy to optimize the patient-level Sequential Organ Failure Assessment (SOFA) score evaluated at 24 h following admission. Our proposed method suggests that all patients should receive a high volume fluid resuscitation strategy ( mL/kg) within 3 h after admission to the Medical Intensive Care Unit (MICU) from the emergency department (ED). If the patient receives high volume fluid resuscitation within the first 3 h following admission in accordance with this strategy, they should receive low volume fluid resuscitation ( mL/kg) between 3 and 24 h following MICU admission. If they did not receive an initial high volume resuscitation strategy in accordance with the estimated guideline, however, they should receive high volume ( mL/kg) fluid resuscitation between 3 and 24 h following MICU admission. mL/kg = milliliters per kilogram.
Baseline covariates considered as candidate tailoring variables for treatment strategies included demographics such as gender, age, weight, and racial/ethnic groups, Elixhauser comorbidity score, and the time of year in which the patient was treated. Stage 1 treatment was defined as either a fluid restrictive ( mL/kg) or a fluid liberal ( mL/kg) strategy within the first 3 h after admission to the MICU. Intermediate variables collected prior to Stage 2 treatment included indicators of treatment with mechanical ventilation and vasopressors within the first 3 h time period, as well as the patient’s Sequential Organ Failure Assessment (SOFA) score evaluated at 3 h post-admission. Stage 2 treatment was defined as either a fluid restrictive ( mL/kg) or a fluid liberal ( mL/kg) strategy between 3 and 24 h after MICU admission. The primary outcome of interest is the SOFA score, used as a predictor of sepsis-related mortality, evaluated at 24 h post-admission.
Four hundred eighty-six (486) patients were included in the analysis cohort. Seventy-six percent (76%) of the study cohort are white and 52% are male, with an average age of 69 years and a majority reporting Elixhauser comorbidities (Table 4). The median length of hospital stay was 7.8 days with an interquartile range (IQR) of 5.0–13.8. The median fluid input received within 0–3 h and 3–24 h post-admission is 41.4 mL/kg (IQR: 22.8–60.9) and 20.2 mL/kg (IQR: 3.6–52.4), respectively. Summary statistics stratified by treatment stage (i.e. 0–3 h and 3–24 h post-MICU admission) demonstrate covariate imbalance for age, gender, and weight across fluid resuscitation strategies in the first treatment stage, for race/ethnicity at the second stage, and for the use of mechanical ventilation and vasopressors across fluid resuscitation strategies for both stages. Additionally, evaluation of our observed data demonstrates strong associations between total fluid intake, as well as the use of ventilation and vasopressors, with SOFA score. Both findings, as well as medical knowledge underpinning this relationship, suggest that confounding is a major issue that must be addressed in our analysis in order to make valid causal inference.
Characteristics of the analysis cohort. Summary statistics of demographics, treatment, and outcomes for MIMIC-III analysis cohort are included. n sample size; Stage 1 0–3 h post-admission; Stage 2 3–24 h post-admission; Restrictive fluid resuscitation ( mL/kg); Liberal fluid resuscitation ( mL/kg); IQR interquartile range; kg kilogram; LOS length of hospital stay; Mech Vent Mechanical Ventilation; Vasos Vasopressors; L liters; mL/kg milliliters per kilogram; SOFA sequential organ failure assessment; Median [IQR] are presented for continuous variables; frequency (percentage) are provided for categorical variables.
Overall (n=486)
Stage 1
Stage 2
Restrictive (n=163)
Liberal (n=323)
Restrictive (n=289)
Liberal (n=197)
Patient characteristics
Age (years)
69 [54–82]
71 [57–81]
68 [53–82]
69 [55–82]
69 [54–82]
Gender
Male
252 (52)
94 (58)
158 (49)
149 (52)
103 (52)
Female
234 (48)
69 (42)
165 (51)
140 (48)
96 (48)
Race/ethnicity
White
371 (76)
123 (76)
248 (77)
226 (78)
145 (74)
Nonwhite
115 (24)
40 (24)
75 (23)
63 (22)
52 (26)
Weight (kg)
77 [65–91]
83 [69–98]
74 [62–87]
80 [68–94]
72 [62–85]
LOS (days)
7.8 [5.0–13.8]
8.5 [5.0–15.7]
7.7 [5.0–12.7]
7.3 [4.8–12.0]
8.2 [5.8–14.8]
0–3 h post-admission
Use of Mech Vent
111 (23)
43 (26)
68 (21)
74 (26)
37 (19)
Use of Vasos
89 (18)
16 (10)
73 (23)
50 (17)
39 (20)
Total input (L)
3.0 [2.0–5.0]
1.2 [0.9–2.0]
4.0 [3.0–5.2]
2.9 [1.5–4.0]
4.0 [2.5–5.3]
Total input (mL/kg)
41.4 [22.8–60.9]
16.7 [11.2–22.9]
53.8 [41.4–71.4]
35.2 [17.4–51.9]
53.3 [34.0–74.4]
SOFA (3 h)
4 [2–6]
4 [2–6]
4 [2–6]
5 [3–6]
4 [2–6]
3–24 h post-admission
Use of Mech vent
197 (41)
69 (42)
128 (40)
99 (34)
98 (50)
Use of Vasos
189 (39)
49 (30)
140 (43)
80 (28)
109 (55)
Total input (L)
2.5 [1.0–4.5]
1.5 [1.0–2.7]
3.2 [1.3–5.5]
1.0 [0.7–1.6]
4.5 [3.3–6.0]
Total input (mL/kg)
20.2 [3.6–52.4]
10.7 [0.0–25.2]
30.3 [9.0–62.7]
7.2 [0.0–16.9]
58.9 [44.0–87.3]
SOFA (24 h)
5 [3–7]
5 [3–7]
5 [3–8]
5 [3–6]
6 [3–9]
As shown in Figure 1, it is recommended that all patients receive liberal fluid resuscitation ( mL/kg) within the first three hours following admission to the MICU for treatment of acute emergent sepsis. If the patient has received the liberal fluid resuscitation by 3 hours post-admission in accordance with this estimated decision rule, the patient should receive restrictive fluid resuscitation ( mL/kg) to follow. If the patient was not given liberal fluid resuscitation within the first three hours following admission, the patient should receive liberal fluid resuscitation within 3–24 h post-admission. Using our proposed method, we found no tailoring variables at the first treatment stage that would result in a meaningful improvement in outcomes overall, while the second stage treatment, however, can be tailored based on the patient’s first-stage treatment in order to optimize counterfactual outcomes overall. After accounting for the transformation of the outcome used in the regression models, this represents a SOFA improvement of roughly , which we consider to be clinically meaningful.
Although the question of how to optimally treat septic patients is complex and multi-faceted, we applied a robust and flexible causal method with interpretable results to determine whether tailoring of fluid resuscitation strategies at each of two stages within the first 24 h after MICU admission can be used to improve outcomes overall. Consistent with the Surviving Sepsis Campaign best practice recommendations, our method suggests that liberal fluid resuscitation should be given as early as possible in this patient population in order to reduce early indicators of organ dysfunction.
Discussion
PenSIT Learning retrofits a decision tree with a novel PenSI purity measure that incorporates the estimated propensity for treatment assignment as a spline predictor rather than a weight.11 Not only does PenSIT Learning retain the flexibility of T-RL and other tree-based optimal DTR estimation methods, but it provides added robustness under conditions of a high degree of confounding. Additionally, the PenSI estimator of the counterfactual mean utilized within the PenSI purity measure fulfills the properties of consistency and double robustness in asymptotia under standard regularity conditions.
There are several distinct advantages of PenSIT Learning. First, PenSIT Learning utilizes a decision tree construct and provides an interpretable estimate of the underlying DTR, which we find of great importance for communicating effectively with clinician-scientists. Secondly, the PenSI estimator of counterfactual outcomes is derived using standard regression models for the treatment assignments and the conditional outcomes at all stages, suggesting that the standard knowledge base surrounding regression models, including model building and selection, model fit diagnostics, etc., can and should be applied liberally. Moreover, although our simulation experiments were designed such that the same variables were used to define both counterfactual outcomes within each stage, PenSIT Learning allows modeling choices for the counterfactual outcomes within a stage to differ, which makes practical sense as there is no reason why we would expect the true mechanisms to be the same. Additionally, although we focused on a continuous outcome that is approximately symmetric, the simplicity of PenSIT Learning makes it straightforward to make adjustments to regression models based on the scale of the outcome, for example, using a generalized linear model framework. Finally, we observed distinctly improved performance of PenSIT Learning over competing methods when the level of confounding is high and when the sample size was small, but reduced performance under a low degree of confounding. Empirically, under low confounding the estimated propensity scores for both stages are symmetric and centered at , compared to U-shaped distributions favoring the extremes under high confounding. A higher-confounding relationship appears to offer an enhanced ability to predict counterfactual outcomes, as we had conjectured—perhaps because the rank order of the treatment effect is maintained but the larger differentiation across the logit of the estimated propensity scores facilitates improved prediction.
Several cautions should be heeded when applying PenSIT Learning. First, whereas using a spline function of the estimated propensity score in the conditional mean model allows for a flexible relationship between the propensity score and the outcome, the propensity model should, in theory, be correctly specified. And, in fact, it is often extremely challenging to devise a correctly-specified propensity model given the complexities that abound in a medical treatment setting. In practice, however, simulation results suggest that the assumption of a correctly-specified propensity model is perhaps less critical. Second, our method is based on fulfilling the assumptions of consistency, positivity, and ignorability. The consistency assumption is reasonable in many experimental settings. Positivity can also generally be justified using content knowledge and by covariate balance diagnostics conditional on the propensity score; however, in the setting with a large number of covariates and/or multiple treatments and stages, some accommodations to ensure positivity may be needed.27–30 Ignorability, conversely, may not always be reasonable in an observational data setting, although investigators may be willing to proceed under an assumption of ignorability due to the fact that optimal DTR estimation is largely exploratory and should be challenged in confirmatory studies. Lastly, supplemental simulation studies reveal that PenSIT Learning may be particularly sensitive to the choice of , especially when there is a lower degree of confounding of the relationship between the treatment and the outcome. In practice, minimal association between a covariate and both the treatment and outcome of interest—indicative of a low degree of confounding—could be identified based on scientific judgment, and statistical evaluation of these association using the observed data. However, we do generally expect a moderate to high degree of confounding when using observational data to evaluate causal effects in a medical setting and we maintain the view that the choice of and other tuning parameters should be guided by scientific knowledge.
There are several extensions to PenSIT Learning that we believe would be of interest to the research community. The first would be to explore the performance of PenSIT Learning using more flexible modeling approaches for the treatment assignment and/or the conditional counterfactual outcome models (e.g. using BART, random forests, kernel-based methods, etc.). It is possible that these methods may improve performance under more complex data generation settings, and may be more desirable for accommodating an abundance of data from which there are few a priori patterns or insights. Second, although we believe that the choice of should be driven based on scientific knowledge, in the absence of information about a suitable , data-driven approaches for selecting the tuning parameters could be explored. Finally, extensions of PenSIT Learning to account for potential overfitting inherent in decision tree-type constructs, for example, with stochastic tree search or incorporating a lookahead procedure, can be considered.
Supplemental Material
sj-pdf-1-smm-10.1177_09622802221122397 - Supplemental material for Penalized Spline-Involved Tree-based (PenSIT) Learning for estimating an optimal dynamic treatment regime using observational data
Supplemental material, sj-pdf-1-smm-10.1177_09622802221122397 for Penalized Spline-Involved Tree-based (PenSIT) Learning for estimating an optimal dynamic treatment regime using observational data by Kelly A Speth, Michael R Elliott, Juan L Marquez and Lu Wang in Statistical Methods in Medical Research
Footnotes
Acknowledgements
The authors would like to thank the MIMIC-III and Physionet groups for providing the dataset and supporting materials.
Authors’ note
Thank you to Statistical Methods in Medical Research for considering this manuscript.
Author contributions
KAS, MRE, and LW developed the statistical methods. JLS developed the research question and provided medical guidance for the data analysis. KAS drafted the manuscript. All authors reviewed and edited the manuscript.
Declaration of conflicting interests
The authors declare that there are no conflicts of interest.
Funding
The authors received no financial support for the research, authorship, and/or publication of this article.
Data accessibility statement
The Medical Information Mart for Intensive Care III (MIMIC-III) data are available through Physionet (mimic.physionet.org). Relevant code used to perform simulations and data application will be available on github.
ORCID iD
Kelly A Speth
Supplemental material
Supplemental material for this article is available online.
References
1.
ChakrabortyBMoodieEM. Statistical methods for dynamic treatment regimes: reinforcement learning, causal inference, and personalized medicine. Statistics for Biology and Health. New York: Springer, 2013.
2.
TsiatisAADavidianMHollowayST, et al. Dynamic treatment regimes: statistical methods for precision medicine. Boca Raton, FL: CRC Press, 2020.
3.
LaberEBZhaoYQ. Tree-based methods for individualized treatment regimes. Biometrika2015; 102: 501–514.
4.
ZhangYLaberEBTsiatisAA, et al. Using decision lists to construct interpretable and parsimonious treatment regimes. Biometrics2015; 71: 895–904.
5.
ZhaoYZengDLaberEB, et al. New statistical learning methods for estimating optimal dynamic treatment regimes. J Am Stat Assoc2015; 110: 583–598.
6.
SunYWangL. Stochastic tree search for estimating optimal dynamic treatment regimes. J Am Stat Assoc2021; 116: 421–432.
7.
TaoYWangLAlmirallD. Tree-based reinforcement learning for estimating optimal dynamic treatment regimes. Ann Appl Stat2018; 12: 1914–1938.
8.
ZhangYLaberEBDavidianM, et al. Interpretable dynamic treatment regimes. J Am Stat Assoc2018; 113: 1541–1549.
9.
LittleRAnH. Robust likelihood-based analysis of multivariate data with missing values. Stat Sin2004; 14: 949–968.
10.
ZhangGLittleR. Extensions of the penalized spline of propensity prediction method of imputation. Biometrics2009; 65: 911–918.
11.
ZhouTElliottMRLittleRJA. Penalized spline of propensity methods for treatment comparison. J Am Stat Assoc2019; 114: 1–19.
12.
RubinDB. Estimating causal effects of treatments in randomized and nonrandomized studies. J Educ Psychol1974; 66: 688–701.
13.
RobinsJMHernanMA. Estimation of the causal effects of time-varying exposures. In: Fitzmaurice and Davidian (eds) Longitudinal Data Analysis. Boca Raton, FL: Chapman and Hall/CRC, 2009.
14.
RosenbaumPRRubinDB. The central role of the propensity score in observational studies for causal effects. Biometrics1983; 70: 41–55.
BreimanLFreidmanJHOlshenRA, et al. Classification and regression trees. Belmont, CA: Wadsworth, 1984.
17.
HastieTTibshiraniRFriedmanJ. The elements of statistical learning: data mining, inference, and prediction. 2nd ed. New York, NY: Springer, 2009.
18.
EilersPHCMarzBD. Flexible smoothing with B-splines and penalties. Statist Sci1996; 11: 89–102.
19.
WandMP. Smoothing and mixed models. Comput Stat2003; 18: 223–249.
20.
RuppertD. Selecting the number of knots for penalized splines. J Comput Graph Stat2002; 11: 735–757.
21.
BoehmkeBGreenwellB. Decision Trees. In: Hands-On machine learning with R 2020.
22.
TherneauTMAtkinsonEJ and Mayo Foundation. An introduction to recursive partitioning using the RPART routines. CRAN R Network, 2019.
23.
MantovaniRGHorvathTCerriR, et al. An empirical study on hyperparameter tuning of decision trees. arXiv 2019; 1812.02207v2.
24.
MarinoPL. Marino’s The ICU Book. 4th ed. Philadelphia, PA: Wolters Kluver Health/Lippincott Williams & Wilkins, 2014.
25.
RhodesAEvansLEWaleedA, et al. International guidelines for management of sepsis and septic shock: 2016. Crit Care Med2017; 45: 486–552.
26.
JohnsonAPollardTShenL, et al. MIMIC-III, a freely accessible critical care database. Sci Data2016; 3: 160035.
27.
CrumpRKHotzVJImbensGW, et al. Dealing with limited overlap in estimation of average treatment effects. Biometrika2009; 96: 187–199.
28.
GutmanRRubinDB. Estimation of causal effects of binary treatments in unconfounded studies. Stat Med2015; 34: 3381–3398.
29.
HoDEImaiKKingG, et al. Matching as nonparametric preprocessing for reducing model dependence in parametric causal inference. Polit Anal2007; 15: 199–236.
30.
RosenbaumPR. Optimal matching of an optimally chosen subset in observational studies. J Comput Graph Stat2012; 21: 57–71.
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.