Policy Learning for Optimal Dynamic Treatment Regimes with Observational Data

Sequential Learning of Optimal Dynamic Treatment Regimes with Observational Data

\setstretch{1.2}

\setstretch{1.3}

Introduction

Many public policies and medical interventions involve dynamics in their treatment assignments. In public policy, for example, some job training programs assign a series of trainings to each individual across multiple stages (e.g., Lechner_2009,Rodriguez_et_al_2022). In clinical medicine, physicians often sequentially assigns medical treatments to each patient depending on the patient's medical history (e.g., Wang_et_al_2012).

\tikzstyle{level 1}=[level distance=3.5cm, sibling distance=3.5cm] \tikzstyle{level 2}=[level distance=3.5cm, sibling distance=2cm]

\tikzstyle{bag} = [text width=4em, text centered] \tikzstyle{end} = [circle, minimum width=3pt,fill, inner sep=0pt]

We consider the setting of sequential treatment assignment (Robins_1986) in which treatments are sequentially assigned to each individual across multiple stages. In this setting, the effect of treatment at each stage is usually heterogenous with respect to the past treatments and associated characteristics. Hence an effective design of sequential treatment assignment should account for such treatment effect heterogeneity at each stage. That is, the optimal decision of treatment to an individual at each stage should depend upon his/her accumulated information at the corresponding stage.

This paper studies statistical decision of sequential treatment assignment using data from an observational study. We assume that the assumption of sequential ignorability (Robins_1997) holds, meaning that the treatment assignment at each stage is independent of potential outcomes conditioning on the history of treatment assignments and observed characteristics. Under this assumption, we construct an approach to learn the optimal Dynamic Treatment Regime (DTR) that is the sequence of the stage-specific policies (treatment rules) and guides the optimal treatment assignment for each individual at each stage based on the individual’s history (Murphy_2003,Chakraborty_Murphy_2014).

We propose two step-wise doubly robust (DR) approaches to learn the optimal DTR. Each approach sequentially learns the optimal stage-specific policy from the final to first stages through backward induction. At each step of the backward induction, each approach constructs an augmented inverse probability weighting (AIPW) estimator of the policy value function for the corresponding stage. The first approach does so by combing estimates of the propensity scores and conditional mean of the potential outcomes as well as the estimated policies for the future stages. The second approach uses the estimates of the state-action value functions (called Q-function (Murphy_2005,Sutton_et_al_2018)) rather than the conditional mean of the potential outcomes. Each approach then estimates the optimal policy for the corresponding stage by maximizing the estimated policy value function over the pre-specified class of stage-specific policies. The optimal DTR is estimated as the sequence of the estimated policies for all stages. Throughout the procedure of each approach, we also use the cross-fitting (chernozhukov_et_al_2018) to make the estimation of the policy value function and the learning of the optimal policy independent. These approaches have the DR property in the sense that the optimal DTR is consistently estimated if at least one of the propensity scores and the conditional mean of the potential outcome or the Q-function are consistently estimated, respectively.

We evaluate the statistical performance of each approach in terms of the regret that is the welfare loss of an estimated DTR relative to the optimal one. For each approach, our main result shows the rate of convergence of the regret of the resulting DTR in relation to rates of convergences of the estimators of the nuisance components and complexity of the class of DTRs. The result also reveals conditions on the class of DTRs and estimators of the nuisance components under which the DTR obtained by each approach achieves the minimax optimal convergence rate $O_{p}(n^{-1/2})$ of the regret. For example, for each approach, if all the nuisance components are estimated with rate $O_{p}(n^{-1/4})$ of convergence in the root-mean-squared error, which is achievable for many machine learning methods under structured assumptions on the nuisance parameters, then the regret of the resulting DTR converges to zero at the rate $O_{p}(n^{-1/2})$. This result is consistent with those of Athey_Wager_2020 and zhou2022offline who study the policy learning in the static case.

Related Work

There is a large literature on the statistical decision of treatment choice, but many works in that literature focus on the static treatment assignment problem; e.g., Manski_2004, Zhao_et_al_2012, Kitagawa_Tetenov_2018a, Athey_Wager_2020, zhou2022offline, and so on. Among these works, this paper is most relevant to Athey_Wager_2020 and zhou2022offline. They propose DR policy learning using the cross-fitting in the static case, and and show that $n^{-1/2}$-upper bound on the regret can be achieved even in the observational data setting. This paper seeks to extend their approaches and results to the dynamic treatment setting.

This work is related to the literature of estimation of the optimal DTRs.\footnote{ Laber_et_al_2014, Chakraborty_Murphy_2014, and Kosorok_et_al_2019 review this literature.} The problem of estimating dynamic sequential decision rules is also called adaptive strategies (Lavori_et_al_2000) or the batch offline policy learning in the reinforcement learning literature (Sutton_et_al_2018).\footnote{Most works in the literature of the batch offline reinforcement learning assume Markov decision process whereas this work does not. Many dynamic treatment cases in social and biomedical sciences are non-Markovian process.} In terms of the statistical setting, this paper is close to classification-based methods of Zhao_et_al_2015 and Sakaguchi_2021 who use propensity weighted outcomes to estimate the value function of a DTR and then estimate the optimal DTR by maximizing it over a pre-specified class of DTRs. Zhao_et_al_2015 develop methods to estimate the optimal DTR using the Support Vector Machine with propensity weighted outcomes, but their focus are on to use experimental data (i.e., the propensity scores are known). Sakaguchi_2021 proposes backward and simultaneous estimating approaches using inverse weights of estimated propensity scores in the observational data setting. These approaches are simple to implement; however, such approach sometimes has problematically high variance and leads to inferior DTRs (see, .e.g., Doroudi_et_al_2018). Our approach improves power of propensity weighting type methods by leveraging outcome regression.

The batch offline Q-learning (Watkins_Dayan_1992) is a dominant approach for estimating the optimal DTRs (e.g., Murphy_2005,Moodie_et_al_2012,Zhang_et_al_2018). This approach sequentially estimates the Q-function (i.e., optimal state-action value function) for each stage, and then estimates the optimal policy from the final to first stage through backward induction. Murphy_2005 shows that the performance of the DTR obtained by the Q-learning depends on how accurately the Q-function is estimated. The resulting DTR may lead to much smaller welfare than the optimal one if the estimated Q-functions are not close to the true ones. Our approaches also use the regression model for the outcome variable. However, leveraging the propensity score models in the AIPW estimations of the policy value functions, our approaches improve the robustness and accuracy of the approaches based on the outcome regression.

DR estimations of the optimal DTRs are also proposed by Zhang_et_2013, Wallace_Moodie_2015, and Nie_et_al_2021. Zhang_et_2013 propose estimating the optimal DTR by maximizing the DR estimator of the welfare function of a DTR over a pre-specified class of DTRs. However, this approach is computationally challenging for two reasons: (i) nuisance components are to be estimated for each specific DTR; (ii) the approach maximizes the estimated welfare function over the whole class of DTRs. Our approach solve these computationally obstacles because (i) nuisance components to be estimated do not depend on any specific DTR and (ii) the optimal DTR is estimated through stage-wise optimization steps rather than simultaneous optimization across all stages. Wallace_Moodie_2015 develops DR estimation of the optimal DTRs building on the Q-learning and G-estimation (Robins_2004). In the problem of learning DTRs for when to start and stop treatment, Nie_et_al_2021 propose DR learning approach with feasible computation property. Jiang_Li_2016 and Thomas_et_al_2016 propose DR evaluation of a fixed DTR, rather than learning of the optimal DTRs.

Structure of the Paper

The remainder of the paper is structured as follows. Section (ref) describes the dynamic treatment framework and defines the dynamic treatment choice problem. Section (ref) presents the two step-wise approaches to learn the optimal DTRs. Section (ref) shows statistical properties of each approach. Section (ref) makes comparison of the proposed approaches to existing ones. Section (ref) presents a simulation study to evaluate their finite sample performances. Section (ref) concludes this paper. Appendix presents proofs of main theorems of this paper and auxiliary lemmas.

Setup

We introduce the dynamic treatment framework, following Robins’s dynamic counterfactual outcomes framework (Robins_1986,Robins_1997), in Section (ref). Subsequently, we define the dynamic treatment choice problem in Section (ref).

Dynamic Treatment Framework

We suppose that there are $T$ $(< \infty)$ stages of multiple treatment assignment. Let $\mathcal{A}_{t} \equiv \{a_{t}^{(1)},\ldots,a_{t}^{(d_t)}\}$ ($t=1,\ldots,T$) denote a set of possible treatment arms in stage $t$, where $d_t$ denotes the number of possible treatment arms in stage $t$ and is allowed to vary across stages. Let $A_{t} \in \mathcal{A}_{t}$ be the treatment in stage $t$. Throughout this paper, for any variable $V_{t}$, we denote by $\text{\ensuremath{\text{\ensuremath{\underline{V}}}}}_{t}\equiv \left(V_{1},\ldots,V_{t}\right)$ the history of the variable up to stage $t$, and denote by $\text{\ensuremath{\underline{V}}}_{s:t} \equiv \left(V_{s},\ldots,V_{t}\right)$, for $s\leq t$, the partial history of the variable from stage $s$ up to stage $t$. Let $Y\left(\text{\ensuremath{\underline{a}}}_{T}\right)$ be the potential outcome that is realized when $\underline{A}_T = \underline{a}_{T}\in \underline{\mathcal{A}}_{T}$, where $\underline{\mathcal{A}}_{T}\equiv \mathcal{A}_{1} \times \cdots \times \mathcal{A}_{T}$. Depending on the history $\underline{A}_T$ of treatments up to the last stage $T$, we observe an outcome $Y \equiv Y(\underline{A}_{T})$ after the end of the sequential treatment assignment. Let $X_{t}$ be a vector of covariates that are observed prior to treatment assignment in stage $t$. The distribution of $X_{t}$ may depend on the past treatments and past covariates.\footnote{In line with Zhang_et_2013 and Nie_et_al_2021, we can define potential covariates $X_{t}(\underline{a}_{t-1})$ that is realized when $\underline{A}_{t-1}=\underline{a}_{t-1}$. Then the observed covariates $X_t$ correspond to $X_{t}(\underline{A}_{t-1})$. We here assume that the covariates are not influenced by the future treatment (i.e., no-anticipation assumption). } We define the history in stage $t$ by $H_{t}\equiv(\underline{A}_{t-1},\underline{X}_{t})$, which is available information for the policy-maker when she chooses a treatment assignment in stage $t$. Let $Z\equiv(\underline{A}_{T},\underline{X}_T,Y)$ be the vector of all the observed variables. We denote the support of $H_{t}$ and $Z$ by ${\cal H}_{t}$ and $\mathcal{Z}$, respectively. Let $P$ denote the distribution of all the defined variables $\left(\underline{A}_T,\underline{X}_T,\{Y(\underline{a}_T)\}_{\underline{a}_T \in \underline{\mathcal{A}}_T}\right)$.

From an observational study, we observe $Z_{i}\equiv \left(\underline{A}_{i,T},\underline{X}_{i,T},T\right)$ for individuals $i=1,\ldots,n$, where $Y_{i} \equiv \sum_{\text{\ensuremath{\underline{a}}}_{T}\in\underline{\mathcal{A}}_{T}}\left(1\left\{ \text{\ensuremath{\underline{A}}}_{i,T}=\text{\ensuremath{\underline{a}}}_{T}\right\}\cdot Y_{i}\left(\text{\ensuremath{\underline{a}}}_{T}\right)\right)$ with $Y_{i}\left(\text{\ensuremath{\underline{a}}}_{T}\right)$ being a potential outcome for individual $i$ that is realized when $\underline{A}_{i,T} = \underline{a}_{T}$. We suppose that the vectors of random variables $\left(\underline{A}_{i,T},\left\{ Y_{i}\left(\text{\ensuremath{\underline{a}}}_{T}\right)\right\} _{\text{\ensuremath{\underline{a}}}_{T}\in\underline{\mathcal{A}}_{T}},\underline{X}_{i,T}\right)$, $i=1,\ldots,n$, are independent and identically distributed (i.i.d) under the distribution $P$. We denote by $P^{n}$ the joint distribution of $\left\{\left(\underline{A}_{i,T},\left\{ Y_{i}\left(\text{\ensuremath{\underline{a}}}_{T}\right)\right\} _{\text{\ensuremath{\underline{a}}}_{T}\in\underline{\mathcal{A}}_{T}},\underline{X}_{i,T}\right):i=1,\ldots,n\right\}$ generated by $P$.

Define $e_{t}\left(h_{t},a_{t}\right)\equiv \Pr\left(A_{t}=a_{t}\mid H_{t}=h_{t}\right)$, the propensity score of the treatment $a_t$ in stage $t$ given the history $h_t$. In the observational data setting we study, the propensity scores are unknown to the analyst and may need to be estimated. This is contrast to the experimental data setting studied by Zhao_et_al_2015 where the propensity scores are known to the analyst.

Throughout the paper, we suppose that the underlying distribution $P$ satisfies the following assumptions.

Assumption (ref) is what is called a dynamic unconfoundedness assumption or sequential ingnorability assumption elsewhere, and is commonly used in the literature of dynamic treatment effect analysis (Robins_1997; Murphy_2003). This assumption means that the treatment assignment at each stage is independent of the contemporaneous and future potential outcomes conditional on the history up to that point. In observational studies, this assumption is satisfied when a sufficient set of confounders is controlled at each stage. Assumptions (ref) and (ref) are standard assumptions in the literature of treatment analysis.

Dynamic Treatment Choice Problem

The aim of this paper is to provide methods to learn the optimal DTRs using data from an observational study. We define a policy at each stage $t$ as $\pi_{t}:\mathcal{H}_{t}\mapsto \mathcal{A}_{t} $, a map from the history space at stage $t$ to the treatment space at stage $t$. A policy $\pi_t$ decides which treatment to be assigned to each individual at stage $t$ based on the history $H_t$. We then define the DTR by $\pi \equiv\left(\pi_{1},\ldots,\pi_{T}\right)$, a sequence of stage-specific policies. The DTR sequentially guides a treatment choice for each individual from the first to last stage depending on the individual’s history up to that point.

We define the welfare of a DTR $\pi$ by

which is the population mean of the outcome that is realized when the sequential treatments are assigned subject to $\pi$. We consider to choose a DTR from a pre-specified class of DTRs, which we denote by $\Pi \equiv \Pi_{1}\times\cdots\times \Pi_{T}$, where $\Pi_{t}$ is a class of policies at stage $t$ (i.e., a pre-specified class of measurable functions $\pi_{t}:\mathcal{H}_{t}\rightarrow \mathcal{A}_{t}$). Zhang_et_al_2018 uses the class of list forms of policies for $\Pi_t$; Laber_Zhao_2015, Tao_et_al_2018, Sun_Wang_2021, and Blumlein_et_al_2022 use the class of trees for $\Pi_t$. These simple policy classes enhance interpretability of the resulting DTR.

The ultimate goal of this paper is to choose an optimal DTR that maximizes the welfare $W\left(\cdot\right)$ over $\Pi$. We are especially interested in learning the optimal DTR from observation data which satisfies the sequential ignorability assumption (Assumption (ref)). The following section introduces sequential approaches to learn the optimal DTR from the observational data.

Sequential Learning of the Optimal DTR

We discuss two sequential approaches to learn the optimal DTRs through backward-induction (dynamic programming). Given a DTR $\pi$ and the class $\Pi$ of DTRs, we denote their partial sequence by $\pi_{s:t}:=(\pi_s,\ldots,\pi_t)$ and $\pi_{s:t}:=\Pi_s \times \cdots \times \Pi_t$, respectively, for $s\leq t$.\footnote{Throughout the paper, for any object $v_{s:t}$ and $\underline{w}_{s:t}$ ($s\leq t$), $v_{t:t}$ and $\underline{w}_{t:t}$ correspond to $v_t$ and $w_t$, respectively.} For each stage $t$, we define the policy value function of $\pi_{t:T}$ as

and the action-value function of $a_t$ for $\pi_{(t+1):T}$ as

We consider to sequentially estimate the optimal DTR by constructing AIPW estimators of the policy value functions and using cross-fitting. In what follows, for any function $f(\cdot,\cdot):\mathcal{H}_{t} \times \mathcal{A}_{t} \rightarrow \mathbb{R}$ and policy $\pi_{t}(\cdot):\mathcal{H}_{t} \rightarrow \mathcal{A}_{t}$, we denote $f(h_{t},\pi_{t}(h_{t}))$ shortly by $f(h_{t},\pi_{t})$ (e.g., $e_{t}(h_t,\pi_t)=e_{t}(h_t,\pi_t(h_t))$).

To sequentially identify and estimate the optimal DTRs through backward-induction, we suppose that the policy class $\Pi_t$ for each stage $t=2,\ldots,T$ contains the first-best policy in the following sense.

We call $\pi_{t}^{\ast,FB}$ that satisfies Assumption (ref) the first-best policy in stage $t$. The first-best policy $\pi_{t}^{\ast,FB}$ always chooses the best treatment arm for any history $h_t$ given that the first-best policies are followed in the future stages. Assumption (ref) is satisfied when $\Pi_t$ ($t=2,\ldots,T$) are flexible enough or correctly specified.

Sakaguchi_2021 illustrates how the backward induction approach does not work to estimate the optimal DTRs when the fist-best policy is not available in the supposed policy classes. Zhang_et_2013 discuss how to correctly specify $\Pi$ so as to satisfy Assumption (ref) depending on models relevant to the treatment effect heterogeneity. Zhao_et_al_2015 use a class of policies represented by a reproducing kernel Hilbert space for each $\Pi_t$, which is flexible enough to approximate the first-best policy. Note that the policy class $\Pi_1$ for the first stage does not need to contain the first-best policy.

We will present two sequential learning approaches of the optimal DTRs in the following two subsections.

First Approach

For each stage $t$, we define $Q_t(\cdot,\cdot):\mathcal{H}_t \times \underline{\mathcal{A}}_{t:T}\rightarrow \mathbb{R}$ as

$Q_t(h_t,\underline{a}_{t:T})$ is the mean of the potential outcome of the sequence of treatments $\underline{a}_{t:T}$ conditional on the history $h_t$ in stage $t$. As we will see, $\{Q_t(\cdot,\cdot)\}_{t=1}^{T}$ and $\{e_{t}(\cdot,\cdot)\}_{t=1}^{T}$ are nuisance components to be estimated in the approach presented in this subsection. Several methods to estimate $Q_{t}(h_t,\underline{a}_{t:T})$ under Assumption (ref) have been proposed (see, e.g., Tsiatis_et_al_2019 (Tsiatis_et_al_2019, Ch. 5.5)). In this paper, we employ the following recursive approach to estimate $\{Q_t(\cdot;\underline{a}_{t:T})\}_{t=1}^{T}$ (with a fixed $\underline{a}_{T}$):

We can apply a flexible regression method (e.g., random forest, lasso, neural network learning) in the regression at each step.

Following the doubly robust policy learning of Athey_Wager_2020 and zhou2022offline, we employ the cross-fitting (Schick_1986,chernozhukov_et_al_2018) to make the estimation of the policy value function and learning of the optimal policy independent; whereby, to reduce the over-fitting. We randomly divide the data set $\{Z_i\}_{i=1,\ldots,n}$ into $K$ evenly-sized folds (e.g., $K=5$). Let $I_k$ be a set of indices of the data in the $k$-th fold and $I_{-k}$ be a set of indices of the data excluded from the $k$-th fold. In what follows, for any statistics $\hat{f}$, we denote by $\hat{f}^{-k}$ the corresponding statistics calculated using data not contained in the $k$-th fold.

The proposed approach consists of multiple steps. As a preliminary step, we estimate the propensity scores $\{e_t(\cdot,\cdot)\}_{t=1,\ldots,T}$ and conditional means of the potential outcome $\{Q_t(\cdot,\cdot)\}_{t=1,\ldots,T}$ by using the data excluded in each fold. For each index $k$ of the fold, we denote by $\hat{e}_{t}^{-k}(\cdot,\cdot)$ and $\widehat{Q}_{t}^{-k}(\cdot,\cdot)$, respectively, the estimators of $e_t(\cdot,\cdot)$ and $Q_t(\cdot,\cdot)$ using data not contained in the $k$-th fold. Flexible regression methods (e.g., random forest, lasso, neural network learning) can be applied to estimate $e_t(\cdot,\cdot)$ and $Q_{t}(\cdot,\cdot)$.

Given the estimators $\{\hat{e}_{t}^{-k}(\cdot,\cdot)\}_{t=1}^{T}$ and $\{\widehat{Q}_{t}^{-k}(\cdot,\cdot)\}_{t=1}^{T}$ for each $k=1,\ldots,K$, the optimal DTR is estimated in the following sequential way. In the first step, regarding the last stage $T$, we make the score function of the treatment $a_T$ in stage $T$ as follows:

where $k(i)$ denotes the fold that contains the $i$-th observation. Note that its sample mean $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T}(a_T)$ is an AIPW estimator of $E\left[Y_{T}\left(\underline{A}_{1:(T-1)},a_{T}\right)\right]$, implying that $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T}(a_T)$ consistently estimates $E\left[Y_{T}\left(\underline{A}_{1:(T-1)},a_{T}\right)\right]$ if either $\widehat{Q}_{T}^{-k(i)}$ or $\hat{e}_{T}^{-k(i)}$ is consistent.

We then find the best candidate policy in stage $T$ by solving

where the objective function $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T}^{\dag}\left(\pi_T(H_{i,T})\right)$ is an AIPW estimator of the policy value function $V_{T}(\pi_T)$.

In the following step, we consider stage $T-1$. Given $\check{\pi}_{T}$, we make the score function of $a_{T-1}$ as

Its sample mean $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T-1}^{\check{\pi}_{T}}(a_{T-1})$ is an AIPW estimator of the action value function $V_{T-1}^{\check{\pi}_T}(a_{T-1})$. We then find the best candidate policy in stage $T-1$ by solving

where the objective function $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T-1}^{\check{\pi}_{T}}\left(\pi_{T-1}(H_{i,T-1})\right)$ is an AIPW estimator of the policy value $V_{T-1}(\pi_{T-1},\check{\pi}_{T})$.

Recursively, for $t=T-2,\ldots,1$, we learn the optimal policy $\pi_{t}^{\ast}$ as follows. We first make the score function of $a_t$ as

We then find the best candidate policy in stage $t$ by solving

where the objective function $(1/n)\sum_{i=1}^{n} \widehat{\Gamma}_{i,t}^{\check{\pi}_{(t+1):T}}\left(\pi_{t}(H_{i,t})\right)$ is an AIPW estimator of the policy value $V_{t}\left(\pi_{t},\check{\pi}_{t+1},\ldots,\check{\pi}_{T}\right)$.

Throughout this procedure, we eventually obtain the sequence $\check{\pi} \equiv (\check{\pi}_{1},\ldots,\check{\pi}_{T})$, which is the resulting estimator of the DTR. We will clarify statistical properties of $\check{\pi}$ with respect to its welfare regret in Section (ref). The following remark gives some discussion about the definition of the score function ((ref)).

Second Approach

The second approach we propose is also based on the backward induction and similar to the first approach. Given a DTR $\pi$, for stage $t=1,\ldots,T-1$, we define the Q-function for the sequence of the policies $\pi_{(t+1):T}$ as follows:

which is the conditional mean of the potential outcome at stage $t$ when the pair of history and treatment is $(h_t,a_T)$ and the future treatments are subject to $\pi_{(t+1):T}$. For $t=T$, we define $Q_{T}(h_T,a_T)\equiv E\left[Y\left(\underline{a}_{T}\right) \middle| H_T = h_T\right]$, which is equal to $Q_{T}(h_T,a_T)$.

The approach presented here also consists of multiple steps. As a preliminary step, we estimate the propensity scores $\{e_t(\cdot,\cdot)\}_{t=1,\ldots,T}$ and Q-function $Q_T(\cdot,\cdot)$ for the last stage by using the data excluded in each fold. For each index $k$ of the fold, we denote by $\hat{e}_{t}^{-k}(\cdot,\cdot)$ and $\widehat{Q}_{T}^{-k}(\cdot,\cdot)$, respectively, the estimators of $e_t(\cdot,\cdot)$ and $Q_T(\cdot,\cdot)$ using data not contained in the $k$-th fold. Any regression method can be used to estimate $e_t(\cdot,\cdot)$ and $Q_{T}(\cdot,\cdot)$.

Given $\{\hat{e}_{t}^{-k}(\cdot,\cdot)\}_{t=1}^{T}$ and $\widehat{Q}_{T}^{-k}(\cdot,\cdot)$ for each $k=1,\ldots,K$, the optimal DTR is estimated in the following sequentially way. In the first step, regarding the last stage $T$, we make the score function of the treatment $a_T$ for stage $T$ as follows:

Note that its sample mean $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T}^{Q}(a_T)$ is an AIPW estimator of $E\left[Y_{T}\left(\underline{A}_{1:(T-1)},a_{T}\right)\right]$.

We then find the best candidate policy in stage $T$ by solving

where the objective function $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T}^{Q}\left(\pi_T(H_{i,T})\right)$ is an AIPW estimator of the policy value function $V_{T}(\pi_T)$. Note that this step is not different from the first of the approach presented in Section (ref) because $\widehat{Q}_{T}^{-k}(\cdot,\cdot)$ does not substantially　differ from $\widehat{Q}_{T}^{-k}(\cdot,\cdot)$.

In the following step, we consider stage $T-1$. Given $\check{\pi}_{T}$, for each index $k$ of the cross-fitting fold, we estimate $Q_{T-1}^{\check{\pi}_{T}}(\cdot,\cdot)$ by regressing $\widehat{Q}_{T}^{-k}(H_{i,t},\check{\pi}_{T})$ on $(H_{i,T-1},A_{i,T-1})$ using the observations whose indices are not contained in $I_k$.\footnote{This step is similar to the regression step in the batch offline Q-learning (Murphy_2005). The difference is that the Q-learning estimates the optimal policy using the regression equation of the outcome whereas our approach estimates the optimal policy by maximizing the AIPW estimator of the policy value function.} We denote by $\widehat{Q}_{T-1}^{\check{\pi}_{T},-k}(\cdot,\cdot)$ the resulting estimator of $Q_{T-1}^{\check{\pi}_{T}}(\cdot,\cdot)$ for each fold $k$. Any regression method can be used in this step.

We next make the score function of $a_{T-1}$ as

Its sample mean $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T-1}^{Q,\check{\pi}_{T}}(a_{T-1})$ is an AIPW estimator of the action value function $V_{T-1}^{\check{\pi}_{T}}(a_{T-1})$. We then find the best candidate policy in stage $T-1$ by solving

where the objective function $(1/n)\sum_{i=1}^{n}\widehat{\Gamma}_{i,T-1}^{Q,\check{\pi}_{T}}\left(\pi_{T-1}(H_{i,T-1})\right)$ is an AIPW estimator of the policy value $V_t(\pi_{T-1},\check{\pi}_T)$.

Recursively, for $t=T-2,\ldots,1$, we learn the optimal policy as follows. For each cross-fitting index $k$, we first estimate $Q_{t}^{\check{\pi}_{(t+1):T}}$ by regressing $\widehat{Q}_{t+1}^{\check{\pi}_{(t+2):T}}(H_{i,t+1},A_{i,t+1})$ on $(H_{i,t},A_{i,t})$ using the observations whose indices are not in $I_k$. Any regression method can be used. Define the score function as

We find the best candidate policy in stage $t$ by solving

Eventually, we obtain the resulting DTR $\check{\pi}=(\check{\pi}_{1},\ldots,\check{\pi}_T)$. The following section clarifies the statistical property of $\check{\pi}$.

Statistical Properties

Given a DTR $\pi \in \Pi$, we define the regret of $\pi$ by $R(\pi) \equiv \sup_{\tilde{\pi} \in \Pi}W(\tilde{\pi}) - W(\pi)$, the loss of the welfare of $\pi$ relative to the maximum welfare achievable in $\Pi$. We study statistical properties of $\check{\pi}$ and $\check{\pi}$ with respect to their regrets $R(\check{\pi})$ and $R(\check{\pi})$, respectively. This section shows the doubly-robust property of the the two approaches and derives the rate of convergences of $R(\check{\pi})$ and $R(\check{\pi})$ depending on the rates of convergence of the estimators of the nuisance components and complexity of $\Pi$.

We first focus on the DTR $\check{\pi}$ obtained by the approach in Section (ref). Let $\widehat{Q}_{t}^{(n)}(\cdot,\cdot)$ and $\hat{e}_{t}^{(n)}(\cdot,\cdot)$, respectively, denote estimators of the conditional mean of the potential outcome $Q_{t}(\cdot,\cdot)$ and propensity score $e_t(\cdot,\cdot)$ using size $n$ sample randomly drawn from the population $P$. We suppose that $\{\widehat{Q}_{t}^{(n)}(\cdot,\cdot)\}_{t=1}^{T}$ and $\{\hat{e}_{t}^{(n)}(\cdot,\cdot)\}_{t=1}^{T}$ satisfy the following assumption.

As we will see later, the $\sqrt{n}$-consistency of the regret $R(\check{\pi})$ to zero can be achieved when Assumption (ref) (i) holds with $\tilde{\tau}_{B}=1$. This is not very strong or restrictive. For example, Assumption (ref) (i) is satisfied with $\tilde{\tau}_B = 1$ when

hold for all $t=1,\ldots,T$ and $s=1,\ldots,t$. Note also that Assumption (ref) (i) encompasses the doubly-robustness property; that is, Assumption (ref) (i) is satisfied if at least one of $\widehat{Q}_{t}^{(n)}(\cdot,\cdot)$ and $\prod_{s=t}^{T}\hat{e}_{s}^{(n)}(\cdot,\cdot)$ is consistent in terms of the mean squared error.

We next consider complexity of the class $\Pi$ of DTRs and each policy class $\Pi_t$ ($t=1,\ldots,T$). Following zhou2022offline, we use the $\epsilon$-Hamming covering number to measure the complexity of the sequence of policies $\pi_{t:T} \in \Pi_{s}\times \cdots \times \Pi_{t}$ for each $s \leq t$.

When $s=t$, note that $N_{d_h}(\epsilon,\Pi_t) = N_{d_h}(\epsilon,\Pi_{t:t})$ and $\kappa(\Pi_t)=\kappa(\Pi_{t:t})$. We suppose that $\Pi_t$ for each $t$ is not too complex in terms of the covering number.

This assumption implies that the covering number of $\Pi_t$ does not grow too quickly, but allows that $\log N_{d_h}(\epsilon,\Pi_{t})$ grows at a rate of $1/\epsilon$. The assumption is satisfied, for example, by the policy class of finite-depth trees (see zhou2022offline (zhou2022offline, Lemma 4)). zhou2022offline (zhou2022offline, Remark 4) shows that the entropy integral $\kappa(\Pi_t)$ is finite under Assumption (ref). As for the class $\Pi$ of whole DTRs, $\kappa(\Pi)$ is finite as well under Assumption (ref).

The following theorem is a main result of this paper and show the rate of convergence for the regret of the approach proposed in Section (ref).

When Assumption (ref) (i) holds with $\tilde{\tau}_B = 1$, the approach proposed in Section (ref) achieves the minimax optimal rate $O_ {p}(1/\sqrt{n})$ of convergence of the regret with respect to the sample size $n$.\footnote{Sakaguchi_2021 shows that the minimax optimal convergence rate of the regret is $O_ {p}(1/\sqrt{n})$ in the case that treatment assignment at each stage is binary.} This result is consistent with those of Athey_Wager_2020 and zhou2022offline who study static policy learning and that of Nie_et_al_2021 who study learning problem of when-to-treat policy.

We next consider the regret property of the DTR $\check{\pi}$ obtained by the approach presented in Section (ref). Let $\widehat{Q}_{t}^{\pi_{(t+1):T},(n)}(\cdot,\cdot)$ denote the estimator of the Q-function $Q_{t}^{\pi_{(t+1):T}}$ for the sequence of fixed policies $\pi_{(t+1):T} \in \Pi_{(t+1):T}$ using size $n$ sample randomly drawn from the population $P$. $\widehat{Q}_{t}^{\pi_{(t+1):T},(n)}(\cdot,\cdot)$ is obtained by the sequential regression procedure in Section (ref) with $\check{\pi}_s$ replaced by $\pi_s$ for $s=t+1,\ldots,T$. We suppose that $\left\{\widehat{Q}_{t}^{\pi_{(t+1):T},(n)}(\cdot,\cdot):\pi_{(t+1):T} \in \Pi_{(t+1):T}\right\}_{t=1,\ldots,T-1}$, $\widehat{Q}_{T}^{(n)}(\cdot,\cdot)$, and $\{\hat{e}_{t}(\cdot,\cdot)\}_{t=1,\ldots,T}$ satisfy the following assumption.

The result for convergence rate of the estimated Q-function $\widehat{Q}_{t}^{\pi_{(t+1):T},(n)}(\cdot,\cdot)$ is not a standard result because $Q_{t}^{\pi_{(t+1):T}}(\cdot,\cdot)$ is not smooth due to the existence of the policies $\pi_{(t+1):T}$ in its definition. However, some results might be applicable. For example, Zhang_et_al_2018 derive the rate of convergence result for the batch offline Q-learning using Support Vector Machine regression with the policy search over the class of list forms of policies.

The following theorem shows the statistical property of the DTR $\check{\pi}$ in terms of the rate of the convergence of the regret.

The theorem implies that the minimax optimal rate of the regret, which is $O_p(n^{-1/2})$, can be achieved by $\check{\pi}$ when Assumption (ref) holds with $\tilde{\tilde{\tau}}_{B}=1$. This condition holds, for example, when

hold for all $t=1,\ldots,T-1$, $t^\prime=1,\ldots,T$, and $s=1,\ldots,t^\prime$. Note also that Assumption (ref) (i) encompasses DR property, meaning that Assumption (ref) (i) holds if at least one of $\widehat{Q}_{t}^{\pi_{(t+1):T}}(\cdot,\cdot)$ and $\hat{e}_t(\cdot,\cdot)$ is consistent in terms of the mean squared error.

Comparison to Existing Approaches

We make comparison of the proposed approaches to several existing approaches to lean the optimal DTRs with observational data. One approach to learn the optimal DTRs is to use the inverse probability weighting (IPW) to estimate the welfare function $W(\pi)$ or the policy value function $V_t(\pi_{t:T})$ for each stage $t$, and then maximize the estimated welfare function of policy valued function over $\Pi$. Zhao_et_al_2015 and Sakaguchi_2021 employ this approach. In the sequential learning of this approach, we can learn the optimal DTR by $\check{\pi}^{B} = (\check{\pi}_{1}^{B},\ldots,\check{\pi}_{T}^{B})$ which recursively solves

for $t=T-1,\ldots,1$. A drawback of this approach is that the estimated DTR is inconsistent when the propensity scores are inconsistently estimated. One way to avoid this is applying nonparametric estimation methods to estimate the propensity scores. In this case, the rate of convergence of the regret of the resulting DTR, which depends on that of the estimated propensity scores, would be slow. The two approaches proposed in this paper improves this rate of convergence by leveraging the models relevant to the outcome variable.

An alternative approach is the batch offline Q-learning (Murphy_2005,Sutton_et_al_2018). This approach is also not robust to the missspecification of the optimal state-action value function called optimal Q-function. Nonparametric regression methods can be applied to the Q-learning; however, the convergence rate of the regret of the resulting DTR depends on that of the estimated Q-function and would be slow (see Murphy_2005). Our approach improves the convergence property of the Q-learning by leveraging the AIPW estimators of the policy value function.

Another approach is estimating the welfare function $W(\pi)$ of a whole DTR $\pi$ by using an AIPW estimator as follows (Zhang_et_2013, Jiang_Li_2016, Thomas_et_al_2016):

where $\hat{\gamma}_{i,t}(\pi)$ and $\hat{q}_{t}^{\pi}(h_{t})$ are estimators of $\gamma_{i,t}(\pi)=\left(\prod_{s=1}^{t}1\left\{ A_{i,s}=\pi_{s}(H_{i,s})\right\}\right)/\left( \prod_{s=1}^{t}e_{t}\left(H_{i,s},\pi_{s}\right)\right)$ and $q_{t}^{\pi}\left(h_{t}\right)=E\left[Y\left(\underline{A}_{t-1},\pi_{t}(H_{t}),\ldots,\pi_{T}(H_{T})\right)\middle| H_{t}=h_{t}\right]$, respectively. The estimator ((ref)) generalizes the AIPW of Robins_et_al_1994 beyond the static case, and is consistent to $W(\pi)$ if either the propensity weights $\{\hat{\gamma}_t(\cdot)\}_{t=1}^{T}$ or the conditional value estimators $\{\widehat{q}_{t}^{\pi}(\cdot)\}_{t=1}^{T}$ are consistent.

Having $\widehat{W}^{DR}(\pi)$, we can find the candidate of the optimal DTRs as a solution of the following maximization problem\footnote{Jiang_Li_2016 and Thomas_et_al_2016 study the evaluation of a fixed DTR $\pi$ rather than estimation of the optimal DTR.}:

However, this maximization problem is computationally challenging due to the following two reasons. First, the nuisance components $\gamma_{i,t}(\pi)$ and $q_{t}^{\pi}(\cdot)$ depend on the DTR $\pi$, requiring to estimate $\gamma_{i,t}(\pi)$ and $q_{t}^{\pi}(\cdot)$ for each specific $\pi$. Second, the problem ((ref)) maximizes $\widehat{W}^{DR}(\pi)$ simultaneously over the whole class $\Pi$ of DTRs rather than sequentially maximizing it. The simultaneous optimization is computationally demanding when $\Pi$ is not a very simple class or $T$ is not very small. Our learning approaches solves these optimization issues. First, as for the approach presented in Section (ref), the nuisance components $\{e_t(\cdot,\cdot)\}_{t=1}^{T}$ and $\{Q_t(\cdot,\cdot)\}_{t=1}^{T}$ do not depend on any specific policy, making the estimation of the nuisance components independent from the searching for the optimal policy. As for the approach presented in Section (ref), the nuisance components $\{Q_{t}^{\check{\pi}_{(t+1):T}}(\cdot,\cdot)\}_{t=1}^{T-1}$ depend only on the sequence of the estimated policies $\check{\pi}_{(t+1):T}$; hence, the computation is much easier than estimating $Q_{t}^{\pi_{(t+1):T}}(\cdot,\cdot)$ for each $\pi_{(t+1):T}$. Second, the proposed approaches sequentially solve the optimization problem from the last to first stage, which is computationally easier than optimizing the estimated welfare function over a whole DTR.

Simulation Study

We conduct a simulation study to examine the finite sample performance of the proposed learning approach. We consider data generating processes (DGPs) that consist of two stages of binary treatment assignment $(A_{1},A_{2}) \in \{0,1\}^2$, associated potential outcomes $\left\{Y_{1}\left(a_{1},a_{2}\right)\right\}_{\left\{ a_{1},a_{2}\right\} \in\left\{ 0,1\right\} ^{2}}$, and $p$ covariates $(X_{1}^{(1)},\ldots,X_{1}^{(p)})$ observed at the first stage, where $p=5$ or $20$ and one covariate $X_2$ observed at the second stage. The data is generated as follows:

The treatment $a_1$ in the first stage has influence on the outcome through the multiple paths.

We compare the performance of the proposed learning approach proposed (labeled “DR”), the IPW-based approach (labeled “IPW”) defined as ((ref)) and ((ref)), and Q-learning (labeled “Q-learn”). For each method, we use the random forest of Athey_et_al_2019 to estimate the nuisance components. When $p$ is large, the estimations of the nuisance components become noisy.

For the proposed approach and IPW-based approach, we use the following class of DTRs: $\Pi=\Pi_{1} \times \Pi_{2}$ with

This class of DTRs contains the first-best policy of the second stage, satisfying Assumption (ref). For each of the proposed and IPW-based approach, the optimization problem at each step can be formulated as Mixed Integer Linear Programming (see Sakaguchi_2021(Sakaguchi_2021, Appendix B)) (MILP) problems, for which some efficient softwares (e.g., CPLEX; Gurobi) are available . We set $K=5$ for the proposed approach.

Figure \textcolor{blue}{1} shows the results of 200 simulations with sample sizes $n=$ 200, 500, 800, 1100, and 1400, where we calculate the mean welfare achieved by each estimated DTR with 3,000 observations randomly drawn from the same DGP used in the estimation. The results show that the B-DR performs better than B-DEWM and Q-learning in terms of the mean welfare especially when the sample size is not large.

Conclusion

We studied statistical learning of the optimal DTRs. We proposed two doubly-robust learning approaches using observational data under the assumption of sequential ignorability. The first approach solves the treatment assignment problem at each stage through backward induction, where we combine the batch Q-learning and inverse propensity score estimators to construct the doubly-robust estimator of the treatment score at each stage. The second approach solves the whole dynamic treatment assignment problem simultaneously across all stages, where we use the doubly-robust estimator of the welfare function. Our main results show that based on the doubly-robust estimators of the treatment scores or welfare function and using cross-fitting, each of the approaches has the property of the doubly-robust learning of the optimal DTRs. Furthermore, under relatively mild conditions on the MSE convergence rate for the estimators of the nuisance components, we show that the proposed approaches achieve the minimax optimal convergence rate $O_ {p}(1/\sqrt{n})$ of welfare regret even if nuisance components are non-parametrically estimated. The simulation study confirms this theoretical property in finite sample setting.