Who will benefit from antidepressants in the acute treatment of bipolar depression? A reanalysis of the STEP-BD study by Sachs et al. 2007, using Q-learning

Background There is substantial uncertainty regarding the efficacy of antidepressants in the treatment of bipolar disorders. Methods Traditional randomized controlled trials and statistical methods are not designed to discover if, when, and to whom an intervention should be applied; thus, other methodological approaches are needed that allow for the practice of personalized, evidence-based medicine with patients with bipolar depression. Results Dynamic treatment regimes operationalize clinical decision-making as a sequence of decision rules, one per stage of clinical intervention, that map patient information to a recommended treatment. Using data from the acute depression randomized care (RAD) pathway of the Systematic Treatment Enhancement Program for Bipolar Disorder (STEP-BD) study, we estimate an optimal dynamic treatment regime via Q-learning. Conclusions The estimated optimal treatment regime presents some evidence that patients in the RAD pathway of STEP-BD who experienced a (hypo)manic episode before the depressive episode may do better to forgo adding an antidepressant to a mandatory mood stabilizer. Electronic supplementary material The online version of this article (doi:10.1186/s40345-014-0018-5) contains supplementary material, which is available to authorized users.

. Furthermore, there is concern that antidepressants can induce abnormal mood elevation (Licht et al. 2008). We use data from the Systematic Treatment Enhancement Program for Bipolar Disorder (STEP-BD) (Sache et al. 2003(Sache et al. , 2007 to estimate an optimal dynamic treatment regime (DTR) (Chakraborty and Murphy 2014;Murphy 2003, Robins 2004, Schulte et al. 2014, for bipolar depression. A DTR is a sequence of decision rules, one per stage of intervention, that map up-to-date patient information to a recommended treatment; thus, an estimated optimal DTR can be used to generate hypotheses about how patient history and outcomes should dictate treatment selection. The estimated optimal DTR for bipolar depression constructed from the STEP-BD study suggests the hypothesis that standard antidepressants should not be used to supplement mood stabilizers for patients with a prior hypomanic episode. A DTR aims to select if, when, what, and to whom treatment should be assigned and thereby fits into the paradigm of personalized medicine. Because DTRs select treatment according to the uniquely evolving health status of each patient, they are suited to manage chronic illnesses with patient response heterogeneity; thus, DTRs have tremendous potential for personalizing and improving treatment strategies for bipolar disorder (Leboyer and Kupfer 2010;Nierenberg et al. 2013). Optimal DTRs have been estimated for wide range of chronic illnesses including major depressive disorder Chakraborty and Moodie 2013), attention deficit hyperactivity disorder Lei et al. 2012, Nahum-Shani et al. 2012a), schizophrenia Shortreed et al. 2011), HIV/AIDS (Moodie et al. 2007;Sterne et al. 2009), and cigarette addiction (Strecher et al. 2006). Estimation of an optimal DTR is typically done as a secondary, exploratory analysis and viewed as a method of generating hypotheses for follow-up confirmatory experiments . This is the perspective we take here; nevertheless, we show that an estimated optimal DTR appears to perform markedly better than any fixed treatment strategy.
In the "STEP-BD study" section, we review the STEP-BD study. In the "Dynamic treatment regimes and Q-learning" section, we formalize DTRs and introduce the Q-learning estimation algorithm. In the "Analysis of STEP-BD" section, we present an analysis of STEP-BD.

Methods
The study on which our analyses are based was approved by the institutional review board at each study site and was overseen by a data and safety monitoring board (for more details, see http://www.ncbi.nlm.nih.gov/pubmed/ 17392295).

STEP-BD study
STEP-BD is a long-term study of bipolar disorder funded by the National Institute of Mental Health (NIMH). Its aim was "to generate externally valid answers to treatment effectiveness questions related to bipolar disorder" (Sachs et al. 2003). Patients of age older than 15 years fulfilling DSM-IV criteria for any subtype of bipolar disorders could enter the study registry. In total, 4,360 patients from 22 sites in United States enrolled. The study lasted for 7 years (2001)(2002)(2003)(2004)(2005)(2006)(2007). In STEP-BD, there are two different treatment pathways: standard care pathway (SCP) and randomized care pathway (RCP). SCP is open to all participants with a diagnosis of bipolar disorders. Each treatment delivered is open and will follow treatment guidelines. Decisions are made on the basis of shared decision-making. After patients signed informed consent for entry into the STEP-BD study registry, all patients enter the SCP. If a patient's status meets the eligibility criteria at one of the follow-up visits during the SCP for a study within the RCP, additional consent is requested for entry into that RCP. The RCP utilizes methods appropriate for efficacy studies, and random assignment is needed to provide answers to clinical questions. In the RCP, there are three different pathways each addressing unmet needs in the treatment of bipolar disorder: acute depression randomized (RAD) pathway, acute depression psychosocial intervention (PAD) pathway, and refractory depression (REFD) pathway. If patients are unwilling to consent to one of the RCPs, they remain in the SCP. In general, the decision of pathway (SCP versus RCP) is based on both the doctor's and patient's opinion. In STEP-BD, patients could switch pathways based on doctor's or their own preference as well as inclusion and exclusion criteria. Figure 1 shows the diagram of STEP-BD study. Our analysis utilizes patients enrolled in RAD.

Acute depression randomized pathway
As mentioned above, the RAD pathway is one of the three RCPs in STEP-BD. In addition to satisfying the general entry criteria for STEP-BD study registry, patients had to be at least 18 years old and fulfill the DSM-IV criteria for a major depressive episode in the context of bipolar I or bipolar II disorder. All patients with a history of intolerance or non-response to both bupropion and paroxetine were excluded, as well as those requiring current shortterm treatment for a coexisting substance-abuse disorder or requiring the addition of antipsychotic medication or a change in the dose of a long-term antipsychotic medication (Sachs et al. 2007). In addition, patients had to take a mood stabilizer at the time of randomization or agree to begin treatment with a mood stabilizer. Moreover, they had to agree to have all non-study antidepressants Figure 1 STEP-BD Design. There were 4,360 patients who participated in SCP and 461 patients who participated in RCP. Patients in SCP satisfying RCP criteria entered RCP. In RCP, there are three different pathways: RAD, PAD, and REFD. Patients can enter more than one random pathways. tapered after initiation of study drug, with the antidepressant discontinued by the end of week 2. The purpose of RAD was to explore the effectiveness of adjunctive antidepressant treatment, in addition to a mood stabilizer. Initially, the mood stabilizers were limited to lithium, valproate, the combination of lithium and valproate, or carbamazepine. However, later on, any FDA-approved antimanic agent could be used as mood stabilizers. Table 1 shows the percentages of different mood stabilizers used in RAD. At week 0, patients were randomly assigned to one antidepressant (150 mg of a sustained release formulation of bupropion or 10 mg of paroxetine to begin with) or placebo. After 6 weeks, patients with non-response on the placebo were randomized to either paroxetine or bupropipon; patients with non-response on the antidepressant were assigned to either openly increase the dose of their current antidepressant or add another antidepressant. At weeks 8, 10, or 12, clinicians will make final decision for patients based on their clinical status collected from clinical monitoring form (CMF). During the study, patients need to visit their doctors every week to fill in the CMF (Sachs et al. 2002). Patients were allowed to switch to SCP (opt out) at any time by their preference or doctor's opinion. Patients who had severe adverse effects or met criteria for hypomania or mania discontinued the antidepressant or placebo and received open treatment while remaining in STEP-BD. Since after 6 weeks, only one patient with non-response on the antidepressant was assigned to add another antidepressant, we ignored this one observation and supposed patients with non-response to antidepressant after 6 weeks were only assigned to increase the dose of their current antidepressant. Figure 2 shows a schematic for the RAD protocol. Response for a given subject was defined as at least 50% improvement over their initial SUMD (scale scores for depression) and not meeting the DSM-IV criteria for hypomania or mania. Scores on the continuous symptom subscale for depression (SUMD) range from 0 to 22, with higher Note that some patients received more than one mood stabilizer.
scores indicating more severe symptoms. Both SUMD and SUMM (symptom subscale for mood elevation, SUMM scores range from 0 to 16) are part of the modified clinical monitoring form for mood disorders (Sachs et al. 2002). Because subjects in RAD are potentially randomized multiple times with randomizations occurring at crucial points of the disease process, RAD is an example of a Sequential Multiple Assignment Randomized Trial (SMART) (Lavori and Dawson 2004;Murphy 2003). Data collected in SMARTs can be used to efficiently estimate and evaluate DTRs. In the next section, we formalize the notion of an optimal DTR and introduce a regressionbased approach called Q-learning for estimating an optimal DTR from a SMART.

Dynamic treatment regimes and Q-learning
The effective management of a chronic illness requires ongoing personalized treatment (Wagner et al. 2001). DTRs formalize clinical decision-making as sequence of decision rules, one per treatment decision, that map patient information to a recommended treatment. An optimal DTR yields the optimal expected outcome when applied to assign treatment to a population of interest. One method for estimating an optimal DTR from observational or randomized study data is Q-learning Schulte et al. 2014). Q-learning is an approximate dynamic programming algorithm that can be viewed as an extension of regression to multi-stage decision problems (Nahum-Shani et al. 2012b). As our focus is the application of Q-learning to the RAD study within the RCP pathway, we focus on data from a two-stage randomized trial with a terminal continuous outcome; however, Q-learning applies in much more general settings (Goldberg and Kosorok 2012;Moodie et al. 2014;Schulte et al. 2014;Sutton and Barto 1998;Watkins and Dayan 1992).
Q-learning estimates an optimal regime using backward induction. For simplicity, we assume that the entire treatment period contains two stages with a distal outcome measured after completion of the second stage; treatment decisions are made in the beginning of each stage. Q-learning proceeds in two steps. In the first step, it estimates an optimal treatment rule for the second stage of treatment given patient-level data accumulated up to and immediately preceding this second treatment assignment. This information includes each patient's baseline information, stage 1 treatment assignment and intermediate, i.e., proximal, outcomes measured during the course of the first stage of treatment. These inputs to the secondstage rule are treated as "independent variables" with no attempt to infer what decision at stage 1 would be optimal for a given patient. This first step is achieved by regressing the distal outcome on patient information up to decision stage 2 and manipulating the obtained analytic Eighty-five patients take Bupropion, 93 patients take Paroxetine, and 187 patients take placebo. After 6 weeks, 104 patients' information are lost. Only 78 patients are tracked with non-response at the end of stage 1. At stage 2, patients with non-response are assigned to secondary treatment intervention. Patients taking Bupropion or Paroxetine at stage 1 will increase current doses. But patients taking placebo at stage 1 will be assigned Bupropion or Paroxetine. expression to find for each patient which treatment at stage 2 optimizes the expected distal outcome.
At the second step, Q-learning looks for treatment assignment at stage 1 that would result in optimal distal outcome, assuming that subsequent stage 2 treatment will be determined by the rule constructed in step 1 of the procedure. Such backward reasoning allows Q-learning to factor in future decisions when making treatment decisions at earlier stages. This can be contrasted with a myopic strategy that only looks at intermediate (proximal) outcomes of a current treatment assignment. For example, treatments at stage 1 may lead to temporary alleviation of symptoms and therefore appear beneficial for a patient; however, the long-term benefits may become questionable after the later (e.g., second) stage decisions are factored in.

Formal mathematical description of Q-learning
We now present formal mathematical description of Qlearning. We assume that data available to estimate a DTR are in the form of n independent, identically distributed , one for each subject where: X 1 ∈ R p 1 denotes baseline (pre-randomization) subject information; A 1 ∈ A 1 denotes the first-stage treatment assignment; X 2 ∈ R p 2 denotes information collected during the course of the first-stage treatment including information dictating first-stage responder status; A 2 ∈ A 2 denotes the second-stage treatment assignment; and Y ∈ R denotes a continuous outcome measured at the end of the study coded so that lower values are better. To match the RAD study, we assume that responders are not re-randomized. In the RAD study, X 1 contains a subject's age, race, gender, marital status, annual household income, employment status, education level, nine different side effect measures, medical insurance type, as well as baseline measures of bipolar type, clinical status prior to depressive episode, scale scores for mood elevation (SUMM), and scale scores for depression (SUMD); A 1 denotes low-dose Bupropion, low-dose Paroxetine, or placebo; X 2 contains responder status at the end of stage 1, as well as SUMM and SUMD at the end of stage 1; A 2 denotes either high-dose Bupropion or highdose Paroxetine; Y is SUMD measured at the end of stage 2.
We use uppercase letters to denote random variables and lowercase letters to denote instances of these ran-dom variables. Define H 1 = X 1 and H 2 = X 1 , A 1 , X 2 , so that H j denotes information available to the decisionmaker at stage j = 1, 2. For any random variable U, let dom U denote the set of possible values for U. A DTR is a pair of functions π = (π 1 , π 2 ) where π j : dom H j → dom A j so that a patient presenting with H j = h j at stage j is assigned treatment π j (h j ). For any h j ∈ dom H j , let F j (h j ) denote the set of feasible treatments for a patient presenting at stage j with H j = h j . In the RAD study F 1 (h 1 ) = {Bupropion, Paroxetine, and placebo}, for all h 1 . At the second stage, responders are not re-randomized; feasible second-stage treatments for non-responders are: denote the class of feasible DTRs (for a more formal discussion of feasibility see (Schulte et al. 2014)). An optimal DTR, say π opt , satisfies E π opt Y ≥ E π Y for all π ∈ , where E π denotes expectation under the restriction that The function Q 2 (h 2 , a 2 ) measures the "quality" of assigning treatment a 2 to a patient presenting at stage 2 with H 2 = h 2 ; the function Q 1 (h 1 , a 1 ) measures the quality of assigning treatment a 1 to a patient presenting at stage 1 with H 1 = h 1 assuming optimal subsequent treatment. It follows from dynamic programming (Bellman 1957) that π opt j (h j ) = arg min a j ∈F j (h j ) Q j (h j , a j ). In practice, dynamic programming cannot be applied because the true Q-functions are not known; instead, estimation of π opt must rely on the observed data. Q-learning is an approximate dynamic programming algorithm which mimics the dynamic programming solution by replacing the conditional expectations required by dynamic programming with regression models fit to the observed data. Let Q j (h j , a j ; θ j ) denote a postulated working model for Q j (h j , a j ) indexed by unknown parameter θ j .
In RAD, only patients who receive placebo as their firststage treatment and failed to respond are randomized at the second stage. Thus, we only use these subjects to estimate θ 2 . Let R denotes a subject first-stage responder status so that R = 1 for responders and R = 0 for nonresponders. Define 1 u to be equal to one if the condition u is true and zero otherwise. A version of the Q-learning algorithm that applies to data from RAD is: Algorithm 1: Q-learning for RAD ; and subsequently estimator Q 2 h 2 , a 2 ; θ 2 of Q 2 (h 2 , a 2 ).
The Q-learning estimator of the optimal regime is π j h j = arg min a j ∈F j (h j ) Q j h j , a j ; θ j . To estimate π opt using data from the RAD study, we posit linear models for the Q-functions. For Q 1 (h 1 , a 1 ), we posit a model of the form Q 1 (h 1 , a 1 ; θ 1 ) = h 10 β 10 + a 11 h 11 β 11 + a 12 h 12 β 12 , where θ 1 = β 10 , β 11 , β 12 , h 1k , k = 0, 1, 2 are known summary vectors of h 1 , and a 1k , k = 1, 2 are dummy variables coding two of the three possible treatments at the first stage. For Q 2 (h 2 , a 2 ), we posit a model of the form Q 2 (h 2 , a 2 ; θ 2 ) = h 20 β 20 + a 2 h 21 β 21 , where θ 2 = β 20 , β 21 , h 2k , k = 0, 1 are known summary vectors of h 2 , and a 2 is a dummy variable coding one of the two possible treatments at the second stage.

Analysis of STEP-BD
The Q-learning algorithm stated in the preceding section assumes (i) complete data and (ii) working models for the Q-functions. However, in RAD, as in most clinical trials, a non-trivial amount of covariate and outcome information are missing. Furthermore, there is no strong theory to suggest working models for the Q-functions. So we must use the data to assist in the choice of these models. We combine multiple imputation with stepwise variable selection to estimate the Q-functions and subsequently the optimal treatment regime. Figure 3 shows the proportions of missing values for the variables under consideration in our analysis of the RAD data. There is a significant amount of missing covariate information at both stages; thus, discarding subjects with missing information is inefficient and may introduce bias (Little and Rubin 2002).

Missing data
One approach in dealing with missing data is multiple imputation (MI) (Rubin 2004). MI creates multiple complete datasets and is thereby suited for conducting a series of exploratory and secondary analyses including estimation of an optimal treatment regime (Shortreed et al. 2011). We use Bayesian MI to "fill in" the missing values which draws from the posterior predictive distribution of the missing values given the observed data (for details and underlying assumptions, see (Little and Rubin 2002;Van Buuren 2012)). Implementation of Bayesian MI requires specification of a prior and likelihood for the observed data. We specify the joint likelihood through the conditional distribution of each variable on all other variables (for discussion of this approach, see (Raghunathan et al. 2001 Buuren 2007)). Thus, the likelihood is determined implicitly through a series of regression models, one for each variable that contains missing information. For continuous variables, we use predictive mean matching, and for binary variables, we use logistic regression models. To reduce variance, we use forward stepwise variable selection applied to the complete data to select predictors for each conditional model. Flat improper priors were used for all parameters. Imputations were carried out using the freely available and open-source mice package with the default settings (http://cran.r-project.org/web/packages/ mice/index.html). Complete R code implementing our imputation model is provided in Additional file 1. Using the procedure described above, we impute m complete datasets. For a given choice of h 1,k , k = 0, 1, 2 and h 2,k , k = 0, 1, we can apply the Q-learning algorithm to each imputed dataset to obtain estimated Qfunctions Q j h j , a j ; θ ( ) j , j = 1, 2, = 1, . . . , m. The final estimated optimal decision rule is obtained as the minimizer of the averaged imputed Q-functions π j h j = arg min a j ∈F j( h j) m −1 m =1 Q j h j , a j ; θ ( ) j .

Estimated optimal treatment regime and empirical results
We use a version of stepwise variable selection to optimize the Bayesian information criteria (BIC); a complete description of this procedure is given in the Appendix section. The variables included in the model for the second-stage Q-function are SIDE3, SUMD1, and SUMM1. The variables included in the model for the first-stage Q-function are AGE, PRONSET, SUMD0, and SUMM0. Thus, the second-stage Q-functions has the form Q 2 (h 2 , a 2 ; θ 2 ) = β 20 h 20 + a 2 β 21 h 21 , where h 20 = (1, SUMM1, SUMD1, SIDE3) , h 21 = (1, SUMM1, SIDE3) , and A 2 is indicator variable for stage 2 treatment coded so that A 2 = 1 denotes high-dose Bupropion and A 2 = 0 denotes high-dose Paroxetine. The estimated coefficients β 20 , β 21 along with 90% bootstrap confidence intervals are shown in Table 2. The table shows that the main effect of A 2 and interaction between second A 2 and SUMM1 is significant at the 90% level. The estimated optimal decision rule is shown in Figure 4. As anticipated by estimated second-stage Q-function, SUMM1 (mood severity) and SIDE3 (sedation side effect) dictate treatment selection; subjects with sedation side effects and low mood severity are recommended to Bupropion, and all others are recommended to Paroxetine. The first-stage Q-function has the form Q 1 (h 1 , a 1 ; θ 1 ) = β 10 h 10 + a 11 β 11 h 11 + a 12 β 12 h 12 , where:  a 11 = 1 if a 1 = Bupriopion otherwise a 11 = 0; a 12 = 1 if a 1 = Paroxetine otherwise a 12 = 0; PRONSET1 = 1 if PRONSET = remission otherwise PRONSET1 = 0; and PRONSET2 = 1 if PRONSET = manic or hypomanic otherwise PRONSET2 = 0. The estimated coefficients and 90% bootstrap intervals (corrected for non-regularity as suggested by ) are listed in Table 3. Figure 5 shows the first-stage optimal decision rule implied by the estimated Q-function. An interesting feature of the first-stage decision rule is that subjects with a (hypo)manic episode immediately preceding the current major depressive episode are recommended to receive placebo. This supports the hypothesis that subjects with (hypo)manic episodes immediately preceding a major depressive episode might not benefit from an adjuvant antidepressant. Figure 5 also shows that among the subjects experiencing remission or mixed/cycling before the Figure 4 Estimated optimal second-stage decision rule. As anticipated by the estimated second-stage Q-function, SUMM1 (scale score for mood elevation) and SIDE3 (sedation side effect) are used to dictate treatment. current major depressive episodes, Bupropion is recommended to older patients and Paroxetine is recommended to younger patients.
Recall that the optimal treatment regime minimizes the expected depression score SUMD measured at week 12. Thus, it is of interest to compare the estimated expected 12-week SUMD under the estimated optimal treatment regime and other potential treatment regimes of interest. Table 4 shows the estimated depression score under the estimated regime and four static treatment regimes. Estimates were computed using the inverse-probabilityweighted estimator (IPWE) (Zhang et al. 2013) and confidence intervals using the non-parametric bootstrap. The estimated optimal regime performs significantly better at the 90% level than any fixed regime under consideration. We note, however, that the confidence intervals must be interpreted with caution as the IPWE is not a smooth estimator which may cause the bootstrap to perform poorly (note that normal approximations do not hold either without strong assumptions ).

Discussion
We estimated an optimal DTR for patients presenting with bipolar depression using data from the RAD pathway in the STEP-BD study. The estimated treatment regime suggests the hypothesis that bipolar-depression patients with (hypo)mania immediately preceding a major depressive episode may do better to forgo adjunctive antidepressant treatment with either paroxetine or bupropion, whereas the opposite is true for who were in remission or experienced a mixed episode (manic episode with mixed Figure 5 Estimated optimal first-stage decision rule. Note that subjects with (hypo)manic episodes immediately preceding the current major depressive episode are recommended to receive placebo.
features, according to DSM-V) before the current major depressive episode. This is a novel finding, which has not been explored so far. At present, there is a consensus that antidepressants in the acute treatment of bipolar depression may be used when there is a history of previous positive response to antidepressants, while they should be avoided in patients with an acute bipolar depressive episode with two or more concomitant core manic symptoms in the presence of psychomotor agitation, in patients with a high number of previous episodes or with a history of rapid cycling and during depressive episodes with mixed features (Pacchiarotti et al. 2013). Furthermore, the use of antidepressants is discouraged if there is a history of past mania, hypomania, or mixed episodes emerging during antidepressant treatment (Pacchiarotti et al. 2013). However, this consensus is mainly based on clinical wisdom than strong external evidence. In our study, the scale scores for measuring symptoms of depression as well as mania were available for baseline and stage 1 to model the Q-functions but did not turn out to be helpful in building an optimal DTR. So far, there are no reliable data on the differential efficacy of paroxetine and bupropion in younger or older adult patients with bipolar depression. In unipolar depression, a recent meta-analysis suggests that the efficacy of antidepressants in general may be reduced in trials involving patients aged 65 years or older (Tedeschini et al. 2011). Similarly, there have not been any reliable data suggesting that patients with higher scores on mood elevation scales do better on paroxetine than bupropion-and vice versa (Pacchiarotti et al. 2013). What is well known on the other hand is that Paroxetine 20 mg/day does not seem to be associated with an increased risk of switch into (hypo)mania in patients with bipolar depression, even in monotherapy (McElroy et al. 2010). The data for our analyses stem from a double-blind, randomized, placebocontrolled trial (Sachs et al. 2007). Consequently, we do not know whether in clinical practice not adding any medication or intervention to a mood stabilizer is of comparable benefit for those who do best on placebo in our analyses (Severus et al. 2012).

Conclusions
As mentioned in the introduction, estimation of an optimal DTR is typically done as a secondary, exploratory analysis and viewed as a method of generating hypotheses for follow-up confirmatory experiments. The latter is just about to start, using patients with bipolar depression being openly treated within the SCP pathway of STEP-BD using the same rating forms, in particular, the clinical monitoring form for mood disorders.

Variable selection
In order to estimate an optimal treatment regime using Q-learning, we need to select which covariates to include the models for the Q-functions. Recall that the 12-week depression SUMD score was used as the outcome (Y ). We identified 24 potential predictors; these predictors are listed in Table 5. We select a subset of these predictors for each Q-function using stepwise variable selection to minimize the Bayes information criterion (BIC) (Schwarz et al. 1978) averaged over the multiply imputed datasets. Let M 2 denote a subset of predictors dictating the features h 2,k , k = 0, 1, and let θ ( ) 2 (M 2 ) denote the coefficients obtained by applying step (Q1) of the Q-learning algorithm with predictors M 2 to the th imputed dataset. Define Y ( ) i (M 2 ), i = 1, . . . , n to be the predicted outcomes computed in step (Q2) of the Q-learning algorithm using the th imputed dataset and parameter θ ( ) 2 (M 2 ). Let M 1 denotes a subset of predictors dictating h 1,k , k = 0, 1, 2, and let θ ( ) 1 (M 1 , M 2 ) denote the coefficients estimated in step (Q3) of the Q-learning algorithm using predictors M 1 and predicted outcomes Y ( ) i (M 2 ), i = 1, . . . , n. In addition, let BIC ( ) (M 2 ) denote the BIC for a second-stage model Q h 2 , a 2 ; θ ( ) 2 (M 2 ) calculated on the th imputed dataset. Similarly, and let BIC ( ) (M 1 , M 2 ) denote the BIC for the first-stage model Q h 1 , a 1 ; θ ( ) 1 (M 1 , M 2 ) calculated on the th imputed dataset. This procedure that we use to construct models for the Q-learning algorithm is: