Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
133,606 characters · 19 sections · 79 citation commands
Estimating and Improving Dynamic Treatment Regimes With a Time-Varying Instrumental Variable
\sectionfont
\subsectionfont
{\bf Abstract}: Estimating dynamic treatment regimes (DTRs) from retrospective observational data is challenging as some degree of unmeasured confounding is often expected. In this work, we develop a framework of estimating properly defined “optimal” DTRs with a time-varying instrumental variable (IV) when unmeasured covariates confound the treatment and outcome, rendering the potential outcome distributions only partially identified. We derive a novel Bellman equation under partial identification, use it to define a generic class of estimands (termed IV-optimal DTRs), and study the associated estimation problem. We then extend the IV-optimality framework to tackle the policy improvement problem, delivering IV-improved DTRs that are guaranteed to perform no worse and potentially better than a pre-specified baseline DTR. Importantly, our IV-improvement framework opens up the possibility of strictly improving upon DTRs that are optimal under the no unmeasured confounding assumption (NUCA). We demonstrate via extensive simulations the superior performance of IV-optimal and IV-improved DTRs over the DTRs that are optimal only under the NUCA. In a real data example, we embed retrospective observational registry data into a natural, two-stage experiment with noncompliance using a time-varying IV and estimate useful IV-optimal DTRs that assign mothers to high-level or low-level neonatal intensive care unit based on their prognostic variables.
{\bf Keywords}: Causal inference, Dynamic treatment regime, Instrumental variable, Offline reinforcement learning, Retrospective observational data
\addtocontents{toc}{\setcounter{tocdepth}{0}}
Estimating single-stage individualized treatment rules (ITRs) and the more general multiple-stage dynamic treatment regimes (DTRs) has attracted a lot of interest from diverse disciplines. Estimating optimal policies (ITRs or DTRs) can be challenging when data come from retrospective observational databases where some degree of unmeasured confounding is often expected. In these scenarios, an instrumental variable (IV) is a useful tool to infer the treatment effect. Motivated by recent works on estimating optimal ITRs using an instrumental variable cui2019semiparametric, qiu2020optimal, pu2020estimating and literature on policy improvement (kallus2018interval,kallus2020confounding,kallus2020minimax), we study in this article how to leverage information contained in a time-varying instrumental variable to estimate properly-defined “optimal” DTRs and improve upon pre-specified baseline DTRs.
An instrumental variable is valid if it is associated with the treatment, affects the outcome only through its association with the treatment, and is independent of the unobserved treatment-outcome confounding variables, possibly conditional on a rich set of observed covariates. One subtlety in IV-based analysis lies in that even a valid IV cannot always identify the mean potential outcome; rather, a valid IV along with appropriate, application-driven IV identification assumptions places certain restrictions on the potential outcome distributions. This line of research is known as partial identification of probability distributions (manski2003partial). This subtlety is inherited by the policy estimation problem with an IV. In particular, when the conditional average treatment effect (CATE) is not point identified from data, the optimal policy that maximizes the value function cannot be identified either, necessitating researchers to target alternative optimality criteria. While such criteria have been proposed in the single-stage setting from different perspectives (see, e.g., cui2019semiparametric,cui2021machine, pu2020estimating), the literature on the more complicated, multiple-stage setting is scarce.
Our first primary interest in this article is to extend optimality criteria in single-stage settings (murphy2003optimal,cui2019semiparametric,cui2021machine,pu2020estimating) and develop an optimality criterion that is tailored to general sequential decision problems and incorporates the rich information contained in a time-varying IV. Our optimality criterion, termed IV-optimality, is based on a carefully weighted version of the partially identified $Q$-function and value function subject to the distributional constraints imposed by the IV. This criterion is distinct from the framework of han2019optimal, who endows the collection of partially identified DTRs with a partial order and characterizes the set of maximal elements (see zhang2020selecting for similar ideas). It also distinct from the recent work in the reinforcement learning literature liao2021instrumental that directly models the transition dynamics and uses instrumental variables to identify the relevant causal parameters. In particular, we do not impose Markovian assumptions and do not pose parametric models a priori. We then take a hybrid approach of $Q$-learning watkins1992q,schulte2014q and weighted classification zhang2012estimating,zhao2012estimating to target this optimality criterion and establish non-asymptotic rate of convergence of the proposed estimators.
The IV-optimality framework also motivates a conceptually simple yet highly informative variant framework that allows researchers to leverage a time-varying IV to improve upon a baseline DTR. The policy improvement problem was first considered in a series of papers by kallus2018interval,kallus2020confounding,kallus2020minimax under a “Rosenbaum-bounds-type” sensitivity analysis model rosenbaum2002observational. Despite its novelty and usefulness in a range of application scenarios, a sensitivity-analysis-based policy improvement framework does have a few limitations. First, each improved policy is indexed by a sensitivity parameter $\Gamma$ that controls the degree of unmeasured confounding. Since the sensitivity parameter is not identified from the observed data, it is often unclear which improved policy best serves the purpose. Second, as pointed out by heng2020interactions, “Rosenbaum-bounds-type” sensitivity analysis model does not fully take into account unmeasured confounding heterogeneity; see also bonvini2019sensitivity. More importantly, one major limitation of a sensitivity-analysis-based policy improvement framework is that it cannot improve upon the NUCA-optimal policy, i.e. the policy that is optimal under the no unmeasured confounding assumption rosenbaum1983central, Robins1992. Intuitively, this is because a sensitivity analysis model and a fixed sensitivity parameter only introduce a partial order (rather than a total order) among all candidate policies and the NUCA-optimal policy always remains a maximal element in this partial order zhang2020selecting. Therefore, their framework provably cannot improve upon the NUCA-optimal policy, a policy of major interest in many application scenarios. See Section (ref) for a detailed discussion. As we will demonstrate in this article, an IV-based policy improvement framework solves all aforementioned problems simultaneously.
The rest of the article is organized as follows. We describe a real data application in Section (ref). Section (ref) reviews alternative optimality criteria in the single-stage setting and provides a general IV-optimality framework that is amenable to being extended to the multiple-stage setting. Section (ref) considers improving upon a baseline ITR with an IV. Building upon the preparations in Sections (ref) and (ref), we describes the IV-optimality framework for policy estimation in the multiple-stage setting in Section (ref) and extends this framework to tackle the policy improvement problem in Section (ref). Section (ref) studies the theoretical properties of the proposed methods. We conducted extensive simulations in Section (ref) and revisited the application in Section (ref). Section (ref) concludes with a brief discussion. For brevity, all proofs are deferred to the Supplementary Material.
lorch2012differential constructed a retrospective cohort study to investigate the effect of delivery hospital on premature babies (gestational age between $23$ and $37$ weeks) and found a significant benefit to neonatal outcomes when premature babies were delivered at hospitals with high-level neonatal intensive care units (NICU) compared to those without NICUs. lorch2012differential used the differential travel time to the nearest high-level versus low-level NICU as an IV so that the outcome analysis is less confounded by mothers' self-selection into high-level NICUs. Put another way, the differential travel time creates a natural experiment with noncompliance: mothers who live relatively close to a high-level NICU were encouraged to deliver, although not necessarily delivered, at a high-level NICU. More recently, michael2020instrumental considered mothers who delivered exactly two babies from 1996 to 2005 in Pennsylvania, and investigated the cumulative effect of delivering at high-level NICUs on neonatal survival probability using the same differential travel time IV.
Currently, there is still limited capacity at high-level NICUs, so it is not practical to direct all mothers to these high-technology, high-volume hospitals. Understanding which mothers would most significantly benefit from delivering at a high-level NICU helps design optimal perinatal regionalization systems that designate hospitals by the scope of perinatal service provided and designate where infants are born or transferred according to the level of care they need at birth (lasswell2010perinatal; kroelinger2018comparison). Indeed, previous studies seemed to suggest that although high-level NICUs significantly reduced deaths for babies of small gestational age, they made little difference for almost mature babies like 37 weeks (yang2014estimation).
Mothers who happened to relocate during two consecutive deliveries of babies constitute a natural two-stage, two-arm randomized controlled trial with noncompliance and present a unique chance to investigate an optimal dynamic treatment regime. See Figure (ref) for the directed acyclic graph (DAG) illustrating this application. We will revisit the application after developing theory and methodology precisely suited for this purpose.
We first define the estimand of interest, an optimal ITR, under the potential outcome framework (neyman1923application, rubin1974estimating), and briefly review how to estimate the optimal ITR under the no unmeasured confounding assumption.
Consider a single-stage decision problem where one observes the covariates $\bsX = \bsx \in \calX \subseteq\bbR^{d}$, takes a binary action $a\in\{\pm 1\}$, and receives the response $Y (a)\in \calY \subseteq\bbR$. Here, $Y(a)$ denotes the potential outcome under action $a$. This decision-making process is formalized by the notion of individualized treatment rule (or a single-stage policy), which is a map $\pi$ from available prognostic variables $\bsx$ to a treatment decision $a = \pi(\bsx)$. Denote by $p^\star_a(\cdot|\bsx)$ the conditional law of $Y(a)$ given the realization of $\bsX = \bsx$. The quality of $\pi$ is quantified by its value:
where the outer expectation is taken with respect to the law of the covariates $\bsX$ and the inner expectation with respect to the potential outcome distribution $p^\star_{a}(\cdot|\bsX)$ with $a$ set to $\pi(\bsX)$. Intuitively, $V^\pi$ measures the expected value of the response if the population were to follow $\pi$. Given a class of candidate ITRs $\Pi$, which we refer to as a policy class, let $V^\star = \max_{\pi \in \Pi} V^\pi$ denote the maximal value of an ITR when restricted to $\Pi$. An ITR is said to be optimal with respect to $\Pi$ if it achieves $V^\star$.
A well-known result of zhang2012estimating (see also zhao2012estimating) asserts the {duality} between value maximization and risk minimization, in the sense that any optimal ITR that maximizes the value $V^\pi$ also minimizes the following risk (and vice versa):
where $\mathscr{C}^\star(\bsx) = \bbE_{Y\sim p^\star_{+1}(\cdot| \bsx)}\left[Y\right] - \bbE_{Y\sim p^\star_{-1}(\cdot|\bsx)}\left[Y\right]$ is the CATE.
It is not hard to show that the sign of the CATE, $\sgn(\mathscr{C}^\star(\bsx))$, is the Bayes ITR (i.e., the optimal ITR when $\Pi$ consists of all Boolean functions). Thus, the risk $\calR^\pi$ admits a natural interpretation as a weighted misclassification error: if the decision made by $\pi$ disagrees with the Bayes ITR, then $|\mathscr{C}^\star(\bsx)|$-many units of loss are incurred.
Suppose we have a dataset consisting of i.i.d. samples from the law of the triplet $(\bsX^\textnormal{\texttt{obs}}, A^{\textnormal{\texttt{obs}}}, Y^\textnormal{\texttt{obs}})$. In parallel to the potential outcome distribution $p^\star_{a}(\cdot|\bsx)$, let $p^\textnormal{\texttt{obs}}_a(\cdot|\bsx)$ denote the conditional law of $\{Y^\textnormal{\texttt{obs}}|A^\textnormal{\texttt{obs}}=a, \bsX^\textnormal{\texttt{obs}} = \bsx\}$. One can then define counterparts of the value and the risk as in (ref) and (ref) respectively, but with $p^\star_a(\cdot|\bsx)$ replaced by $p^\textnormal{\texttt{obs}}_a(\cdot|\bsx)$. For instance, we can define
where $\mathscr{C}^\textnormal{\texttt{obs}}(\bsx) = \bbE_{Y\sim p^\textnormal{\texttt{obs}}_{+1}(\cdot| \bsx)}\left[Y\right] - \bbE_{Y\sim p^\textnormal{\texttt{obs}}_{-1}(\cdot|\bsx)}\left[Y\right]$. Unlike $p^\star_{a}(\cdot|\bsx)$, which is defined on the potential outcomes, the distribution $p^{\textnormal{\texttt{obs}}}_{a}(\cdot|\bsx)$ and hence $\calR^{\textnormal{\texttt{obs}}, \pi}$ are always identified from the observed data, thus rendering the task of minimizing $\calR^{\textnormal{\texttt{obs}}, \pi}$ feasible using the observed data. Under a version of the no unmeasured confounding assumption, the two distributions $p^\star_a(\cdot|\bsX)$ and $p^\textnormal{\texttt{obs}}_a(\cdot|\bsX)$ agree, and thus a minimizer of $\calR^{\textnormal{\texttt{obs}}, \pi}$ is indeed an optimal ITR in that it also minimizes $\calR^\pi$. To make the distinction clear, we refer to any minimizer of $\calR^{\textnormal{\texttt{obs}}, \pi}$ as a {NUCA-optimal} ITR (i.e., it is only optimal in the conventional sense under the NUCA) and denote it as $\pi_\texttt{NUCA}$. Many estimation strategies targeting the NUCA-optimal ITR have been proposed in the literature; see, e.g., structural equation models and its variants (murphy2001marginal,murphy2003optimal), outcome weighted learning and its variants (zhao2012estimating, zhao2019efficient, athey2020policy), tree-based methods (laber2015tree, zhang2018interpretable), among others.
The no unmeasured confounding assumption is often a heroic assumption when data come from retrospective observational studies and should be made with caution. Indeed, pu2020estimating demonstrates via extensive simulations that a NUCA-optimal ITR could have poor generalization performance when the NUCA fails.
In classical causal inference literature, an instrumental variable is a widely-used tool to infer the causal effect from noisy observational data (AIR1996, rosenbaum2002covariance, Imbens2004, hernan2006instruments). A random variable $Z$ is said to be a valid IV if it satisfies the {core IV assumptions}: Stable Unit Treatment Value Assumption (SUTVA), correlation between IV and treatment, exclusion restriction (ER), and IV unconfoundedness conditional on the observed covariates (AIR1996, baiocchi2014instrumental). In general, even a valid IV cannot point identify the mean conditional potential outcomes $\mathbb{E}[Y(\pm 1)|\bsX=\bsx]$ or the CATE $\mathscr{C}^\star(\bsx)$ (robins1996identification, balke1997bounds, manski2003partial,swanson2018partial); hence, neither the value nor the risk of an ITR can be point identified with an IV without additional assumptions. wang2018bounded, cui2019semiparametric and qiu2020optimal establish a set of sufficient conditions, under which the CATE can be point identified and the optimal ITR can be estimated with an IV.
Despite the incapability of point identifying the causal effects, a valid IV can still be useful in that even under minimal identification assumptions, it can produce meaningful partial identification intervals/bounds for the CATE. That is, one can construct two functions $L(\bsx)$ and $U(\bsx)$, both of which are functionals of the observed data distribution (and thus estimable from the data), such that the CATE satisfies $\mathscr{C}^\star(\bsx)\in [L(\bsx), U(\bsx)]$ almost surely. Important examples include the Balke-Pearl bounds balke1997bounds and the Manski-Pepper bounds manski1998monotone.
Given the partial identification interval $[L(\bsx), U(\bsx)]$, the risk of an ITR $\pi$ defined in (ref) is bounded between
and
pu2020estimating argued that a sensible criterion is to minimize the expected worst-case risk $\overline{\mathcal{R}}^{\pi}$, and the resulting minimizer is termed an IV-optimal ITR, as this ITR is “worst-case risk-optimal” with respect to the partial identification interval induced by an IV and its associated identification assumptions. Note that an IV-optimal ITR is not an optimal ITR without further assumptions.
cui2021machine proposed an alternative set of optimality criteria from the perspective of the partially identified value function. Instead of constructing two functions $L(\bsx), U(\bsx)$ that contains the CATE, one may alternatively construct functions $\underline{Q}(\bsx, a)$, $\overline{Q}(\bsx, a)$ that sandwiches the conditional mean potential outcome $\bbE[Y(a)|\bsX=\bsx]$, and the value defined in (ref) satisfies $\bbE[\underline{Q}(\bsX, \pi(\bsX))]\leq V^\pi \leq \bbE[\overline{Q}(\bsX, \pi(\bsX))]$. cui2021machine advocated maximizing some carefully-chosen middle ground between the lower and upper bounds of the value:
where $\lambda(\bsx, a)$ is a pre-specified function that captures a second-level individualism, i.e., how optimistic/pessimistic individuals with covariates $\bsX = \bsx$ are.
As pointed out by cui2021machine, minimizing the maximum risk $\overline{R}^\pi$ is not equivalent to maximizing the minimum value $V^\pi_{\lambda=0}$, even when the bounds $L(\bsx), U(\bsx)$ of the CATE are obtained via bounds on the conditional values, i.e., $L(\bsx) = \overline{Q}(\bsx, +1) - \underline{Q}(\bsx, -1)$ and $U(\bsx) = \underline{Q}(\bsx, +1) - \overline{Q}(\bsx, -1)$. Rather, minimizing $\overline{R}^\pi$ is equivalent to maximizing the midpoint of the minimum and maximum values, namely $V^\pi_{\lambda= 1/2}$.
In this section, we present a general framework that incorporates the extra information in IVs for better policy estimation. Conceptually, a valid IV and the associated identification assumptions impose distributional constraints on potential outcome distributions $p^\star_a(\cdot|\bsx)$. For example, under assumptions leveraged in cui2019semiparametric and qiu2020optimal, $p^\star_a(\cdot|\bsx)$ can be expressed as functionals of the observed data distribution. As another example, if less stronger assumptions are imposed so that point identification is impossible, the partial identification results assert that $p^\star_a(\cdot|\bsx)$ is “weakly bounded”, in the sense that for a sufficiently regular function $f$, we can find two functions $\underline{Q}(\bsx, a; f)$ and $\overline{Q}(\bsx, a; f)$ such that
If $f$ is the identity function, then the above display is precisely the partial identification intervals of the mean conditional potential outcome $\bbE[Y(a)|\bsX=\bsx]$ that appeared in (ref). In both examples, a valid IV along with the identification assumptions allows us to specify a collection of distributions $\calP_{\bsx, a}$, so that $p^\star_a(\cdot|\bsx) \in \calP_{\bsx, a}$. In the first example, the set $\calP_{\bsx, a}$ is a singleton consisting of the ground truth potential outcome distribution $p^\star_a(\cdot|\bsx)$, whereas in the second example, the set $\calP_{\bsx, a}$ consists of all distributions that are weakly bounded in the sense of (ref). In words, $\calP_{\bsx, a}$ is the collection of all possible potential outcome distributions $p^\star_{a}(\cdot|\bsx)$ that are compatible with the putative IV and the associated IV identification assumptions, and we refer to it as an IV-constrained set.
We may further equip an IV-constrained set $\calP_{\bsx, a}$ with a prior distribution $\mathscr{P}_{\bsx, a}$. Here, $\mathscr{P}_{\bsx,a}$ is a probability distribution on $\calP_{\bsx, a}$, the latter of which itself is a collection of probability distributions. For readers with a machine learning background, this is reminiscent of Baxter's model of inductive bias learning baxter2000model: in his language, $\calP_{\bsx, a}$ is called an environment, whose elements are called tasks, and it is assumed that nature can sample a task from $\mathscr{P}_{\bsx, a}$, which is a probability distribution on the environment. To have a fully rigorous treatment, we equip $\calP_{\bsx, a}$ with a metric (e.g., Wasserstein metric) and work with the induced Borel sigma algebra (parthasarathy2005probability).
Given the prior distributions $\{\mathscr{P}_{\bsx, a}: \bsx\in\calX, a =\pm 1\}$, we can define the IV-constrained value of a policy $\pi$ as
where we write $p_{\pi(\bsx)} = p_{\pi(\bsx)}(\cdot|\bsx)$ for notational simplicity. From now on, we will refer to the collection of criteria given by maximizing $V^\pi_{\mathscr{P}}$ as IV-optimality. An ITR is said to be IV-optimal with respect to the prior distributions $\{\mathscr{P}_{\bsx, a}: \bsx\in\bbR^d, a =\pm 1\}$ and a policy class $\Pi$ if it maximizes $V^\pi_\mathscr{P}$ among all $\pi\in\Pi$. Flow chart in Figure (ref) summarizes this conceptual framework.
It is clear that the formulation in (ref) can be recovered by carefully choosing the prior distributions. In particular, let $\underline{p}_a(\cdot|\bsx)$ and $\overline{p}_a(\cdot|\bsx)$ be distributions that witness the partial identification bounds $\underline{Q}(\bsx, a)$ and $\overline{Q}(\bsx, a)$, respectively:
Then the criterion (ref) is recovered by considering the following two-point priors:
where $\delta_p$ is a point mass at $p$. As discussed near the end of Section (ref), setting $\lambda(\bsx, a)$ uniformly equal to $1/2$ recovers the original IV-optimality criterion considered in pu2020estimating that minimizes the worst-case risk (ref). In fact, under certain regularity conditions, one can show that the reverse is also true: the formulation (ref) for a specified collection prior distributions can also be recovered from (ref) by a careful choice of $\lambda$. In view of such an equivalence, the criterion (ref) should not be regarded as a generalization of (ref). Rather, it is a convenient tool amenable to being generalized to the multiple-stage setting. A proof of this equivalence statement will appear in Section (ref). We defer how to estimate an IV-optimal ITR to Section (ref).
Compared to estimating an optimal ITR, a less ambitious goal it to improve upon a baseline ITR $\pi^{\textnormal{b}}$, so that the improved ITR is no worse and potentially better than $\pi^{\textnormal{b}}$. In this section, we show how to achieve this goal with an IV in the single-stage setup, and prepare readers for our main results concerning policy improvement in the general multiple-stage setup in Section (ref).
The goal of “never being worse” is reminiscent of the min-max risk criterion in (ref). In view of this, it is natural to consider minimizing the maximum excess risk with respect to the baseline ITR $\pi^{\textnormal{b}}$, subject to the IV-informed partial identification constraints. This strategy is summarized in the following definition.
When $\Pi$ consists of all Boolean functions, (ref) admits the following explicit solution.
The ITR in (ref) admits a rather intuitive explanation: it takes action $+1$ when the partial identification interval is positive (i.e., $0 < L(\bsx) \leq U(\bsx)$), and it takes action $-1$ when the interval is negative (i.e., $L(\bsx) \leq U(\bsx) < 0$), and it follows the baseline ITR $\pi^{\textnormal{b}}$ otherwise.
A closely related criterion to that appeared in Definition (ref) is to maximize the minimum “excess value” with respect to $\pi^{\textnormal{b}}$, subject to the distributional constraints imposed by the putative IV and its associated identification assumptions. This strategy is detailed as follows.
The above formulation is known as distributionally robust optimization in the optimization literature delage2010distributionally. When the IV-constrained sets are derived from partial identification intervals $[\underline{Q}(\bsx, a), \overline{Q}(\bsx, a)]$ and $\Pi$ consists of all Boolean functions, the optimization problem (ref) admits the explicit solution below, analogous to the one given in Proposition (ref).
Propositions (ref) and (ref) together reveal an interesting duality between worst-case excess risk minimization and worst-case excess value maximization. When the partial identification interval for the CATE is derived directly from the partial identification intervals for the mean conditional potential outcomes, so that $L^{\textnormal{\texttt{V}}}= L$ and $U^{\textnormal{\texttt{V}}} = U$, then the two ITRs defined in (ref) and (ref) agree, and they simultaneously satisfy the two IV-improvement criteria presented in Definitions (ref) and (ref). Hence, we will not distinguish between two types of IV-improved ITRs.
Such a duality is reminiscent of the duality between risk minimization and value maximization in the classical policy estimation problem under the NUCA discussed in Section (ref). Moreover, as discussed in Section (ref), minimizing maximum risk and maximizing minimum value subject to IV-informed partial identification intervals are not equivalent. Thus, it is curious to see that such a duality is restored in the policy improvement problem. We defer a discussion of how to estimate IV-improved ITRs to the more general dynamic treatment regimes setting studied in Section (ref).
We conclude this section with a comparison between the IV-improved ITR and the improved ITR derived in a series of works by kallus2018interval,kallus2018confounding,kallus2020minimax. In this line of work, the authors considered minimizing the maximum excess risk subject to a “Rosenbaum-bounds-type” sensitivity analysis model, and the expression of their improved ITR is similar to (ref) and (ref), but with $L(\bsx), U(\bsx)$ replaced by the bounds under the sensitivity analysis model. Despite the apparent similarity, there is a profound, practical difference between sensitivity-analysis-based and IV-based policy improvement. Since a sensitivity analysis model only relaxes the NUCA, the CATE derived under NUCA (i.e., $\mathscr{C}^\textnormal{\texttt{obs}}$) is always contained in a sensitivity analysis model. Hence, their framework can never improve upon the NUCA-optimal rule $\pi^\texttt{NUCA}$, defined as the minimizer of (ref). More explicitly, let $\textnormal{Imp}^\textnormal{\texttt{IV}}$ be a (potentially multi-valued) policy improvement operator that sends a baseline ITR $\pi^{\textnormal{b}}$ to its IV-improved counterparts, and let $\textnormal{Imp}^{\textnormal{\texttt{sens}}}$ be the corresponding policy improvement operator considered in kallus2018confounding. We necessarily have $\pi^\texttt{NUCA} \in \textnormal{Imp}^\textnormal{\texttt{sens}}(\pi^\texttt{NUCA})$, meaning that the NUCA-optimal rule can never be improved by $\textnormal{Imp}^\textnormal{\texttt{sens}}$. On the other hand, we will demonstrate via extensive simulations in Section (ref) that $\textnormal{Imp}^{\textnormal{\texttt{IV}}}(\pi^\texttt{NUCA})$ may yield a strictly better ITR than $\pi^\texttt{NUCA}$ as an IV contains additional information.
We now consider a general $K$-stage decision making problem murphy2003optimal,schulte2014q. At the first stage we are given baseline covariates $\bsX_1 = \bsx_1 \in \calX$, make a decision $a_1 = \pi_1(\bsx_1)$, and observe a reward $R_1(a_1) = r_1$. At stage $k = 2, \cdots, K$, let $\vec{a}_{k-1} = (a_1, \cdots, a_{k-1})$ denote treatment decisions from stage $1$ up to $k-1$, and $\vec{R}_{k - 1}(\vec{a}_{k-1}) = (R_1(a_1), R_2(\vec{a}_2), \cdots, R_{k-1}(\vec{a}_{k-1}))$ the rewards from stage $1$ to $k-1$. Meanwhile, denote by $\bsX_k(\vec{a}_{k - 1})$ new covariate information (e.g., time-varying covariates, auxiliary outcomes, etc) that would arise after stage $k-1$ but before stage $k$ if the patient were to follow the treatment decisions $\vec{a}_{k-1}$, and $\Vec{\bsX}_k(\vec{a}_{k-1}) = (\bsX_1, \bsX_2(a_1), \cdots, \bsX_k(\vec{a}_{k-1}))$ the entire covariate history up to stage $k$. Given the realizations of $\vec R_{r-1} = \vec r_{k-1}, \vec\bsX_{k} = \vec\bsx_k$, a possibly data-driven decision $a_k = \pi_k(\vec{a}_{k - 1}, \vec{r}_{k-1}, \vec{\bsx}_k)$ is then made and we observe reward $R_k(\vec a_k) = r_k$. Our goal is to estimate a dynamic treatment regime $\pi = \{\pi_k\}_{k=1}^K$, such that the expected value of cumulative rewards $\sum_{k=1}^K r_k$ is maximized if the population were to follow $\pi$.
To simplify the notations, the historical information available for making the $k$-th decision at stage $k$ is denoted by $ \vec{\bsH}_1 = \bsX_1 $ for $k = 1$ and $ \vec\bsH_k = \big(\vec a_{k-1}, \vec R_{k-1}(\vec a_{k-1}), \vec\bsX_k(\vec a_{k-1})\big) $ for any $ 2\leq k \leq K. $ The information $\vec\bsH_k$ consists of $\vec \bsH_{k-1}$ and $\bsH_k $: $\vec\bsH_{k-1}$ is the information available at the previous stage and $\bsH_k = \big(a_{k-1}, R_{k-1}(\vec{a}_{k-1}), \bsX_k(\vec a_{k-1})\big)$ is the new information generated after decision $a_{k-1}$ is made. Let $p^\star_{a_K}(\cdot|\vec\bsh_K)$ be the conditional law of $R(\vec a_K)$ given a specific realization of the historical information $\vec\bsH_K = \vec\bsh_k$. We define the following action-value function (or $Q$-function) at stage $K$:
where we emphasize that the expectation is taken over the potential outcome distribution $p^\star_{a_K}(\cdot|\vec\bsh_K)$. For a policy $\pi$, its value function at stage $K$ is then defined as
Note that the value function at stage $K$ depends on $\pi$ only through $\pi_K$.
Next, we define $Q$-functions and value functions at a generic stage $k\in[K]$ recursively (denote $[K] = \{1, \hdots, K\}$). In particular, for $k = K-1, \hdots, 1$, we let $p^\star_{a_k}(\cdot, \cdot | \vec\bsh_k)$ denote the joint law of the reward $R_k(\vec a_k)$ and new covariate information $\bsX_{k+1}(\vec a_k)$ observed immediately after decision $a_k$ has been made conditional on $\vec\bsH_k = \vec\bsh_k$. The $Q$-function of $\pi$ at stage $k$ are then defined as:
where $\vec \bsH_{k+1} = (\vec\bsh_k , a_{k}, R_{k}, \bsX_{k+1})$ is a function of $R_k$ and $\bsX_{k+1}$. The corresponding value function of $\pi$ is taken to be $V^\pi_k(\vec\bsh_k) = Q^\pi_k\big(\vec\bsh_k, \pi_k(\vec\bsh_k)\big).$ Similar to (ref), the expectation is taken over the potential outcome distribution $p^\star_{a_k}(\cdot, \cdot|\vec\bsh_k)$. We interpret the $Q$-function $Q^\pi_k(\vec\bsh_k, a_k)$ as the cumulative rewards collected by executing $a_k$ at stage $k$ and follow $\pi$ from stage $k+1$ and onwards. In contrast, the value function $V^\pi_k(\vec\bsh_k)$ is the cumulative rewards collected by executing $\pi$ from stage $k$ and onwards. Therefore, $Q^\pi_k$ depends on $\pi$ only through $\pi_{(k+1):K} = (\pi_{k+1},\hdots, \pi_K)$ and $V^\pi_k$ depends on $\pi$ only through $\pi_{k:K}=(\pi_k, \hdots, \pi_K)$.\footnote{For notational simplicity we interpret quantities whose subscripts do not make sense (e.g., $a_0, r_0$ and $\bsx_{K+1}$) as “null” quantities and their occurrences in mathematical expressions will be disregarded.}
With a slight abuse of notation, we let $\Pi = \Pi_1\times \cdots \Pi_K$ be a policy class. A DTR is said to be optimal with respect to $\Pi$ if it maximizes $V^\pi_1(\bsx_1)$ for all fixed $\bsx_1$ (and hence $\bbE[V^\pi_1(\bsX_1)]$ for an arbitrary covariate distribution) over the policy class $\Pi$.
Assume for now that $\Pi$ consists of all DTRs (i.e., each $\Pi_k$ consists of all Boolean functions). A celebrated result from control theory states that the dynamic programming approach below, also known as backward induction, yields an optimal DTR $\pi^\star$ (murphy2003optimal, sutton2018reinforcement):
More explicitly, $\pi^\star$ satisfies $V^{\pi^\star}_k(\vec\bsh_k) = \max_{\pi} V^\pi_k(\vec\bsh_k)$ for any $k\in[K]$ and any configuration of the historical information $\vec\bsh_k$, and is always well-defined as $Q^{\pi^\star}_k$ depends on $\pi^\star$ only through $\pi^\star_{(k+1):K}$.
The foregoing discussion is based on the potential outcome distributions. Suppose that we have collected i.i.d. data from the law of the random trajectory $(\bsX_k^\textnormal{\texttt{obs}}, A_k^\textnormal{\texttt{obs}}, R_k^\textnormal{\texttt{obs}})_{k = 1}^K$. Let $\{p^\textnormal{\texttt{obs}}_{a_k}(\cdot, \cdot|\vec\bsh_k)\}_{k=1}^{K-1}$ be the conditional laws of the observed rewards and new covariate information identified from the observed data:
for $k\leq K-1$. Here, $\vec\bsH_k^\textnormal{\texttt{obs}}$ denotes the observed historical information observed up to stage $k$. Suppose that we obtain a DTR using the dynamic programming approach described in (ref) with $p^\textnormal{\texttt{obs}}_{a_k}$ in place of $p^\star_{a_k}$. Under a version of the sequential randomization assumptions (robins1998marginal) and with additional assumptions of consistency and positivity, we have $ p^\star_{a_k} = p^\textnormal{\texttt{obs}}_{a_k}, $ from which it follows that the DTR obtained is in fact an optimal DTR murphy2003optimal, schulte2014q. In view of this fact, we will refer to this policy as an SRA-optimal DTR and denote it as $\pi_\textnormal{\texttt{SRA}}$. Methods that estimate SRA-optimal DTRs have been well studied in the literature. See murphy2003optimal, zhao2015new,tao2018tree,zhang2018c, among many others.
Suppose that in addition to the observed trajectory $(\bsX_k^\textnormal{\texttt{obs}}, A_k^\textnormal{\texttt{obs}}, R_k^\textnormal{\texttt{obs}})_{k = 1}^K$, we have access to a time-varying instrumental variable $\{Z_k\}_{k = 1}^K$. Similar to the single-stage setting in Section (ref), the IV, along with its associated identification assumptions, imposes distributional constraints on the potential outcome distributions $p^\star_{a_k}(\cdot, \cdot|\vec\bsh_k)$. That is, we can specify, for each action $a_k$ and historical information $\vec\bsh_k$ at stage $k$, an IV-constrained set $\calP_{\vec\bsh_k, a_k}$, which contains the ground truth potential outcome distribution $p^\star_{a_k}(\cdot, \cdot | \vec\bsh_k)$. Again, two primary examples are that $\calP_{\vec\bsh_k, a_k}$ is a singleton under point identification (michael2020instrumental), and that $\calP_{\vec\bsh_k, a_k}$ contains weakly bounded distributions by the partial identification intervals in the sense of (ref). For ease of exposition, we treat $\calP_{\vec\bsh_k, a_k}$ as a generic set of distributions compatible with the IV and identification assumptions for now; examples of $\calP_{\vec\bsh_k, a_k}$ will be given in Section (ref) when we formally describe estimation procedures.
It is essential to have a time-varying IV (e.g., daily precipitation) rather than a time-independent IV (e.g., sickle cell trait) in the multiple-stage setting; see Supplementary Materials (ref) for details.
Given the IV-constrained set $\calP_{\vec\bsh_k, a_k}$, we impose a prior distribution $\mathscr{P}_{\vec\bsh_k, a_k}$ on it. This notion generalizes the single-stage setting in Section (ref). We use $p_{a_k}(\cdot,\cdot|\vec\bsh_k) \sim \mathscr{P}_{\vec\bsh_k, a_k}$ to denote sampling a distribution $p_{a_k}(\cdot, \cdot |\vec\bsh_k )\in \calP_{\vec\bsh_k, a_k}$ from $\mathscr{P}_{\vec\bsh_k, a_k}$. We will use the shorthand $p_{a_k}$ for $p_{a_k}(\cdot,\cdot|\vec\bsh_k)$ where there is no ambiguity.
Under these notation, we introduce IV-constrained counterparts of the conventional $Q$- and value functions defined in (ref)--(ref).
Definition (ref) generalizes the notion of IV-constrained value in the single-stage setting. Similar to the conventional $Q$- and value functions, the IV-constrained $Q$-function $Q^\pi_{\mathscr{P},k}$ depends on $\pi$ only through $\pi_{(k+1):K}$ and the IV-constrained value function $V^\pi_{\mathscr{P}, k}$ depends on $\pi$ only through $\pi_{k:K}$.
Definition (ref) follows from Definition (ref) and defines the notion of IV-optimality for DTRs.
Note that a priori, we do not know if an IV-optimal policy $\pi^\star$ exists, as the above definition requires $\pi^\star$ to maximize the IV-constrained value function for every fixed $\bsx_1$ (and thus $\bbE_{\bsX_1}[V^\pi_{\mathscr{P}, 1}(\bsX_1)]$ for any law of $\bsX_1$). However, as we will see in our first main result below, the optimization problem (ref) can be solved via a modified dynamic programming algorithm if $\Pi$ consists of all policies.
The above result is a generalization of the classical dynamic programming algorithm (ref) under partial identification. Note that the policy $\pi^\star$ in the above theorem is always well-defined as $Q^{\pi^\star}_{\mathscr{P}_k}$ depends on $\pi^\star$ only through $\pi^\star_{(k+1):K}$.
Estimating an IV-optimal DTR based on Theorem (ref) alone is difficult, primarily because the definitions of the IV-constrained $Q$- and value functions involve integration over the prior distributions $\{\mathscr{P}_{\vec\bsh_k, a_k}\}$. In this subsection, we provide an alternative characterization of an IV-optimal DTR which illuminates a practical estimation strategy.
Recall that in the single-stage setting, there is an equivalence between the IV constrained value (ref) and the convex combination of the lower and upper bounds of the value (ref). The latter expression (ref) is easier to compute as long as the weights of the convex combination are specified. Below, we extend such an equivalence relationship to the multiple-stage setting and use it to get rid of the intractable integration over prior distributions.
We start by defining the worst-case and best-case IV-constrained $Q$- and value functions as well as their weighted versions as follows.
Definition (ref) generalizes the criterion (ref). According to Definition (ref), the worst-case and best-case $Q$-functions $\underlineQ^\pi_{\vec{\lambda}, k}$ and $\overlineQ^\pi_{\vec{\lambda}, k}$ depend on the weighting functions only through $\lambda_t(\vec\bsh_{t}, a_t)$ for $k+1 \leq t \leq K\}$, whereas the corresponding value functions depend on the weighting functions only through $\lambda_t(\vec\bsh_{t}, a_t)$ for $k \leq t \leq K$. For notational simplicity, in the rest of the paper, we add superscript $\pi$ and subscript $\vec\lambda$ in $Q$-functions at stage $K$ (e.g., we write $\underline{Q}_K = Q^\pi_{\vec\lambda, K}$), although they have no dependence on $\pi$ and the weighting functions.
Proposition (ref) below establishes the equivalence between the weighted $Q$-functions (resp. value functions) and their IV-constrained counterparts in Definition (ref).
The above theorem effectively translates the problem of specifying a collection of prior distributions to specifying a collection of weighting functions, thus allowing one to bypass the integration over the prior distributions. This theorem, along with the dynamic programming algorithm presented in Theorem (ref), leads to the following alternative characterization of the IV-optimal DTR.
Compared to Theorem (ref), Corollary (ref) gives an alternative, analytic and constructive characterization of the IV-optimal DTR. The proof is an immediate consequence of Theorem (ref) and Proposition (ref), and omitted.
\paragraph{An illustrative example when $K=2$.} Figure (ref) illustrates the decision process when $K = 2$. At the second stage, given $\vec{\bsh}_2^+ = ({\bsh}_1, a_1=+1, R_1(+1) = r_1^+, \bsX_2(+1) = \bsx_2^+)$ and $\lambda_2(\vec\bsh_2^+, \pm 1)$, an IV-optimal action is made based on comparing two weighted $Q$-functions $Q_{\vec\lambda, 2}(\vec\bsh_2^+, +1) = 0.8$ and $Q_{\vec\lambda, 2}(\vec\bsh_2^+, -1) = 0.5$. Since $0.8>0.5$, we have $\pi^\star_2(\vec\bsh_2^+) = +1$. The decision given $\vec\bsh_2^- = (\bsh_1, a_1 = -1, R_1(-1) = r_1^-, \bsX_2(-1) = \bsx_2^-)$ and $\lambda_2(\vec\bsh_2^-, \pm 1)$ is similar, and we have $\pi^\star_2(\vec\bsh_2^-) = -1$ in this case.
Given the knowledge of $\pi^\star_2$, we can compute the weighted value function $V^{\pi^\star}_{\vec\lambda, 2}(\cdot)$ of $\pi^\star$ at the second stage. This allows us to construct a “pseudo-outcome” at the first stage, denoted as $PO_1(r_1, \bsx_2|{\bsh}_1, a_1) = r_1 + V_2^{\pi^\star}({\bsh}_1, a_1, r_1, \bsx_2)$. Importantly, $PO_1$ depends on $\pi^\star_2$, as it is the cumulative rewards if we observe $\bsh_1$ at the first stage, take an immediate action $a_1$, and then act according to the IV-optimal decision $\pi^\star_2$ at the second stage.
The partial identification interval for the expected value of $PO_1(R_1(a_1), \bsX_2(a_1)|\bsh_1, a_1)$, where the expectation is taken over the potential outcome distribution of $R_1(a_1)$ and $\bsX_2(a_1)$, is precisely the worst-case and best-case $Q$-functions $\underline{Q}^{\pi^\star}_{\vec\lambda, 1}(\bsh_1, a_1)$ and $\overline{Q}^{\pi^\star}_{\vec\lambda, 1}(\bsh_1, a_1)$. To this end, we can compute the weighted $Q$-function $Q^{\pi^\star}_{\vec\lambda,1}(\bsh_1, a_1)$ if $\lambda_1(\bsh_1, a_1)$ is specified, from which we can decide $\pi^\star_1$. Reading the numbers off Figure (ref), we conclude that $\pi^\star_1(\bsh_1) = -1$.
\paragraph{Choosing weighting functions.} The specification of weighting functions reveals one's level of optimism. Suppose that the future weighting functions $\{\lambda_t: t\geq k+1\}$ has been specified. At stage $k$, if one adopts a worst-case perspective and would like to maximize the worst-case gain at this stage (fixing the weighting function specifications at all future stages), then it suffices to compare the two worse-case $Q$-functions at stage $k$, namely $\underline{Q}_{\vec{\lambda} ,k}^{\pi^\star} (\vec\bsh_k, +1)$ and $\underline{Q}_{\vec{\lambda} ,k}^{\pi^\star} (\vec\bsh_k, -1)$. And the pessimistic action is taken to be sign of the difference of the two worst-case $Q$-functions, which corresponds to taking $\lambda_k(\vec\bsh_k, \pm 1) = 1$ (see Figure (ref)). Alternatively, if one adopts a best-case perspective at stage $k$ and would like to maximize the best-case gain at this stage, then one shall compare the two best-case $Q$-functions, namely $\overline{Q}_{\vec{\lambda} ,k}^{\pi^\star} (\vec\bsh_k, +1)$ and $\overline{Q}_{\vec{\lambda} ,k}^{\pi^\star} (\vec\bsh_k, -1)$. The optimistic action is then taken to be the difference between the two best-case $Q$-functions, which corresponds to taking $\lambda_k(\vec\bsh_k, \pm 1) = 0$ (see Figure (ref)). Recall that in the single-stage setting, the min-max risk criterion (ref) corresponds to maximizing the weighted value (ref) with weights set to be $1/2$. This criterion can be seamlessly generalized to the current multiple-stage setting by setting $\lambda_k(\vec\bsh_k, \pm 1) = 1/2$ (see Figure (ref)). Other choices of weighting functions can be made to incorporate domain knowledge and user preference, as suggested by cui2021machine.
We discuss how to estimate an IV-optimal DTR given i.i.d. samples from the law of the random trajectory $(\bsX_k^\textnormal{\texttt{obs}}, A_k^\textnormal{\texttt{obs}}, R_k^\textnormal{\texttt{obs}})_{k = 1}^K$ and a time-varying instrument variable $\{Z_k\}_{k=1}^K$. The IV-constrained sets $\{\calP_{\vec\bsh_k, a_k}\}$ are specified via partial identification intervals for the $Q$-functions. Specifically, let
The boundedness of $R_k$ ensures that partial identification intervals (e.g, Manski-Pepper bounds) have finite width. This is a plausible assumption in many applications swanson2018partial.
At the final stage $K$, we take $\calF_K = \{\textnormal{Id}\}$, a singleton containing the identity function that sends $r_K$ to itself (note that $\bsX_{x+1}$ is a null quantity and thus disregarded). Then the lower and upper bounds, namely $\underline{Q}(\vec\bsh_k, a_K; \textnormal{Id}) = \underline{Q}_K(\vec\bsh_k, a_K)$ and $\overline{Q}(\vec\bsh_k, a_K; \textnormal{Id}) = \underline{Q}_K(\vec\bsh_k, a_K)$, are precisely the endpoints of the partial identification intervals (constructed using the IV $Z_K$) of the expected final-stage reward $R_K$, where the expectation is taken over the potential outcome distribution $R_K\sim p^\star_{a_K}(\cdot|\vec\bsh_k)$. By construction, we have $p^\star_{a_K}(\cdot|\vec\bsh_k) \in \calP_{\vec\bsh_k, a_K}$.
For $k\leq K-1$, we will take $\calF_k$ to be a class of properly defined functions (the specific forms to be specified later) of the current reward $r_k$ and the next-stage covariate $\bsx_{k+1}$. The lower and upper bounds $\underline{Q}(\vec\bsh_k, a_k; f), \overline{Q}(\vec\bsh_k, a_k; f)$ constitute the partial identification intervals (constructed using the IV $Z_k$) of the expected value of $f(R_k, \bsX_{k+1})$, where the expectation is taken over the potential outcome distribution $(R_k, \bsX_{k+1})\sim p^\star_{a_k}(\cdot, \cdot|\bsh_k)$. By construction, we have $p^\star_{a_k}(\cdot, \cdot|\vec\bsh_k) \in \calP_{\vec\bsh_k, a_k}$.
Let $\{(\bsx_{k, i}, a_{k, i}, r_{k, i})_{k=1}^K \}_{i=1}^n$ denote the observed data. Let $\bsh_{1, i}= \bsx_{1, i}$ be the baseline covariates for the $i$-th sample. For $2 \leq k \leq K$, let $\vec\bsh_{k, i}= (\vec a_{k-1, i}, \vec r_{k-1, i}, \vec\bsx_{k, i})$ be the historical information up to stage $k$ for the $i$-th sample. In addition, let $\{z_{k, i}\}_{i=1}^n$ denote the stage-$k$ IV data.
\paragraph{Estimating the contrasts by $Q$-learning.} Corollary (ref) shows that to estimate the IV-optimal DTR $\pi^\star$, it suffices to estimate the contrast functions $\{\mathscr{C}^{\pi^\star}_{\vec\lambda, k}(\vec\bsh_k)\}$ defined in (ref), which in turn calls for estimating the weighted $Q$-functions $\{Q^{\pi^\star}_{\vec\lambda, k}(\vec\bsh_k, \pm 1)\}$. We use a $Q$-learning approach to estimate the weighted $Q$-functions watkins1992q,schulte2014q.
For ease of exposition, we specify (ref) via Manski-Pepper bounds manski1998monotone. Generalization to other types of partial identification intervals is immediate. For Manski-Pepper bounds to hold, we make the mean exchangeability assumption (manski1990nonparametric, hernan2006instruments, swanson2018partial) or a relaxed monotone instrumental variable (MIV) assumption (manski1998monotone). The mean exchangeability assumption is automatically satisfied in a sequential randomized controlled trial with noncompliance.
At Stage $K$, define
Manski-Pepper bounds state that if $R_K \in [\underline{C}_K, \overline{C}_K]$ almost surely (with respect to $p^\star_{a_K}(\cdot|\vec\bsh_k)$), then its conditional mean potential outcome $\bbE_{R_K\sim p^\star_{a_K}(\cdot|\vec\bsh_k)}[R_K]$ is lower and upper bounded by
where $\land$ and $\lor$ are shorthands for $\min$ and $\max$, and both bounds are tight. Therefore, as long as we take $\calF_K \supseteq \{\text{Id}\}$ when defining $\calP_{\vec\bsh_k, a_K}$ in (ref), the worst-case and best-case $Q$ functions at stage $K$ defined in (ref)--(ref) can be set to (ref) and (ref), respectively. Along with $\lambda_K(\vec\bsh_k, a_K)$ specifications, the construction of the stage-$K$ weighted $Q$-function $Q_{\vec\lambda, K}(\vec\bsh_k, a_K)$ is concluded.
Since both (ref) and (ref) are functionals of the observed data distribution, $Q_{\vec\lambda, K}(\vec\bsh_k, a_K)$ can be estimated from the data by fitting parametric models (e.g., linear models) or flexible machine learning models (e.g., regression trees and random forests), and then invoke the plug-in principle.
Given an estimate $\widehat Q_{\vec\lambda, K}^\star$ of $Q_{\vec\lambda, K}$, we can estimate the contrast function $\mathscr{C}_{\vec\lambda, K}(\vec\bsh_k)$ by $\widehat \mathscr{C}^\star_{\vec\lambda, K}(\vec\bsh_k) = \widehat Q_{\vec\lambda, K}^\star(\vec\bsh_k, +1) - \widehat Q_{\vec\lambda, K}^\star(\vec\bsh_k, -1)$. In view of Corollary (ref), the IV-optimal DTR at stage $K$, $\pi^\star_K(\vec\bsh_k)$, can be estimated by $\widehat\pi^{\textnormal{\texttt{Q}}}_K(\vec\bsh_k) = \sgn(\widehat\mathscr{C}_{\vec\lambda, K}(\vec\bsh_k))$. Moreover, the weighted value function of $\pi^\star$ at stage $K$, $V^{\pi^\star}_{\vec\lambda, K}(\vec\bsh_k)$, can be estimated by $\widehat V^{\star}_{\vec\lambda, K} = \widehat Q^\star_{\vec \lambda, K}(\vec\bsh_k, \widehat\pi^{\textnormal{\texttt{Q}}}_K(\vec\bsh_k))$.
Now, assume for any stage $t\geq k+1$, the specification of the IV-constrained sets (ref) has been made, and the weighted value function at stage $t$, $V^{\pi^\star}_{\vec\lambda, t}(\vec\bsh_t)$, has been estimated by $\widehat V^\star_{\vec\lambda, t}(\vec\bsh_t)$. In addition, assume that $R_t \in [\underline{C}_t, \overline{C}_t]$ almost surely (with respect to $p^\star_{a_k}(\cdot, \cdot|\vec\bsh_k)$). At stage $k$, define the pseudo-outcome $PO_k(r_k, \bsx_{k+1}|\vec\bsh_k, a_k) = r_k + V^{\pi^\star}_{\vec\lambda, k+1}(\vec\bsh_k, a_k, r_k, \bsx_{k+1})$. By construction, we have $PO_k(R_k, \bsX_{k+1}|\vec\bsh_k, a_k) \in [\sum_{t\geq k}\underline{C}_t, \sum_{t\geq k} \overline{C}_t]$ almost surely for $(R_k, \bsX_{k+1}) \sim p^\star_{a_k}(\cdot, \cdot|\vec\bsh_k)$. Thus, we can apply Manski-Pepper bounds again to bound the expected value of $PO_k(R_k, \bsX_{k+1}|\vec\bsh_k, a_k)$, and obtain an estimate $\widehat Q_{\vec\lambda, k}^\star$ of the weighted $Q$-function of $\pi^\star$ at stage $k$; see Supplementary Material (ref) for detailed expressions. One nuance is that the pseudo-outcome $PO_k(r_k, \bsx_{k+1}|\vec\bsh_k, a_k)$ depends on the unknown quantity $V^{\pi^\star}_{\vec\lambda, k+1}$. At stage $k$, we have already obtained an estimate $\widehat V^{\star}_{\vec\lambda, k+1}$ of $V^{\pi^\star}_{\vec\lambda, k+1}$. Thus, the pseudo-outcome can be estimated by $\widehat{PO}_k(r_k, \bsx_{k+1}|\vec\bsh_k, a_k) = r_k + \widehat V^{\star}_{\vec\lambda, k+1}(\vec\bsh_k, a_k ,r_k, \bsx_{k+1})$.
Finally, the contrast function $\mathscr{C}_{\vec\lambda, k}^{\pi^\star}(\vec\bsh_k)$ is estimated by $\widehat \mathscr{C}^\star_{\vec\lambda, k}(\vec\bsh_k) = \widehat Q_{\vec\lambda, k}^\star(\vec\bsh_k, +1) - \widehat Q_{\vec\lambda, k}^\star(\vec\bsh_k, -1)$, the IV-optimal DTR at stage $k$ is estimated by $\widehat\pi^{\textnormal{\texttt{Q}}}_k(\vec\bsh_k) = \sgn(\widehat\mathscr{C}_{\vec\lambda, k}(\vec\bsh_k))$, and the weighted value function of $\pi^\star$ at stage $k$ is estimated by $\widehat V^{\star}_{\vec\lambda, k} = \widehat Q_{\vec \lambda, k}(\vec\bsh_k, \widehat\pi^{\textnormal{\texttt{Q}}}_k(\vec\bsh_k))$. In this way, we recursively estimate all contrast functions and obtain an estimated IV-optimal DTR $\widehat\pi^{\textnormal{\texttt{Q}}}$.
\paragraph{Obtaining parsimonious policies via weighted classification.} In many applications, it is desirable to impose additional constrains on the estimated DTR. For example, one may require the DTR to be parsimonious and thus more interpretable. Such constraints are usually encoded by a restricted function class $\Pi = \Pi_1 \times \cdots \times \Pi_K$, such that any $\pi\notin \Pi$ will not be considered.
We next introduce a strategy that “projects” $\widehat\pi^{\textnormal{\texttt{Q}}}$, the DTR obtained by $Q$-learning and is thus not necessarily inside $\Pi$, onto the function class $\Pi$. Recall the formula for the IV-optimal DTR $\pi^\star$ given in (ref). With some algebra, one readily checks that $\pi^\star$ is a solution to the following weighted classification problem:
The above representation illuminates a rich class of strategies to search for a DTR $\pi$ within a desired, possibly parsimonious function class $\Pi$, via sequentially solving the following weighted classification problem:
where we recall that $\widehat\mathscr{C}^\star_{\vec\lambda, k}$ is an estimate of $\widehat\mathscr{C}^{\pi^\star}_{\vec\lambda, k}$ obtained via $Q$-learning, and $\{\vec\bsh_{k, i}\}_{i=1}^n$ is the historical information at stage $k$ in our dataset. This strategy bears similarities with the strategy proposed in zhao2015new. The estimation procedure is summarized in Algorithm (ref).
In this section, we show how the IV-optimality framework developed in Section (ref) can be modified to tackle the policy improvement problem under the multiple-stage setup. Let $\pi^{\textnormal{b}}$ be a baseline DTR to be improved. Some most important baseline DTRs include the standard-of-care DTR (i.e., $\pi^{\textnormal{b}}_k = -1$ for any $k$) and the SRA-optimal DTR defined in Section (ref). The goal of policy improvement, as discussed in Section (ref), is to obtain a DTR $\pi^{\uparrow}$ so that $\pi^{\textnormal{b}}$ is no worse and potentially better than $\pi^{\textnormal{b}}$.
To achieve this goal, we leverage additional information encoded in the collection of IV-constrained sets $\{\calP_{\vec\bsh_k, a_k}\}$. We are to introduce the notion of IV-improved DTR, generalizing the notion of IV-improved ITR in Definitions (ref) and (ref). To start with, we define the following {relative $Q$-functions} and the corresponding {relative value functions} of a DTR $\pi$ with respect to the baseline DTR $\pi^{\textnormal{b}}$.
The relative value function $V^{\pi/\pi^{\textnormal{b}}}_k$ at stage $k$ captures the cumulative worst-case (subject to the IV constraints) excess value of following the DTR $\pi$ over $\pi^{\textnormal{b}}$ from stage $k$ and onwards. Note that the relative value function $V^{\pi/\pi^{\textnormal{b}}}_K$ at stage $K$ is analogous to the objective function in (ref). Maximizing the relative value functions would then deliver a DTR that follows the baseline regime $\pi^{\textnormal{b}}$ unless there is compelling evidence not to, similar to the IV-improved ITR whose explicit form is given in Propositions (ref) and (ref).
We are now ready to define the estimand of interest, generalizing Definitions (ref) and (ref).
Analogous to Theorem (ref), the following result solves the optimization problem (ref) using a dynamic programming approach when $\Pi$ is the collection of all DTRs.
The following corollary is similar to Corollary (ref), and gives an alternative characterization of the IV-improved DTR.
Expressions in (ref) and (ref) admit a rather intuitive explanation. The contrast function $\mathscr{C}^{\pi^{\uparrow}/\pi^{\textnormal{b}}}_k$ measures the cumulative worst-case gain of following $\pi^{\uparrow}$ over $\pi^{\textnormal{b}}$, starting from stage $k$, and one would only flip the decision made by $\pi^{\textnormal{b}}_k$ if this worse-case gain is positive. In fact, with some algebra, one readily checks that $\pi^{\uparrow}_K$ can be expressed in a similar form as (ref) and (ref).
Given Corollary (ref) and the strategy in Section (ref), Algorithm (ref) can be modified mutatis mutandis to yield an algorithm for estimating an IV-improved DTR. For brevity, we refer the readers to Supplementary Material (ref) for details.
In this section, we prove non-asymptotic bounds on the deviance between the estimated IV-optimal as well as the IV-improved DTRs to their population counterparts. To do this, we need several standard assumptions, which we detail below.
First of all, we need a proper control on the complexity of the policy class $\Pi = \Pi_1\times \cdots \Pi_K$. Note that any $\pi_k \in \Pi_k$ is a Boolean functions that sends a specific configuration of the historical information $\vec\bsh_k$ to a binary decision $\pi_k(\vec\bsh_k)$. Let $\calH_k$ be the collection of all possible $\vec\bsh_k$s. A canonical measure of complexity of Boolean functions is the Vapnik–Chervonenkis (VC) dimension vapnik1968uniform. The VC-dimension of $\Pi_k$, denoted as $\textnormal{\texttt{vc}}(\Pi_k)$, is the largest positive integer $d$ such that there exists a set of $d$ points $\{\vec\bsh_{k}^{(1)}, \cdots, \vec\bsh_{k}^{(d)}\}\subseteq \calH_k$ shattered by $\Pi_k$, in the sense that for any binary vector $\bsv \in \{\pm 1\}^d$, there exists $\pi_{k}^{(\bsv)} \in \Pi_k$ such that $\pi_{k}^{(\bsv)}(\vec\bsh_k^{(j)}) = v_j$ for $j\in[d]$. If no such $d$ exists, then $\textnormal{\texttt{vc}}(\Pi_k) = \infty$.
In practice, a DTR is most useful when it is parsimonious. For this reason, the class of linear decision rules and the class of decision trees with a fixed depth are popular in various application fields tao2018tree,speth2020assessment. It is well-known that when the domain is a subset of $\bbR^d$, the VC-dimension of the class of linear decision rules is at most $d+1$ (see, e.g., Example 4.21 of wainwright2019high) and the VC dimension of the class of decision tree with $L$ leaves is $\calO(L \log (Ld))$ leboeuf2020decision.
Note that a DTR $\pi_{1:(k-1)}$, along with the distributions $\{p_{a_t}(\cdot, \cdot|\vec\bsh_t): a_t\in\{\pm 1\}, \vec\bsh_t \in \calH_t\}_{t=1}^{k-1}$, induces a probability distribution on $\calH_k$, the set of stage-$k$ historical information. Let $\scrH_k$ denote the set of all such probability distributions when we vary $\pi_{1:(k-1)} \in \Pi_{1:(k-1)}$ and $p_{a_t}(\cdot, \cdot|\vec\bsh_t) \in \calP_{\vec\bsh_t, a_t}$ for $t\in[k-1]$. An element in $\scrH_k$ will be denoted as $q_k$, and the law of $\vec\bsH_k^\textnormal{\texttt{obs}}$ will be denoted as $q_k^\textnormal{\texttt{obs}}$. The next assumption concerns how much information on $q_k^\textnormal{\texttt{obs}}$ can be generalized to the information on a specific $q_k \in \scrH_k$.
Assumption (ref) is often made in the reinforcement learning literature (see, e.g., munos2003error,szepesvari2005finite,chen2019information) and is closely related to the “overlap” assumption in causal inference. A sufficient condition for the above assumption to hold is that the probability of seeing any historical information (according to the observed data distribution or any distribution in $\scrH_k$) is strictly bounded away from $0$ and $1$.
Recall that our algorithms for estimating IV-optimal and IV-improved DTRs is a two-step procedure, where in Step \RN{1}, the contrast functions are estimated via $Q$-learning, and in Step \RN{2}, a parsimonious DTR is obtained via weighted classification. As mentioned in Section (ref), in the first stage, there are a lot of flexibilities in choosing the specific models and algorithms for estimating the contrast functions. As the result, a fine-grained understanding of this stage would require a case-by-case analysis. In order to not over-complicate the discussion, we make the following assumption, which asserts the existence of a “contrast estimation oracle”.
We will use (ref) when we analyze the estimated IV-optimal DTR and we will invoke (ref) when we analyze the estimated IV-improved DTR.
Assumption (ref) is usually a relatively mild assumption. Indeed, the ground truth contrast functions $C^{\pi^\star}_{\vec\lambda, k}$ are superimpositions of several conditional expectations of the observed data distribution (see (ref)--(ref)), and it is reasonable to assume that they can be estimated at a vanishing rate as the sample size tends to infinity. While this assumption simplifies the analysis by allowing us to bypass the case-by-case analyses of $Q$-learning, the analysis of the weighted classification remains highly non-trivial.
To proceed further, we adopt a form of sample splitting procedure called cross-fitting (chernozhukov2016double, athey2020policy). Specifically, we split the $n$ samples into $m$ equally-sized batches: $[n] = \cup_{j\in[m]} B_j$, where $m\geq 2$ and $m\asymp 1$. Each batch has size $|B_j| = n_j \asymp n/m \asymp n$. For each $i\in[n]$, let $j_i\in[m]$ be the index of its batch, so that $i\in B_{j_i}$. For the $i$-the sample, we apply the contrast estimation oracle in Assumption (ref) on the out-of-batch samples $B_{-j_i} = [n]\setminus B_{j_i}$ to obtain either an estimate $\widehat\mathscr{C}^\star_{\vec\lambda, k}(\vec\bsh_k; B_{-j_i})$ of $\mathscr{C}^{\pi^\star}_{\vec\lambda, k}(\vec\bsh_k)$ for policy estimation, or an estimate $\widehat \mathscr{C}^{\uparrow}_k(\vec\bsh_k; B_{-{j_i}})$ of $\mathscr{C}^{\pi^{\uparrow} /\pi^{\textnormal{b}}}_{\vec\lambda, k}$ for policy improvement. For the policy estimation problem, we solve the following optimization problem:
The corresponding optimization problem for the policy improvement problem is given by
{We emphasize that cross-fitting is mostly for theoretical convenience. The performance of our algorithm with and without cross-fitting is similar; see simulations in Section (ref).}
In practice, given a limited computational budget, we may only solve the optimization problem (ref) up to a certain precision. Our analysis will be conducted for the approximate minimizers $\widehat \pi^\star_k, \widehat{\pi}^{\uparrow}_k\in \Pi_k$ that satisfy
respectively, where $\ep^\textnormal{\texttt{opt}}_k$ is the optimization error when solving (ref) and (ref).
To quantify the loss of information from restriction to the policy class $\Pi$, we define
Note that according to Corollaries (ref) and (ref), when $\Pi_k$ is the set of all Boolean functions, we have $\widetilde \pi_k^\star = \pi^\star_k$ and $\widetilde{\pi}^{\uparrow}_k = \pi^{\uparrow}_k$. Otherwise, the difference between $\widetilde\pi^\star_k$ and $\pi^\star_k$ (resp. $\widetilde{\pi}^{\uparrow}_k$ and $\pi^{\uparrow}_k$) measures the approximation error when we restrict ourselves to $\Pi_k$ instead of the set of all Boolean functions in the policy estimation (resp. policy improvement) problem.
We are now ready to present the main result of this section.
Theorem (ref) shows that the weighted value function of $\widehat\pi^\star$ (resp. relative value function of $\widehat{\pi}^{\uparrow}$) at the first stage converges to that of the IV-optimal DTR (resp. IV-improved) up to three sources of errors: the optimization error that stems from only approximately solving the weighted classification problem, the approximation error that results from the restriction to a parsimonious policy class $\Pi$, and a vanishing generalization error. We defer the proof to Supplementary Material (ref).
We verify that the IV-optimal DTRs indeed have superior performance compared to the baseline DTRs and investigate the performance of the IV-optimal DTRs for assorted choices of $\lambda_k(\Vec{\bsh}_k, a_k)$ via simulations. We consider a data-generating process with two time-independent covariates $X_1$, $X_2 \sim \text{Unif}[-1, 1]$. At stage one, there exists an unmeasured confounder $U_1 \sim \textnormal{Bernoulli}(1/2)$. The instrumental variable $Z_1$ is independent of $U_1$ and follows $\textnormal{Rademacher}(1/2)$, the observed $A_1$ is Rademacher with head probability $\text{expit}\{C_1\cdot (Z_1+1) - \xi U_1 - 2\}$, and the reward $R_1$ is Bernoulli with head probability $\text{expit}\{0.5(\sgn\{X_1 - 1\} - \lambda U_1 + 0.2)\cdot (A_1 + 1)$, where $C_1, \xi, \lambda$ are constants to be specified later, and $\textnormal{expit}\{x\} = 1/(1+e^{-x})$. At stage two, there exists a second unmeasured confounder $U_2 \sim \textnormal{Bernoulli}(1/2)$. The instrumental variable is again independent of $Z_2$ and follows $\textnormal{Rademacher}(1/2)$, the action $A_2$ is Rademacher with head probability $\text{expit}\{C_1\cdot (Z_1+1) + X_1 - 7(R_1 - 0.5) - \xi(1 + X_1)(2U_2 - 1)\}$, and the reward $R_2$ is Bernoulli with head probability $\text{expit}\{0.1(A_1+1) + 0.4[1 - X_1 + R_1 - \lambda (2U_2 - 1)]\cdot (A_2 + 1)\}$.
According to this data-generating process, $Z_1$ and $Z_2$ are valid instrumental variables. There are two unmeasured confounders $U_1$ and $U_2$, one at each stage, and both unmeasured confounders are effect modifiers. We will interpret the treatment option $A_1, A_2 = -1$ as the standard-of-care, e.g., low-level NICU, and $+1$ as the prospective treatment, e.g., high-level NICU. One may check that for large $\xi$, the prospective treatment has a negative effect on $R_1$; however, the unmeasured confounding $U_1$ may create a spurious, positive treatment effect in an analysis that omits $U_1$. On the other hand, the prospective treatment at the second stage may have a positive or negative treatment effect on $R_2$ depending on the baseline covariate $X_1$ and the first stage outcome $R_1$.
Our simulation can be compactly summarized as a $3\times 3 \times 3 \times 2 \times 2$ factorial design with the following five factors:
The observed data consist of $\mathcal{D}_{\text{obs}} = \{(X_1, X_2, A_1, Z_1, R_1, A_2, Z_2, R_2),~i = 1, \cdots, n_{\textnormal{\texttt{train}}}\}$. An SRA-optimal DTR $\pi^{\textnormal{b}}_\textnormal{\texttt{SRA}}$ does not leverage the IV data ($Z_1$ and $Z_2$) while IV-improved DTR does and uses this information to improve upon baseline rules. We estimated the conditional average treatment effect involved in estimating $\pi^{\textnormal{b}}_{\textnormal{\texttt{SRA}}}$ using a robust augmented inverse probability weighting estimator (AIPW), and then applied a weighted classification routine. This is known as C-learning in the literature (zhang2018c), and can also be considered a variant of the backward outcome weighted learning (BOWL) procedure (zhao2015new). All classification problems involved in estimating the $\pi^{\textnormal{b}}_{\textnormal{\texttt{SRA}}}$ and the IV-improved DTRs were implemented using a classification tree with a maximum depth of $2$, which is meant to replicate the real data application where parsimonious rules are more useful and can deliver most insight (laber2015tree, speth2020assessment). Three IV-improved DTRs are denoted by $\pi^{\uparrow}_{\textnormal{\texttt{std}}}$, $\pi^{\uparrow}_{\textnormal{\texttt{prosp}}}$, and $\pi^{\uparrow}_{\textnormal{\texttt{SRA}}}$, respectively. Finally, for each data-generating process, we further estimated three IV-optimal DTRs with $\lambda_k(\Vec{\bsh}_k, \pm 1) = 0$, $1/2$, and $1$ for $k = 1, 2$, and these three IV-optimal DTRs are referred to as $\pi_{\texttt{IV}, 0}$, $\pi_{\texttt{IV}, 1/2}$, and $\pi_{\texttt{IV}, 1}$, respectively. Therefore, we have a total of nine regimes (three baseline regimes, three IV-improved regimes, and three IV-optimal regimes) under consideration. We evaluated each estimated regime $\widehat{\pi}$ by calculating its value function using $1,000,000$ fresh samples $(X_1, X_2) \sim \text{Unif}[-1, 1]$ and integrating out the binary unmeasured confounders $U_1, U_2$ using Monte Carlo.
Table (ref) summarizes the estimated mean and interquartile range of the value functions for baseline DTRs, their corresponding IV-improved DTRs, and three different IV-optimal DTRs, when $n_{\textnormal{\texttt{train}}} = 1000$, all relevant conditional expectations estimated via random forests (breiman2001random) as implemented in the R package randomForest, and for various $(\xi, C_1)$ combinations.
There are three trends consistent with our theory and intuition upon examining the simulation results. First and foremost, we observed that the IV-improved DTRs indeed had superior performance compared to their corresponding baseline DTRs including the SRA-optimal DTR. Although the extent of improvement depends on the specifics of data-generating processes, the improvement was uniform across all data-generating processes. Figure (ref) plots the cumulative distribution functions (CDFs) of three baseline DTRs and their IV-improved DTRs across $1,000$ simulations in two data-generating processes. It is evident that in either data-generating process and for any of the three baseline DTRs, the value functions of IV-improved DTRs always stochastically dominate those of corresponding baseline DTRs. We observed the same stochastic dominance phenomenon in each of the $108$ data-generating processes considered in the simulation studies. Second, when comparing three IV-optimal DTRs corresponding to different choices of weighting functions $\lambda_k(\Vec{\bsh}_k, a_k)$, we observed that the choices of $1/2$ and $1$, corresponding to the min-max and the worst-case perspectives, had better performance compared to the best-case DTR and the SRA-optimal DTR. We further plot CDFs of the value functions of each IV-optimal DTR and the SRA-optimal DTR (Figures (ref) and (ref) in Supplementary Material (ref)) and observed that the min-max and worst-case DTRs stochastically dominated the SRA-optimal DTRs in all sampling situations considered in the simulation studies. Lastly, we found that estimating relevant conditional expectations using simple parametric models (Tables (ref) and (ref) in Supplementary Material (ref)) and the cross-fitting version of the algorithm (Tables (ref) and (ref) in Supplementary Material (ref)) yielded slightly inferior, but qualitatively similar results.
We considered a total of $183,487$ mothers who delivered exactly two births during $1995$ and $2009$ in the Commonwealth of Pennsylvania, and relocated at their second deliveries so that their “excess-travel-time” IVs at two deliveries were different. We considered covariates that measured mothers' neighborhood circumstances including poverty rate, median income, etc, mothers' demographic information including race (white or not), age, years of education, etc, and variables related to delivery including gestational age in weeks and length of prenatal care in months, and eight congenital diseases. The “excess-travel-time" IVs in both stages were then dichotomized: $1$ if above the median and $0$ otherwise. Mothers' treatment choice and their babies' mortality status at the first delivery were included as covariates for studying the second delivery. We used the multiple imputation by chained equations method (buuren2010mice) implemented in the R package MICE to impute missing covariate data, and repeat analysis on $5$ imputed datasets.
We assume that high-level NICUs do no harm compared to low-level NICUs; therefore, all partial identification intervals in this application were estimated under the monotone treatment response (MTR) assumption (manski2003partial). We considered estimating an IV-optimal DTR minimizing the maximum risk at each delivery; see Section (ref), and explored the trade-off between minimizing the maximum risk and the cost/capacity constraint by adding a generic penalty to the value function. We performed weighted classification using a classification tree with maximum depth equal to $3$ so that the resulting DTR is interpretable. Figure (ref) plots three estimated DTRs corresponding to no penalty attached, a moderate penalty, and a large penalty attached to attending a high-level NICU. When there is no penalty attached, all mothers are assigned to high-level NICUs. As we increase the penalty, fewer mothers (albeit mothers who benefit most from attending a high-level NICU) are assigned to high-level NICUs. For instance, Figure (ref) corresponds to sending $68\%$ mothers to a high-level NICU at their first deliveries and $59.9\%$ at their second deliveries. Mothers who are assigned to high-level NICUs according to this DTR either belong to racial and ethnic minority groups or are older and have premature gestational age. Similarly, Figure (ref) plots a regime where less than $10\%$ of mothers are assigned to a high-level NICU. Mothers who are assigned to high-level NICUs according to this DTR belong to racial and ethnicity minority groups and have premature births. Our analysis here seems to suggest that in general race/ethnicity, age, and gestational age are the most significant effect modifiers. Gestational age has long been hypothesized as an effect modifier; see lorch2012differential, yang2014estimation, michael2020instrumental; more recently, yannekis2020differential found a differential effect between different race/ethnic groups. On the other hand, mother's age appears to be a new discovery that worth looking into. Overall, our method both complemented previous published results and generated new insights.
We systematically study the problem of estimating an dynamic treatment regime from retrospective observational data using a time-varying instrumental variable. We formulate the problem under a generic partial identification framework, derive a counterpart of the classical $Q$-learning and Bellman equation under partial identification, and use it as the basis for generalizing a notion of IV-optimality to the dynamic treatment regimes. One important variant of the developed framework is a strategy to improve upon a baseline dynamic treatment regime. As demonstrated via extensive simulations, IV-improved DTRs indeed have favorable performance compared to the baseline DTRs, including baseline DTRs that are optimal under the no unmeasured confounding assumption.
With the increasing availability of administrative databases that keep track of clinical data, it is tempting to estimate some useful, real-world-evidence-based dynamic treatment regimes from such retrospective data. To make any causal/treatment effect statements from non-RCT data, an instrumental variable analysis is often better-received by clinicians. Fortunately, many reasonably good instrumental variables are available, e.g., daily precipitation, geographic distances, service providers' preference, etc. Many of these IVs are intrinsically time-varying and could be leveraged to estimate a dynamic treatment regime using the framework proposed in this article.
In practice, to deliver a most useful policy intervention, it is important to take into account various practical constraints, e.g., those arising from limited facility capacity or increased cost. Our framework can be readily extended to incorporating various constraints.
We conclude this article by mentioning a few open problems. First, our analysis depends on the assumption of bounded concentration coefficients (Assumption (ref)). It will be interesting to see if this assumption can be relaxed by imposing additional structural assumptions. There are some recent process in the reinforcement learning literature (see, e..g, jiang2017contextual,sun2019model,du2021bilinear), but whether those structural assumptions can be adapted to the current setting remains a question. Meanwhile, our proofs bypasses case-by-case analyses of $Q$-learning by assuming the existence of the contrast estimating oracle (Assumption (ref)). It is an interesting future direction to conduct more fine-grained analyses and characterize the optimal rate of convergence of those contrast functions. Finally, the current article considers estimating DTRs from a historical dataset. It might be of interest to extend our framework to the online interactive setup such as the one considered in liao2021instrumental.
{ {0.2pt plus 0.3ex} } \appendixtitleon \appendixtitletocon \addtocontents{toc}{\setcounter{tocdepth}{2}}