# Weighted Q-learning for optimal dynamic treatment regimes with MNAR
coavriates

Jian Sun , Li Su , and Bo Fu School of Data Science, Fudan University,
Shanghai 200433, China. Email<EMAIL_ADDRESS>Biostatistics Unit,
University of Cambridge, Cambridge, CB2 0SR, UK. Email: li.su@mrc-
bsu.cam.ac.uk.School of Data Science, Fudan University, Shanghai 200433,
China. Email<EMAIL_ADDRESS>
###### Abstract

Dynamic treatment regimes (DTRs) formalize medical decision-making as a
sequence of rules for different stages, mapping patient-level information to
recommended treatments. In practice, estimating an optimal DTR using
observational data from electronic medical record (EMR) databases can be
complicated by covariates that are missing not at random (MNAR) due to
informative monitoring of patients. Since complete case analysis can result in
consistent estimation of outcome model parameters under the assumption of
outcome-independent missingness, Q-learning is a natural approach to
accommodating MNAR covariates. However, the backward induction algorithm used
in Q-learning can introduce challenges, as MNAR covariates at later stages can
result in MNAR pseudo-outcomes at earlier stages, leading to suboptimal DTRs,
even if the longitudinal outcome variables are fully observed. To address this
unique missing data problem in DTR settings, we propose two weighted
Q-learning approaches where inverse probability weights for missingness of the
pseudo-outcomes are obtained through estimating equations with valid
nonresponse instrumental variables or sensitivity analysis. Asymptotic
properties of the weighted Q-learning estimators are derived and the finite-
sample performance of the proposed methods is evaluated and compared with
alternative methods through extensive simulation studies. Using EMR data from
the Medical Information Mart for Intensive Care database, we apply the
proposed methods to investigate the optimal fluid strategy for sepsis patients
in intensive care units.

Keywords  Backward-induction induced missing pseudo-outcome; Future-
independent missingness; Missing not at random; Nonresponse instrumental
variable; Q-learning; Sensitivity analysis.

##  1 Introduction

Dynamic treatment regimes (DTRs) formalize medical decision-making as a
sequence of decision rules, each corresponding to a specific decision stage,
that map available information at patient-level to a recommended treatment
(Murphy, 2003). Optimal DTRs have been estimated and evaluated in various
biomedical areas, including cancer (Xu et al., 2019), HIV (Yu et al., 2019),
and sepsis (Yu and Zhang, 2022), to name just a few. In practice, estimating
optimal DTRs becomes challenging when there are missing data in covariates,
treatments, or outcomes.

There is a rich literature on identifying optimal DTRs with complete data.
Various methods have been developed, including Q-learning (Murphy, 2003;
Nahum-Shani et al., 2012), A-learning (Schulte et al., 2014), and their
numerous variants (Goldberg and Kosorok, 2012; Shi et al., 2018; Kosorok and
Laber, 2019). These methods rely on approximate dynamic programming methods.
Direct-optimization methods, such as simultaneous outcome weighted learning
(Zhao et al., 2015), and model-based planning via g-computation (Xu et al.,
2016) have also been employed. In contrast, limited research to date has been
devoted to accommodating missing data when estimating optimal DTRs. For
survival outcomes, a few methods have been proposed to estimate optimal DTRs
by addressing missingness due to right-censoring (Simoneau et al., 2020; Cho
et al., 2023). For handling missing data in covariates, treatments and
outcomes in general, Shortreed et al. (2014) proposed a time-ordered nested
conditional imputation strategy when estimating optimal DTRs in sequential
multiple assignment randomized trials. Dong et al. (2020) applied the
augmented inverse probability weighting method to Q-learning and a generalized
version of outcome-weighted learning when dealing with missing data due to
patients’ dropout. In simpler settings of estimating single-stage optimal
treatment regimes, Shen et al. (2023) further discussed the multiple
imputation method proposed by Shortreed et al. (2014) for estimating and
evaluating optimal individualized treatment rules that were not directly
observed in the design. Also for the single-stage setting, Huang and Zhou
(2020) investigated the performance of an augment inverse probability weighted
estimator in the direct-optimization framework when covariates are missing at
random. However, these existing methods all focus on data missing at random
and could potentially yield sub-optimal DTRs when data for covariates,
treatments or outcomes are missing not at random (MNAR) (Little and Rubin,
2014).

In this research, we aim to address the issue of MNAR covariates in the
estimation of optimal DTRs, which is prominent in settings where patients are
informatively monitored, e.g., in electronic medical records (EMR) databases
or disease clinics. Our research was motivated by a recent investigation on
optimal fluid strategy for sepsis patients in intensive care units (ICUs)
(Speth et al., 2022). Using EMR data from the Medical Information Mart for
Intensive Care (MIMIC-III) (Johnson et al., 2016), Speth et al. (2022)
estimated an optimal two-stage DTR for treating sepsis patients with fluid
resuscitation in ICUs by incorporating baseline covariates such as
demographics, time-varying covariates such as mechanical ventilation usage and
vasopressor administration, and the outcome measure of multi-organ failure at
the first stage. However, according to the recommendations of the ‘Surviving
Sepsis Campaign’, hemodynamic variables should also be taken into account when
reassessing the need for additional fluid administration following the initial
resuscitation (Rhodes et al., 2017). Unfortunately, within the MIMIC-III
database, certain hemodynamic variables, including blood pressure, respiratory
rate, and body temperature, are subject to missingness at the conclusion of
the initial resuscitation. Furthermore, as highlighted by Awad et al. (2017),
patients in the MIMIC-III database were informatively monitored since patients
with more severe conditions, who typically presented with more abnormal
physiological indicators than others, received more intensive monitoring.
Consequently, their measurements of hemodynamic variables were less likely to
be missing. In this case, the missingness probabilities of the hemodynamic
variables are directly related to their own values, indicating that the
missing data mechanism is MNAR. This characteristic could potentially lead to
sub-optimal estimation of DTRs if existing methods are employed.

The identification and estimation of optimal DTRs become much more challenging
when covariates are MNAR. In the literature on handling MNAR covariates, it
has been concluded that the parameters of an outcome regression model can be
consistently estimated using complete case analysis as long as the missingness
of the covariates is not directly related to the outcome (Yang et al., 2019).
Therefore, considering that Q-learning relies on consistently estimated
outcome models and is both flexible and straightforward to implement, it
offers a natural approach for accommodating MNAR covariates when estimating
optimal DTRs. However, when implementing the step-wise outcome regression
algorithm in Q-learning, challenges arise not only from the missingness of the
covariates but also from the missingness of the _pseudo-outcomes_ , even if
outcome variables are fully observed. This is because the Q-learning algorithm
computes the pseudo-outcome in the preceding stage using covariates from the
subsequent stage. As a result, if there are MNAR covariates in the later
stages, it can lead to MNAR pseudo-outcomes in the earlier stages. We provide
a detailed illustration of this unique missing data problem with DTRs caused
by backward-induction-induced MNAR pseudo-outcomes in Section 2.3.

It is well known that with MNAR outcomes, outcome regression models are
generally not identifiable; see examples provided by Miao et al. (2016).
Inspired by the MNAR outcome literature, we propose two weighted Q-learning
approaches to tackling MNAR pseudo-outcomes in optimal DTR estimation. The
first approach is to employ estimating equations to estimate a working model
for missingness probability of the pseudo-outcome with valid nonresponse
instrumental variables (D’Haultfoeuille, 2010; Wang et al., 2014; Miao and
Tchetgen Tchetgen, 2016). The pseudo-outcome estimates based on observed data
are then weighted by the inverse of the estimated missingness probabilities in
the Q-learning algorithm. In scenarios where nonresponse instrumental
variables are not available, we specify a sensitivity parameter to quantify
the residual association between the missing pseudo-outcome at stage $t$ and
its conditional missingness probability given observed information up to stage
$t$. Thus the missingness probabilities of the pseudo-outcomes are functions
of the sensitivity parameter. We then perform a sensitivity analysis to
evaluate the robustness of the estimated optimal DTRs with respect to the pre-
specified sensitivity parameter (Kim and Yu, 2011; Wen and Seaman, 2018).

##  2 Method

###  2.1 Setting and notation

We consider an observational cohort with $n$ patients and finite $T$ treatment
stages. Let subscripts $i=1,\ldots,n$, $t=1,\ldots,T$ denote patients and
stages, respectively. Unless specified otherwise, we use capital letters to
denote random variables and lowercase letters to indicate specific
realizations of random variables. Complete data for the patients are assumed
independent and identically distributed and thus we suppress the patient-
specific subscript $i$. Let $A_{t}\in\\{-1,1\\}$ be the assigned treatment at
stage $t$ for $t=1,\ldots,T$. $\bm{X}_{t}\in\mathbb{R}^{p_{t}}$ is a $1\times
p_{t}$ vector of covariates measured before $A_{t}$, and $Y_{t}$ denotes the
longitudinal outcome measured after $A_{t}$ at stage $t$.

The final outcome of interest is a pre-specified function (e.g., sum or
maximum) of the longitudinal outcome in $T$ stages, $Y=f(Y_{1},\ldots,Y_{T})$.
Higher values of $Y$ are assumed to be better. We define
$\bm{H}_{1}=\bm{X}_{1}$, and
$\bm{H}_{t}=\left(\bm{H}_{t-1},A_{t-1},Y_{t-1},\bm{X}_{t}\right)$ for
$t=2,\ldots,T$. Thus, $\bm{H}_{t}$ represents the information available before
making the treatment decision at stage $t$. Let $\mathcal{H}_{t}$ denote the
support of $\bm{H}_{t}$. A DTR $\bm{d}$ consists of a set of decision rules
$(d_{1},\ldots,d_{T})$, where $d_{t}:\mathcal{H}_{t}\to\\{-1,1\\}$ is a
function that takes the observed history $\bm{h}_{t}$ as input and outputs a
treatment decision at stage $t$. An optimal DTR is the set of decision rules
$\left(d^{\mathrm{opt}}_{1},\ldots,d^{\mathrm{opt}}_{T}\right)$ that maximizes
the expectation of the final outcome. We formalize this definition using the
potential outcome framework.

For $t=1,\ldots,T$, let $Y_{t}^{*}(\bar{\bm{a}}_{t})$ denote the potential
longitudinal outcome at stage $t$ and $\bm{X}_{t+1}^{*}(\bar{\bm{a}}_{t})$
denote the potential covariates at stage $t+1$ if a patient, possibly contrary
to fact, had received the treatment sequence
$\bar{\bm{a}}_{t}=(a_{1},\ldots,a_{t})$ by stage $t$. We then define the set
of potential outcomes under $\bar{\bm{a}}_{t}$ as
$\bm{O}_{t}^{*}(\bar{\bm{a}}_{t})=\\{Y_{1}^{*}(a_{1}),\bm{X}_{2}^{*}(a_{1}),\ldots,Y_{t}^{*}(\bar{\bm{a}}_{t}),\bm{X}_{t+1}^{*}(\bar{\bm{a}}_{t})\\}$
for $t=1,\ldots,T-1$. The potential final outcome under a regime $\bm{d}$ is

$\displaystyle
Y^{*}(\bm{d})=\sum_{\bar{\bm{a}}_{T}}Y^{*}\left(\bar{\bm{a}}_{T}\right)\prod_{t=1}^{T}\mathbb{I}[d_{t}\left\\{(\bm{X}_{1},a_{1},Y_{1}^{*}(a_{1}),\ldots,a_{t-1},Y_{t-1}^{*}(\bar{\bm{a}}_{t-1}),\bm{X}_{t}^{*}(\bar{\bm{a}}_{t-1}))\right\\}=a_{t}]$

where
$Y^{*}\left(\bar{\bm{a}}_{T}\right)=f\left(Y_{1}^{*}\left(a_{1}\right),Y_{2}^{*}\left(\bar{\bm{a}}_{2}\right),\ldots,Y_{T}^{*}\left(\bar{\bm{a}}_{T}\right)\right)$.
Define the value of a regime $\bm{d}$ to be $V(\bm{d})=E[Y^{*}(\bm{d})]$. Let
$\mathcal{D}$ denote all feasible regimes. An optimal DTR,
$\bm{d}^{\mathrm{opt}}\in\mathcal{D}$, satisfies that
$V(\bm{d}^{\mathrm{opt}})\geq V(\bm{d})$ for all $\bm{d}\in\mathcal{D}$.

###  2.2 Q-learning with complete data

In this section, we briefly review the Q-learning algorithm for estimating
optimal DTRs with complete observational data. To estimate the value of
candidate regimes from observed data, it is necessary to express the value of
regimes solely in terms of observables rather than potential outcomes. This
becomes feasible under the following, now standard, causal inference
assumptions:

Assumption 1 (Consistency). $Y_{t}^{*}(\bar{\bm{a}}_{t})=Y_{t}$ for
$t=1,\ldots,T$ and $\bm{X}_{t+1}^{*}(\bar{\bm{a}}_{t})=\bm{X}_{t+1}$ for
$t=1,\ldots,T-1$ when $\bar{\bm{a}}_{t}$ are actually received.

Assumption 2 (Sequential ignorability).
$\\{\bm{O}_{T-1}^{*}(\bar{\bm{a}}_{T-1}),Y_{T}^{*}(\bar{\bm{a}}_{T}):\bar{\bm{a}}_{T}\in\bigotimes_{t=1}^{T}\\{-1,1\\}\\}\upmodels
A_{t}\mid\bm{H}_{t}$ for $t=1,\ldots,T$, where $\bigotimes$ denotes the
Cartesian product.

Assumption 3 (Positivity). $P(A_{t}=a_{t}\mid\bm{H}_{t})>0$ for
$a_{t}\in\\{-1,1\\}$, $\bm{H}_{t}\in\mathcal{H}_{t}$ and $t=1,\ldots,T$.

Under Assumptions 1-3, the value of a DTR can be identified from the observed
data. In particular, for the 2-stage scenario we have that

$\displaystyle
E[Y^{*}(\bm{d})]=E[E\\{E(Y\mid\bm{H}_{2},A_{2}=d_{2}(\bm{H}_{2}))\mid\bm{H}_{1},A_{1}=d_{1}(\bm{H}_{1})\\}]$

This repeated expectation form suggests that the optimal DTR can be computed
with a backward-induction procedure, which is expressed in terms of the
following recursively defined Q-functions in canonical Q-learning:

$\displaystyle Q_{T}\left(\bm{h}_{T},a_{T}\right)$
$\displaystyle=E\left(Y\mid\bm{H}_{T}=\bm{h}_{T},A_{T}=a_{T}\right),$
$\displaystyle Q_{t}\left(\bm{h}_{t},a_{t}\right)$
$\displaystyle=E\left\\{\max_{a_{t+1}\in\\{-1,1\\}}Q_{t+1}\left(\bm{H}_{t+1},a_{t+1}\right)\mid\bm{H}_{t}=\bm{h}_{t},A_{t}=a_{t}\right\\},$

for $t=T-1,\ldots,1$.

The true Q-functions are unknown and must be estimated from the data. Since
Q-functions represent conditional expectations, it is natural to estimate them
using regression models. Note that in these models, the response variable is
$\max_{a_{t+1}\in\\{-1,1\\}}Q_{t+1}\left(\bm{H}_{t+1},a_{t+1}\right)$ instead
of $Y_{t}$ at stage $t=1,\ldots,T-1$. These response variables are commonly
referred to as pseudo-outcomes, and we define $Y_{pse,T}=Y$ and
$Y_{pse,t}=\max_{a_{t+1}\in\\{-1,1\\}}Q_{t+1}\left(\bm{H}_{t+1},a_{t+1}\right)$
for $t=1,\ldots,T-1$. Subsequently,
$Q_{t}(\bm{h}_{t},a_{t})=E\left(Y_{pse,t}\mid\bm{h}_{t},a_{t}\right)$ for
$t=1,\ldots,T$. Since $a_{t}\in\\{-1,1\\}$ , for any specified model for
$Q_{t}(\bm{h}_{t},a_{t})$, we can decompose it as
$Q_{t}(\bm{h}_{t},a_{t};\bm{\theta}_{t})=q_{t,0}(\bm{h}_{t};\bm{\beta}_{t})+a_{t}q_{t,1}(\bm{h}_{t};\bm{\psi}_{t})$,
where $\bm{\theta}_{t}=(\bm{\beta}_{t},\bm{\psi}_{t})$,
$q_{t,0}(\bm{h}_{t})=\frac{Q_{t}(\bm{h}_{t},1)+Q_{t}(\bm{h}_{t},-1)}{2}$, and
$q_{t,1}(\bm{h}_{t})=\frac{Q_{t}(\bm{h}_{t},1)-Q_{t}(\bm{h}_{t},-1)}{2}$ for
$t=1\ldots,T$. As a result, the Q-function is divided into two components: (1)
a stage $t$ treatment-free component $q_{t,0}(\bm{h}_{t};\bm{\beta}_{t})$
which depends on (a subset of) the history $\bm{h}_{t}$ before stage $t$ but
not on the stage $t$ treatment $a_{t}$, and (2) a stage $t$ treatment effect
component $a_{t}q_{t,1}(\bm{h}_{t};\bm{\psi}_{t})$ which depends on (a
potentially different subset of) $\bm{h}_{t}$ and specifically includes the
main effect of treatment $a_{t}$ and its interactions with tailoring variables
(i.e., effect modifiers). The blip function, as defined in Chakraborty and
Moodie (2013), equals to
$a_{t}q_{t,1}(\bm{h}_{t};\bm{\psi}_{t})+q_{t,1}(\bm{h}_{t};\bm{\psi}_{t})$
under this setting. Therefore, we refer to $\bm{\psi}_{t}$ as the blip
function parameters. Using this modelling approach, we have
$d_{t}^{\mathrm{opt}}(\bm{h}_{t})=2\mathbb{I}\left(q_{t,1}(\bm{h}_{t};\bm{\psi}_{t})>0\right)-1$.
Under Assumptions 1-3, the estimator $\hat{\bm{\theta}}_{t}$ of
$\bm{\theta}_{t}$ can be recursively constructed by regressing $Y_{pse,t}$ on
$\bm{H}_{t}$ and $A_{t}$ for $t=T,\ldots,1$. Then, the Q-learning estimator of
$d^{\mathrm{opt}}_{t}$ can be obtained by
$\hat{d}_{t}^{\mathrm{opt}}\left(\bm{h}_{t}\right)=2\mathbb{I}\left(q_{t,1}(\bm{h}_{t};\hat{\bm{\psi}}_{t})>0\right)-1$
for $t=1,\ldots,T$.

###  2.3 Q-learning for DTRs with MNAR covariates: why complete case analysis
fails at earlier stages while working at the last stage?

We focus on the scenario that the treatment and longitudinal outcome,
$(A_{t},Y_{t})$, are fully observed while the covariates, $\bm{X}_{t}$,
contain missing values. This scenario is frequently encountered in practice
because treatments are typically well-documented and the longitudinal outcomes
of interest are usually key clinical disease activity indices which are
regularly measured. Moreover, the probability of a covariate being missing is
likely to depend on the actual value of the covariate itself due to
informative monitoring of patients. For instance, in the MIMIC-III data,
patients who displayed symptoms of abnormal heart rates such as chest pain
were more likely to have their heart rates recorded. Therefore, the absence of
heart rate measurements could indicate MNAR, since patients with normal heart
rates were less likely to be frequently monitored.

Let $R_{t}$ represent the missingness indicator for $\bm{X}_{t}$ such that
$R_{t}=1$ if all variables in $\bm{X}_{t}$ are fully observed and $R_{t}=0$ if
at least one component of $\bm{X}_{t}$ is missing. Define
$\bar{\bm{R}}_{t}=(R_{1},\ldots,R_{t})$ as the history of missingness
indicators up to stage $t$. We use $\bm{1}_{t}$ to denote a $t$-dimensional
vector in which all elements are equal to $1$ for $t=1,\ldots,T$. In ICU
settings, the effect of a treatment is commonly evaluated by the patient’s
longitudinal outcomes (e.g., physiological symptoms) observed several hours
after treatment administration, while the treatment decision for a patient is
based on the patient’s medical history. Because of the temporal ordering of
these covariates, it is plausible that the vector of missingness indicators
$\bar{\bm{R}_{t}}$ does not directly depend on the patient’s variables
measured after receiving the current stage treatment $A_{t}$. Thus we make the
following assumption:

Assumption 4 (Future-independent missingness).
$\bar{\bm{R}_{t}}\upmodels(Y_{t},\bm{X}_{t+1},R_{t+1},\ldots,\bm{X}_{T},R_{T},A_{T},Y_{T})\mid(\bm{H}_{t},A_{t})$
for $t=1,\ldots,T$.

Figure 1 is a directed acyclic graph illustrating Assumptions 2 and 4 in a
2-stage scenario. In this example, $\bm{H}_{1}=\bm{X}_{1}$,
$\bm{H}_{2}=(\bm{X}_{1},A_{1},Y_{1},\bm{X}_{2})$. $R_{1}$ and $R_{2}$ are
missingness indicators for $\bm{X}_{1}$ and $\bm{X}_{2}$, respectively.
Graphically, $A_{1}$ has no common parents with
($Y_{1},\bm{X}_{2},A_{2},Y_{2}$) except for $\bm{H}_{1}$, and $A_{2}$ has no
common parents with $Y_{2}$ except for $\bm{H}_{2}$, encoding Assumption 2.
Besides, $R_{1}$ is d-separated with $(Y_{1},\bm{X}_{2},R_{2},A_{2},Y_{2})$ by
$(\bm{H}_{1},A_{1})$, and $(R_{1},R_{2})$ is d-separated with $Y_{2}$ by
$(\bm{H_{2}},A_{2})$ (VanderWeele and Robins, 2007), encoding Assumption 4.

Figure 1: A direct acyclic graph illustrating Assumptions 2 and 4

In the single-stage scenario, Assumption 4 degenerates to the outcome-
independent missingness assumption (Ding and Geng, 2014; Yang et al., 2019),
which posits that the missingness probability of covariates or confounders is
conditionally independent of the outcome, given all variables except for the
outcome, i.e., $R_{1}\upmodels Y_{1}\mid(\bm{X}_{1},A_{1})$. Under this
assumption, we have
$P(Y_{1}\mid\bm{X}_{1},A_{1})=P(Y_{1}\mid\bm{X}_{1},A_{1},R_{1}=1)$, so that
we can obtain unbiased parameter estimates of the outcome model with complete
case analysis. Consequently, in the single-stage scenario, Q-learning based on
complete case analysis yields consistently optimal individualized treatment
rules. In this research, considering the plausibility of conditional
independence ensured by temporal ordering of the variables, we adapt the
outcome-independent missingness assumption to longitudinal settings with
multiple stages and make Assumption 4.

Given $\bm{H}_{t}$ and $A_{t}$, the pseudo-outcome $Y_{pse,t}$ only depends on
$Y_{t}$ and the variables measured after stage $t$. Therefore, under
Assumption 4, $Y_{pse,t}\upmodels\bar{\bm{R}}_{t}\mid(\bm{H}_{t},A_{t}$) and
we have
$P(Y_{pse,t}\mid\bm{H}_{t},A_{t})=P(Y_{pse,t}\mid\bm{H}_{t},A_{t},\bar{\bm{R}}_{t}=\bm{1}_{t})$
for $t=1,\ldots,T$. Thus, we can obtain unbiased estimates of
$Q_{t}\left(\bm{h}_{t},a_{t}\right)$ by utilizing data from patients with
complete data up to stage $t$ if there are no missing values in $Y_{pse,t}$.
For example, at stage $T$, $Y_{pse,T}=Y$ is fully observed. Thus we can use
complete case analysis to estimate $Q_{T}(\bm{h}_{T},a_{T})$.

However, conditional on $\bar{\bm{R}}_{t}=\bm{1}_{t}$ (i.e., for complete
cases up to stage $t$), although there is no missingness in $\bm{H}_{t}$, the
pseudo-outcome $Y_{pse,t}$ may contain missing values because $Y_{pse,t}$
depends on $\bm{X}_{t+1}$, which can be missing (i.e., $R_{t+1}=0$) when
$\bar{\bm{R}}_{t}=\bm{1}_{t}$. We illustrate this backward-induction-induced
missing pseudo-outcome problem in Figure 2 with a hypothetical example coming
from the same setting of our simulation study in Section 3. In this example,
all variables for Patient 3 at stage one are fully observed. Nevertheless,
$Y_{pse,1}$ for Patient 3 is unavailable because $X_{22}$, a time-varying
covariate required for computing $Y_{pse,1}$, is missing.

Figure 2: A hypothetical example illustrating the relationship between the
missing data pattern of covariates and that of pseudo-outcomes. Gray areas
represent missing information. The setting of the variables is the same as
that in Simulation 1 in Section 3.

Let $R_{pse,t}$ denote the missingness indicator of $Y_{pse,t}$ for
$t=1,\ldots,T$. Because the missingness probabilities of $\bm{X}_{t+1}$ depend
on their own values, and $Y_{pse,t}$ is computed based on $\bm{X}_{t+1}$, it
is likely that $Y_{pse,t}$ is missing not at random. This backward-induction-
induced MNAR pseudo-outcome problem would make identification and consistent
estimation of Q-functions impossible. A complete case analysis for Q-learning
would likely result in suboptimal DTR. In Section 3, we will use simulations
to demonstrate this phenomenon.

###  2.4 Weighted Q-learning for DTR with MNAR pseudo-outcomes: a nonresponse
instrument variable approach

We propose two approaches for inverse probability weighted Q-learning to
tackle the MNAR pseudo-outcome problem. The first involves utilizing valid
nonresponse instrumental variables. Specifically, we consider a semiparametric
model for $R_{pse,t}$:

Assumption 5 (semiparametric working model for missingness probability of the
pseudo-outcomes). For $t=1,\ldots,T-1$,

$\displaystyle
P(R_{pse,t}=1\mid\bm{h}_{t},a_{t},y_{pse,t},\bar{\bm{R}}_{t}=\bm{1}_{t})=\pi(\bm{u}_{t},y_{pse,t};\gamma_{t})=\frac{1}{1+\exp\\{s_{t}(\bm{u}_{t})+\gamma_{t}y_{pse,t}\\}},$
(1)

where $\bm{U}_{t}\subset(\bm{H}_{t},A_{t})$ is a proper subset of
$(\bm{H}_{t},A_{t})$, $P(Y_{pse,t}\mid\bm{H}_{t},A_{t})\neq
P(Y_{pse,t}\mid\bm{U}_{t})$, and $s_{t}(\cdot)$ is an unknown and unspecified
function of $\bm{u}_{t}$.

This semiparametric model assumption does not impose any parametric constraint
on the function of $\bm{u}_{t}$ in working model (1). Instead, it restricts
the parametric form of $y_{pse,t}$ and requires that there is no interaction
between $\bm{u}_{t}$ and $y_{pse,t}$ in working model (1). It is noteworthy
that the variables observed after stage $t$ may exhibit an association with
$R_{pse,t}$ even after conditioning on
($\bm{h}_{t},a_{t},y_{pse,t},\bar{\bm{R}}_{t}=\bm{1}_{t}$) and therefore model
(1) may not be the true data-generating model for $R_{pse,t}$. However,
weighting by the inverse of
$P(R_{pse,t}=1\mid\bm{h}_{t},a_{t},y_{pse,t},\bar{\bm{R}}_{t}=\bm{1}_{t})$
would suffice to balance the distribution of the variables between the
patients with and without missing $y_{pse,t}$ for the purpose of estimating
$Q_{t}(\bm{h}_{t},a_{t})$ and thus lead to unbiased estimates of the
Q-functions after weighting.

Assumption 5 also implies that $(\bm{H}_{t},A_{t})$ can be decomposed as
$(\bm{H}_{t},A_{t})=(\bm{U}_{t},\bm{Z}_{t})$ such that $\bm{Z}_{t}\nupmodels
Y_{pse,t}\mid\bm{U}_{t}$ and $\bm{Z}_{t}\upmodels
R_{pse,t}\mid\bm{U}_{t},Y_{pse,t}$ when $\bar{\bm{R}}_{t}=\bm{1}_{t}$. That
is, conditional on $\bm{U}_{t}$, $\bm{Z}_{t}$ is a predictor of the pseudo-
outcome $Y_{pse,t}$, but given $\bm{U}_{t}$ and $Y_{pse,t}$, $\bm{Z}_{t}$ is
conditionally independent of the missingness indicator $R_{pse,t}$.
$\bm{Z}_{t}$ is referred to as the nonresponse instrumental variables at stage
$t$. The nonresponse instrumental variables, sometimes known as the ‘shadow
variables’, are commonly employed in the literature for MNAR outcomes. It has
been demonstrated in previous studies of MNAR outcomes (Shao and Wang, 2016;
Miao and Tchetgen Tchetgen, 2016) that nonresponse instrumental variables
facilitate the identification of the working model (1) and the outcome model.
In certain scenarios, it may not be difficult to find valid nonresponse
instrumental variables. For example, gender and weight are associated with the
severity of sepsis in ICU patients (Kalani et al., 2020; Sakr et al., 2013).
It is plausible that they are independent of the missing data mechanism of
hemodynamic variables, given other pre-treatment covariates of the patients
(e.g., comorbidities) and the severity of sepsis represented by $Y_{pse,t}$.
These characteristics render gender and weight suitable candidates for
nonresponse instrumental variables.

If Assumptions 5 holds, to estimate the working model (1), we have to overcome
the challenge that $s_{t}(\cdot)$ is nonparametric. Following Shao and Wang
(2016), we consider the idea of profiling and obtain the following kernel
estimate of $s_{t}(\cdot)$ for a given $\gamma_{t}$.

$\displaystyle\exp\\{\hat{s}_{\gamma_{t},t}(\bm{u}_{t})\\}=\frac{\sum_{i=1}^{n}\left(1-r_{pse,t,i}\right)K_{c_{t}}\left(\bm{u}_{t}-\bm{u}_{t,i}\right)\mathbb{I}(\bar{\bm{R}}_{t,i}=\bm{1}_{t})}{\sum_{i=1}^{n}r_{pse,t,i}\exp\left(\gamma_{t}y_{pse,t,i}\right)K_{c_{t}}\left(\bm{u}_{t}-\bm{u}_{t,i}\right)\mathbb{I}(\bar{\bm{R}}_{t,i}=\bm{1}_{t})},$
(2)

where $K_{c_{t}}(\cdot)=c_{t}^{-1}K(\cdot/c_{t})$, with $K(\cdot)$ being a
symmetric kernel function and $c_{t}$ a bandwidth. For a fixed $\gamma_{t}$,
the working model (1) can be expressed as
$\hat{\pi}_{t}(\bm{u}_{t},y_{pse,t};\gamma_{t})=[1+\exp\\{\hat{s}_{\gamma_{t},t}(\bm{u}_{t})+\gamma_{t}y_{pse,t}\\}]^{-1}$.
Then, we can estimate $\gamma_{t}$ with the following estimating equation,

$\displaystyle\widehat{E}\left[\mathbb{I}(\bar{\bm{R}}_{t}=\bm{1}_{t})\bm{l}_{t}(\bm{z}_{t})\left\\{\frac{r_{pse,t}}{\hat{\pi}_{t}(\bm{u}_{t},y_{pse,t};\gamma_{t})}-1\right\\}\right]=\bm{0},$
(3)

where $\widehat{E}[\cdot]$ denotes the empirical mean,
$\bm{l}_{t}(\bm{z}_{t})=(l_{t,1}(\bm{z}_{t}),\ldots,l_{t,L}(\bm{z}_{t}))$ is a
user-specified differentiable vector function of the nonresponse instrumental
variables of dimension $L$, $L\geq 2$.

Note that in estimating equations (3), if $\bm{z}_{t}$ are constants, then the
above equations are not solvable because they are under-identified.
Nevertheless, when $L>2$ and the equations are over-identified, we can employ
the two-step generalized method of moments (GMM) (Hansen, 1982) to estimate
$\gamma_{t}$. Let
$b_{t,j}(r_{pse,t},\bm{u}_{t},y_{pse,t};\gamma_{t})=\mathbb{I}(\bar{\bm{R}}_{t}=\bm{1}_{t})l_{t,j}(\bm{z}_{t})\left(\frac{r_{pse,t}}{\hat{\pi}_{t}(\bm{u}_{t},y_{pse,t};\gamma_{t})}-1\right)$
for $j=1,\ldots,L$ and
$\bm{B}_{t}(\gamma_{t})=(\widehat{E}[b_{t,1}],\widehat{E}[b_{t,2}],\ldots,\widehat{E}[b_{t,L}])$.
Then, the first-step GMM estimator of $\gamma_{t}$ is
$\hat{\gamma}^{(1)}_{t}=\underset{\gamma_{t}}{\arg\min}\;\bm{B}_{t}(\gamma_{t})\bm{B}_{t}(\gamma_{t})^{\mathrm{T}}$.

Based on the first-step GMM estimator, we can build a more efficient estimator
of $\gamma_{t}$ with a weight matrix $\bm{W}_{t}$. Let $\bm{W}_{t}$ be the
inverse matrix of the $L\times L$ matrix with the $\left(j,j^{\prime}\right)$
element as
$\frac{1}{n}\sum_{i=1}^{n}b_{t,j}(r_{pse,t},\bm{u}_{t},y_{pse,t};\hat{\gamma}^{(1)}_{t})b_{t,j^{\prime}}(r_{pse,t},\bm{u}_{t},y_{pse,t};\hat{\gamma}^{(1)}_{t})$.
Then the second-step GMM estimator of $\gamma_{t}$ is
$\hat{\gamma}_{t}=\underset{\gamma_{t}}{\arg\min}\;\bm{B}_{t}(\gamma_{t})\bm{W}_{t}\bm{B}_{t}(\gamma_{t})^{\mathrm{T}}$.

With the estimated working model for the missingness probability, we can
employ the weighted least square method to obtain unbiased estimates of the
Q-function parameters. Subsequently, the optimal treatment at stage $t$ can be
determined based on the estimated parameters. We henceforth refer to this
approach as the estimating equation (EE) method.

###  2.5 Weighted Q-learning for DTRs with MNAR pseudo-outcomes: a sensitivity
analysis approach

In real data applications, it might be challenging to find valid nonresponse
instrument variables for every stage based on domain knowledge. In such
scenarios, we develop a sensitivity analysis approach for the conditional
missingness probability of $Y_{pse,t}$. Specifically, we slightly modify
Assumption 5 and introduce the following assumption.

Assumption 6 (Semiparametric working model for missingness probability without
nonresponse instrumental variables).

$\displaystyle
P(R_{pse,t}=1\mid\bm{h}_{t},a_{t},y_{pse,t},\bar{\bm{R}}_{t}=\bm{1}_{t})=\pi(\bm{h}_{t},a_{t},y_{pse,t};\gamma_{t}^{\prime})=\frac{1}{1+\exp\\{s^{\prime}_{t}(\bm{h}_{t},a_{t})+\gamma^{\prime}_{t}y_{pse,t}\\}},$
(4)

where $s^{\prime}_{t}(\cdot)$ is an unknown and unspecified function of
$(\bm{h}_{t},a_{t})$ and $\gamma^{\prime}_{t}$ is the pre-specified
sensitivity parameter. Note that here we include all historical variables in
(4).

The sensitivity parameter $\gamma^{\prime}_{t}$ quantifies the residual effect
of $Y_{pse,t}$ on $R_{pse,t}$ after adjusting for $\bm{h}_{t}$ and $a_{t}$ for
patients with complete records up to stage $t$. The working model (4) can be
identified with a fixed $\gamma^{\prime}_{t}$ (Kim and Yu, 2011). We can
obtain the kernel estimate of
$\exp\\{s^{\prime}_{\gamma^{\prime}_{t},t}(\bm{h}_{t})\\}$ in a similar manner
as in equation (2).

$\displaystyle\exp\\{\hat{s}^{\prime}_{\gamma^{\prime}_{t},t}(\bm{h}_{t})\\}=\frac{\sum_{i=1}^{n}\left(1-r_{pse,t,i}\right)K_{c^{\prime}_{t}}\left(\bm{h}_{t}-\bm{h}_{t,i}\right)\mathbb{I}(\bar{\bm{R}}_{t,i}=\bm{1}_{t})}{\sum_{i=1}^{n}r_{pse,t,i}\exp\left(\gamma^{\prime}_{t}y_{pse,t,i}\right)K_{c^{\prime}_{t}}\left(\bm{h}_{t}-\bm{h}_{t,i}\right)\mathbb{I}(\bar{\bm{R}}_{t,i}=\bm{1}_{t})},$
(5)

where
$K_{{c^{\prime}}_{t}}(\cdot)={c^{\prime}}_{t}^{-1}K(\cdot/c^{\prime}_{t})$,
with $K(\cdot)$ being a symmetric kernel function and $c^{\prime}_{t}$ a
bandwidth.

Subsequently, we can employ the inverse probability weighting method to obtain
unbiased estimates of the Q-functions and the optimal DTRs. If the working
model (4) is correctly specified (including the sensitivity parameter
$\gamma^{\prime}_{t}$) and certain regularity conditions for
$\exp\\{s^{\prime}_{\gamma^{\prime}_{t},t}(\bm{h}_{t})\\}$ are satisfied,
along with correctly specified Q-functions, the parameter estimators in the
Q-functions are consistent. As the true value of $\gamma^{\prime}_{t}$ is
unknown, we conduct a sensitivity analysis by varying the value of
$\gamma^{\prime}_{t}$ to assess the robustness of the Q-learning results with
respect to $\gamma^{\prime}_{t}$. We refer to this approach as the SA method.

###  2.6 Summary of the estimation steps

We summarize the estimation steps of the proposed weighted Q-learning approach
as follows:

  1. 1.

Compute $y_{pse,T}=f(y_{1},y_{2},\ldots,y_{T})$. Let
$\bm{M}_{T}(\bm{\theta}_{T})=\widehat{E}\left\\{\mathbb{I}(\bar{\bm{R}}_{T}=\bm{1}_{T})\frac{\partial\left[\left\\{Q_{T}(\bm{h}_{T},a_{T};\bm{\theta}_{T})-y_{pse,T}\right\\}^{2}\right]}{\partial\bm{\theta}_{T}}\right\\}$.
Then obtain $\hat{\bm{\theta}}_{T}$ by solving the estimating equation
$\bm{M}_{T}(\bm{\theta}_{T})=\bm{0}$

  2. 2.

For $t=T-1,\ldots,1$, repeat the following steps recursively:

    1. (a)

Construct stage $t$ pseudo-outcome as
$\hat{y}_{pse,t}=\max\limits_{a_{t+1}}\;Q_{t+1}(\bm{h}_{t+1},a_{t+1};\hat{\bm{\theta}}_{t+1})$,
denote the missingness indicator of $\hat{y}_{pse,t}$ as $\hat{r}_{pse,t}$.

    2. (b)

Estimate
$P(\hat{r}_{pse,t}\mid\bm{h}_{t},a_{t},\hat{y}_{pse,t},\bar{\bm{R}}_{t}=\bm{1}_{t})$
using the EE method or the SA method with corresponding working models for the
missingness probability of the estimated pseudo-outcome.

    3. (c)

Let
$\bm{M}_{t}(\bm{\theta}_{t})=\widehat{E}\left\\{\frac{\hat{r}_{pse,t}}{\hat{\pi}_{t}}\mathbb{I}(\bar{\bm{R}}_{t}=\bm{1}_{t})\frac{\partial\left[\left\\{Q_{t}(\bm{h}_{t},a_{t};\bm{\theta}_{t})-\hat{y}_{pse,t}\right\\}^{2}\right]}{\partial\bm{\theta}_{t}}\right\\}$.
Solve the estimating equation
$\bm{M}_{t}(\bm{\theta}_{t})=\bm{0}$ to obtain $\hat{\bm{\theta}}_{t}$.

  3. 3.

The optimal treatment regime at stage $t$ can be obtained by
$\hat{d}_{t}^{\mathrm{opt}}\left(\bm{h}_{t}\right)=\operatorname*{arg\,max}\limits_{a_{t}}Q_{t}(\bm{h}_{t},a_{t};\hat{\bm{\theta}}_{t})$
for $t=1,\ldots,T$.

Remark 1. Recall that, since $a_{t}\in\\{-1,1\\}$,
$Q_{t}(\bm{h}_{t},a_{t};\hat{\bm{\theta}}_{t})$ can be written as
$q_{t,0}(\bm{h}_{t};\hat{\bm{\beta}}_{t})+a_{t}q_{t,1}(\bm{h}_{t};\hat{\bm{\psi}}_{t})$.
Subsequently, $\hat{y}_{pse,t}$ can be calculated by
$q_{t+1,0}(\bm{h}_{t+1};\hat{\bm{\beta}}_{t+1})+|q_{t+1,1}(\bm{h}_{t+1};\hat{\bm{\psi}}_{t+1})|$
and $\hat{d}_{t}^{\mathrm{opt}}\left(\bm{h}_{t}\right)$ is determined by the
sign of $q_{t,1}(\bm{h}_{t};\hat{\bm{\psi}}_{t})$.

Remark 2. The pseudo-outcome $Y_{pse,t}$ is determined by the specified
parametric form of $Q_{t+1}(\cdot)$ and the true values of $\bm{\theta}_{t+1}$
for $t=1,\ldots,T-1$. Since the true value of $\bm{\theta}_{t+1}$ is unknown
in practice, we use $\hat{Y}_{pse,t}$ calculated by the parameter estimates
$\hat{\bm{\theta}}_{t+1}$ in estimation step (2). In this scenario,
$\hat{Y}_{pse,1}$ differs slightly from $Y_{pse,1}$, which leads
$\hat{\pi}_{t}$ to differ slightly from $\pi(\bm{u}_{t},y_{pse,t})$ and
$\pi(\bm{h}_{t},a_{t},y_{pse,t})$, when the EE and SA methods are used
respectively. Nevertheless, if all the models are correctly specified,
$\hat{\bm{\theta}}_{t+1}$ serves as a consistent estimator for
$\bm{\theta}_{t+1}$. Subsequently, $\hat{Y}_{pse,t}$ converges to $Y_{pse,t}$.
It follows that the estimated optimal DTR $\hat{\bm{d}}$ will also converge to
the true optimal DTR $\bm{d}$. We formalize these asymptotic properties in the
following section. In Section 3, we also conduct simulations to assess the
robustness of the proposed EE method against mild model misspecification,
e.g., when the chosen instrumental variable is weakly associated with
$R_{pse,t}$ after conditioning on $\bm{U}_{t}$ and $Y_{pse,t}$.

We provide an R Markdown tutorial for implementing the weighted Q-learning
with the EE method using simulated data in the Supplementary Materials.

###  2.7 Asymptotic properties

The consistency and asymptotic normality of the estimators in Section 2.6 can
be established under suitable regularity conditions for estimating equations
and an additional assumption that the optimal treatment is unique for all
subjects at all stages. Note that the parameters are estimated separately
using estimating equations recursively at each stage. Therefore, it is natural
to establish the asymptotic properties of the stage-specific parameters
recursively as well. For simplicity, we will focus on the two-stage setting,
but extensions to the general $T$-stage setting follow directly.

###### Theorem 1

In the 2-stage scenario, if Assumptions 1-4 hold, the Q-functions at stage 2
are Lipschitz continuous and correctly specified, and
$\bm{M}_{2}(\bm{\theta}_{2})=\bm{0}$ has a unique solution, then

$\displaystyle\hat{\bm{\theta}}_{2}-\bm{\theta}_{2}^{*}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N(\bm{0},\Sigma_{\bm{\theta}_{2}}),$
$\displaystyle\hat{d}_{2}^{\mathrm{opt}}(\bm{h}_{2})$
$\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}d_{2}^{\mathrm{opt}}(\bm{h}_{2})\text{ for all
}\bm{h}_{2}\in\mathcal{H}_{2},$ $\displaystyle\hat{Y}_{pse,1}$
$\displaystyle\stackrel{{\scriptstyle p}}{{\longrightarrow}}Y_{pse,1},$

where $\bm{\theta}_{2}^{*}$ stands or the true values of $\bm{\theta}_{2}$.
$\hat{\bm{\theta}}_{2}$, $\hat{d}_{2}^{\mathrm{opt}}(\cdot)$, and
$\hat{Y}_{pse,1}$ denote the estimators of $\bm{\theta}_{2}$,
$d_{2}^{\mathrm{opt}}(\cdot)$, and $Y_{pse,1}$, respectively. Notations
$\stackrel{{\scriptstyle d}}{{\longrightarrow}}$ and $\stackrel{{\scriptstyle
p}}{{\longrightarrow}}$ stand for converge in distribution and converge in
probability, respectively. The detailed form of $\Sigma_{\bm{\theta}_{2}}$ is
given in Web Appendix A.

###### Theorem 2

In the two-stage scenario, if Assumptions 1-5 hold, all models are Lipschitz
continuous and correctly specified, and the estimating equations for the
Q-functions and the estimating equation (3) each have a unique solution, then
under certain additional regularity conditions,

$\displaystyle\hat{\bm{\theta}}_{1,EE}-\bm{\theta}_{1}^{*}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N(\bm{0},\Sigma_{\bm{\theta}_{1}}),$
$\displaystyle\hat{d}_{1,EE}^{\mathrm{opt}}(\bm{h}_{1})$
$\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}d_{1}^{\mathrm{opt}}(\bm{h}_{1})\text{ for all
}\bm{h}_{1}\in\mathcal{H}_{1},$

where $\bm{\theta}_{1}^{*}$ stands or the true values of $\bm{\theta}_{1}$.
$\hat{\bm{\theta}}_{1,EE}$ and $\hat{d}_{1,EE}^{\mathrm{opt}}(\cdot)$ denote
the estimators of $\bm{\theta}_{1}$ and $d_{1}^{\mathrm{opt}}(\cdot)$ with the
working model (1) estimated by the EE method, respectively. The regularity
conditions and the detailed form of $\Sigma_{\bm{\theta}_{1}}$ are given in
Web Appendix B.

###### Theorem 3

In the two-stage scenario, if Assumptions 1-4 and 6 hold, all models are
Lipschitz continuous and correctly specified, the sensitivity parameter
$\gamma^{\prime}$ is correctly specified, and the estimating equations for the
Q-functions each have a unique solution, then under certain additional
regularity conditions,

$\displaystyle\hat{\bm{\theta}}_{1,SA}-\bm{\theta}_{1}^{*}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N(\bm{0},\Sigma^{\prime}_{\bm{\theta}_{1}}),$
$\displaystyle\hat{d}_{1,SA}^{\mathrm{opt}}(\bm{h}_{1})$
$\displaystyle\stackrel{{\scriptstyle
p}}{{\longrightarrow}}d_{1}^{\mathrm{opt}}(\bm{h}_{1})\text{ for all
}\bm{h}_{1}\in\mathcal{H}_{1},$

where $\bm{\theta}_{1}^{*}$ stands or the true values of $\bm{\theta}_{1}$.
$\hat{\bm{\theta}}_{1,SA}$ and $\hat{d}_{1,SA}^{\mathrm{opt}}(\cdot)$ denote
the estimators of $\bm{\theta}_{1}$ and $d_{1}^{\mathrm{opt}}(\cdot)$ with the
working model (4) estimated by the SA method, respectively. The regularity
conditions and the detailed form of $\Sigma^{\prime}_{\bm{\theta}_{1}}$ are
given in Web Appendix C.

Note that the form of the asymptotic variance of the parameters at stage 1 is
complicated. Therefore, we suggest using the bootstrap method to estimate it
in practice.

##  3 Simulation

We conducted two simulation studies to evaluate the finite-sample performance
of the WQ-EE method. In Simulation 1, we focused on the scenario where
Assumptions 1-5 were satisfied. In Simulation 2, we evaluated the robustness
of the WQ-EE method against the violations of Assumption 5 when (1) the chosen
instrumental variable was weakly associated with $R_{pse,t}$ and (2) there
were interactions between $\bm{u}_{t}$ and $y_{pse,t}$ in the true missingness
probability model.

For comparison, we also assessed three other methods for handling partially
missing covariate data in Q-learning: a naive method which ignored partially
observed covariates (‘naive’), complete case analysis (CC), and multiple
imputation (MI) (Little and Rubin, 2014). Q-learning when all the covariates
were fully observed (All) was used as the benchmark. In the naive method, we
omitted all terms containing the partially observed covariates in the
specified Q-functions. Consequently, the Q-functions were misspecified in this
scenario as Assumption 2 was violated. In the CC method, we excluded all
patients whose covariates were subject to missingness in estimation. In the MI
method, we first imputed missing covariates using predictive mean matching
multiple times, which was implemented through the R package MICE.
Subsequently, we applied the standard Q-learning algorithm to the imputed data
sets. Finally, we calculated the average of the estimated parameters from the
imputed datasets to obtain the final estimates of the Q-function parameters.

All methods were evaluated in two aspects. Firstly, we compared the
performances of the methods for estimating the blip function parameters
$\bm{\psi}_{t}$. Secondly, we compared their abilities to identify the true
optimal dynamic treatment regime (DTR) and the values of the estimated optimal
DTRs.

###  3.1 Simulation 1: comparative performance when Assumptions 1-5 were
satisfied

We simulated data from an observational study with two stages of intervention.
There were two covariates $X_{t,1}$ and $X_{t,2}$ at stage $t=1,2$. Among
these covariates, only $X_{1,2}$ and $X_{2,2}$ were partially observed, and
their missingness indicators were $R_{1}$ and $R_{2}$, respectively. Denote
$\text{expit}(b)=\exp(b)/(1+\exp(b))$ for $b$ in $\mathcal{R}$. The variables
were generated as follows:

$\displaystyle\left(\begin{matrix}X_{1,1}\\\ X_{2,1}\end{matrix}\right)\sim
N\left(\left(\begin{matrix}0\\\
0\end{matrix}\right),\left(\begin{matrix}1&0.5\\\
0.5&1\end{matrix}\right)\right);$ $\displaystyle
X_{1,2}\sim\text{Uniform}(0,2);R_{1}\sim\text{Bernoulli}([1+\exp(-3+X_{1,2})]^{-1});$
$\displaystyle A_{1}\sim
2\times\text{Bernoulli}(\text{expit}(-1+X_{1,1}+X_{1,2}-R_{1}))-1;$
$\displaystyle
Y_{1}=A_{1}(\psi_{10}+2+\psi_{11}X_{1,2})+\beta_{10}+\beta_{11}X_{1,1}+\beta_{12}X_{1,2}+\epsilon_{1};\epsilon_{1}\sim
N(0,3);$ $\displaystyle
X_{2,2}\sim\text{Uniform}(0,2);R_{2}\sim\text{Bernoulli}([1+\exp(-1+0.5X_{1,2}-Y_{1}-X_{2,1}-0.5X_{2,2})]^{-1});$
$\displaystyle A_{2}\sim
2\times\text{Bernoulli}(\text{expit}(-1-X_{1,1}-X_{1,2}+Y_{1}+X_{2,1}-R_{2}))-1;$
$\displaystyle
Y_{2}=A_{2}(\psi_{20}+\psi_{2A}A_{1}+\psi_{22}X_{2,2})+\beta_{2,A}A_{1}+\beta_{22}X_{2,1}+\beta_{23}X_{2,2}+\epsilon_{2};\epsilon_{2}\sim
N(0,1).$

The true values of the parameters
$\bm{\theta}_{1}=(\psi_{10},\psi_{11},\beta_{10},\beta_{11},\beta_{12})$ and
$\bm{\theta}_{2}=(\psi_{20},\psi_{2A},\psi_{22},\beta_{2,A},\beta_{22},\beta_{23})$
were $(1,-1,1.5,0.5,-0.5)$ and $(1,-1,1,-1,1,-0.5)$, respectively. We
generated $1000$ data sets with sample sizes of $500$ and $2000$.

We considered $Y=Y_{1}+Y_{2}$ as the final outcome. Under this setting, the
true optimal treatment rule at stage 2 was
$2\mathbb{I}(\psi_{20}+\psi_{2A}a_{1}+\psi_{22}x_{2,2}\geq 0)-1$. Because
$\bm{\theta}_{2}=(1,-1,1,-1,1,-0.5)$ and $1-a_{1}+x_{2,2}\geq 0$, the optimal
choice of $A_{2}$ was always 1. Then, the pseudo-outcome at stage $1$,
$Y_{pse,1}$, was equal to $1-2A_{1}+Y_{1}+X_{2,1}+0.5X_{2,2}$. Therefore,
$Q_{1}(x_{1,1},x_{1,2},a_{1})=1.5+\beta_{10}+a_{1}(\psi_{10}+\psi_{11}x_{1,2})+(0.5+\beta_{11})x_{1,1}-\beta_{1,2}x_{1,2}$.
And the optimal treatment rule at stage 1 was
$2\mathbb{I}(\psi_{10}+\psi_{11}X_{1,2}\geq 0)-1$. Note that the missingness
indicators of $X_{1,2}$ and $X_{2,2}$ were directly associated with their own
values, thus they were missing not at random (MNAR). $R_{1}$ only depended on
$X_{1,2}$, and $R_{2}$ only depended on $(X_{1,2},Y_{1},X_{2,1},X_{2,2})$, so
the future-independent missingness assumption held. The missingness indicator
for $Y_{pse,1}$ was the same as that for $X_{2,2}$, since
$Y_{pse,1}=1-2A_{1}+Y_{1}+X_{2,1}+0.5X_{2,2}$ and $(A_{1},Y_{1},X_{2,1})$ were
fully observed. Besides, since
$P(R_{2}=1\mid\bm{H}_{2})=\\{1+\exp(-1+0.5X_{1,2}-Y_{1}-X_{2,1}-0.5X_{2,2})\\}^{-1}$,
it follows that
$P(R_{pse,1}=1\mid\bm{H}_{1},A_{1},Y_{pse,1})=\left\\{1+\exp(0.5X_{1,2}-2A_{1}-Y_{pse,1})\right\\}^{-1}$,
which indicates that $Y_{pse,1}$ was MNAR and $X_{1,1}$ was conditional
independent of $R_{pse,1}$ given $(X_{1,2},A_{1},Y_{pse,1})$. Since $X_{1,1}$
was also associated with $Y_{pse,1}$ after conditioning on
$(X_{1,2},A_{1},Y_{pse,1})$, we employed $X_{1,1}$ as a nonresponse
instrumental variable.

Figure 3: Boxplots of blip function parameter estimates from the compared
methods in Simulation 1. In each boxplot, white and grey boxes are the results
for the sample sizes 500 and 2,000, respectively. The horizontal line marks
the true values of the blip function parameter estimates. Note that the naive
method did not generate estimates of $\psi_{11}$ and $\psi_{22}$ as
$\psi_{11}$ and $\psi_{22}$ are the coefficients of partially missing
$X_{1,2}$ and $X_{2,2}$ in the true treatment rules.

Figure 3 shows the empirical distributions of the estimates of
$\bm{\psi}_{t}=(\psi_{10},\psi_{11},\psi_{20},\psi_{2A},\psi_{22})$, which
directly influence the estimated optimal DTR. At stage 2, the proposed WQ-EE
method and the CC method yielded unbiased estimates of the blip function
parameters since the future-independent missingness assumption held and the
outcome at stage 2 was fully observed. In contrast, the naive and MI methods
led to non-negligible biases of the blip function parameters. At stage 1, the
estimators in the WQ-EE method remained unbiased while the estimators in the
naive, CC, and MI methods displayed relatively large biases, which were not
mitigated as the sample size increased. Note that the naive estimates of
$\psi_{11}$ and $\psi_{22}$ were unavailable, as we excluded the terms
$A_{1}X_{1,2}$ and $A_{2}X_{2,2}$ from the specified Q-functions when
implementing this method.

The proportion of patients whose optimal treatments in both stages were
correctly identified was high when employing the WQ-EE method, as demonstrated
in Web Table 1. With the smaller sample size ($n=500$), the WQ-EE method
achieved an average correct classification rate of $89.9\%$ across the 1000
simulated datasets, with a mean final outcome of $2.98$. In general, larger
sample sizes led to higher proportions of patients with their true optimal DTR
correctly identified. In contrast, the DTRs obtained from CC and MI methods
had much lower correct classification rates and the mean final outcomes under
the estimated DTRs were lower than that from the WQ-EE method.

Table 1: Performance of estimated DTRs in Simulation 1

| n=500 | n=2000
---|---|---
| Value | Opt% | Value | Opt%
All | 2.988 (0.057) | 0.928 (0.052) | 2.997 (0.028) | 0.964 (0.024)
naive | 2.501 (0.053) | 0.500 (0.022) | 2.500 (0.027) | 0.499 (0.011)
CC | 2.879 (0.128) | 0.773 (0.111) | 2.898 (0.059) | 0.777 (0.057)
MI | 2.832 (0.139) | 0.729 (0.108) | 2.849 (0.071) | 0.730 (0.058)
WQ-EE | 2.974 (0.069) | 0.899 (0.073) | 2.994 (0.029) | 0.949 (0.036)

  * •

Note: Value, the value of the estimated regime; Opt%, the percentage of
patients correctly classified to their true optimal treatments in both stages;
Std, the sample standard error.

###  3.2 Simulation 2: comparative performance when Assumption 5 is violated

To investigate the robustness of the proposed WQ-EE method when Assumption 5
was violated, we kept all other data generating mechanisms unchanged, while
altering the true model for $R_{pse,1}$ to
$P(R_{pse,1}=1)=\left\\{1+\exp(\gamma_{z}X_{1,1}+0.5X_{1,2}-2A_{1}-Y_{pse,1})\right\\}^{-1}$.
When $\gamma_{z}=0$, the data generating mechanism was the same as that in
Simulation 1 and Assumption 5 was satisfied. However, when $\gamma_{z}\neq 0$,
$X_{1,1}$ was associated with $R_{pse,1}$ conditional on
$(X_{1,2},A_{1},Y_{pse,1})$. Thus there were no nonresponse instrumental
variables. We set $\gamma_{z}\in\\{-0.4,-0.2,0.2,0.4\\}$ such that $X_{1,1}$
was weakly associated with $P(R_{pse,1}=1)$ conditional on $X_{1,2}$ and
$Y_{pse,1}$. We generated 1000 Monte Carlo samples with sample size $n=500$.

Since the future-independent missingness assumption still held in Simulation
2, the estimators for $(\psi_{20},\psi_{2A},\psi_{22})$ in the WQ-EE method
were still consistent. We report the performance of the estimators for
$(\psi_{10},\psi_{11})$ and the value and correct classification rate of the
estimated optimal DTRs. The results were calculated by averaging over 1000
simulated data sets and are shown in Web Table 2. We see that the proposed WQ-
EE method was robust when $X_{11}$ was weakly related to $R_{pse,1}$. Despite
the declining performance of the corresponding optimal DTR as the absolute
value of $\gamma_{z}$ increased, the WQ-EE method consistently outperformed
the other three methods with higher mean final outcomes and correct
classification rates.

Table 2: Performance of the WQ-EE method without valid nonresponse
instrumental variable in comparison with alternative methods in Simulation 2

| $\psi_{10}$ | $\psi_{11}$ | Value | Opt%
---|---|---|---|---
Method | Bias | Std | Bias | Std | Mean | Std | Mean | Std
$\gamma_{z}=-0.4$
naive | -1.005 | 0.168 | — | — | 2.498 | 0.054 | 0.499 | 0.022
CC | 0.068 | 0.188 | -0.447 | 0.277 | 2.881 | 0.131 | 0.773 | 0.113
MI | 0.166 | 0.169 | -0.588 | 0.242 | 2.838 | 0.144 | 0.738 | 0.116
WQ-EE | -0.042 | 0.237 | 0.113 | 0.318 | 2.971 | 0.078 | 0.896 | 0.080
$\gamma_{z}=-0.2$
naive | -1.008 | 0.166 | — | — | 2.498 | 0.054 | 0.500 | 0.023
CC | 0.062 | 0.184 | -0.443 | 0.266 | 2.885 | 0.126 | 0.781 | 0.113
MI | 0.161 | 0.165 | -0.591 | 0.231 | 2.838 | 0.142 | 0.733 | 0.113
WQ-EE | -0.043 | 0.243 | 0.090 | 0.315 | 2.972 | 0.081 | 0.898 | 0.076
$\gamma_{z}=0.2$
naive | -1.000 | 0.155 | — | — | 2.502 | 0.052 | 0.499 | 0.022
CC | 0.073 | 0.185 | -0.482 | 0.253 | 2.871 | 0.125 | 0.766 | 0.110
MI | 0.174 | 0.164 | -0.627 | 0.227 | 2.816 | 0.144 | 0.718 | 0.113
WQ-EE | -0.040 | 0.239 | 0.031 | 0.310 | 2.974 | 0.070 | 0.899 | 0.072
$\gamma_{z}=0.4$
naive | -1.007 | 0.151 | — | — | 2.501 | 0.053 | 0.499 | 0.023
CC | 0.088 | 0.185 | -0.495 | 0.274 | 2.862 | 0.132 | 0.755 | 0.114
MI | 0.183 | 0.163 | -0.639 | 0.238 | 2.808 | 0.146 | 0.705 | 0.114
WQ-EE | -0.020 | 0.233 | 0.053 | 0.321 | 2.968 | 0.074 | 0.895 | 0.078

  * •

Note: Value, the value of the estimated regime; Opt%, the percentage of
patients correctly classified to their true optimal treatments in both stages;
Std, the sample standard error.

We also considered scenarios where there were interactions between $y_{pse,t}$
and $\bm{u}_{t}$ in the true missingness probability model. The true model for
$R_{pse,1}$ was modified as
$P(R_{pse,1}=1)=\left[1+\exp\\{0.5X_{1,2}-2A_{1}-Y_{pse,1}+\gamma_{uy}(X_{1,2}-1)Y_{pse,1}\\}\right]^{-1}$,
while the data generating mechanism for other covariates was the same as that
in Simulation 1. The specific form for the interaction,
$(X_{1,2}-1)Y_{pse,1}$, was chosen to ensure that the missing data proportions
were largely similar between Simulations 1 and 2. We set
$\gamma_{uy}\in\\{-0.2,-0.1,0.1,0.2\\}$ by allowing weak interactions between
$X_{1,2}$ and $Y_{pse,1}$. We again generated 1000 simulated data sets with
sample size $n=500$. The results of the compared methods are shown in Web
Table 3. The WQ-EE method exhibited superior performance compared to the other
three methods despite the introduction of interaction terms involving
$Y_{pse,1}$, showcasing its robustness in the presence of mild model
misspecification.

Table 3: Performance of the WQ-EE method when interaction terms between
$y_{pse,t}$ and $\bm{u}_{t}$ were incorrectly omitted, in comparison with
alternative methods in Simulation 2

| $\psi_{10}$ | $\psi_{11}$ | Value | Opt%
---|---|---|---|---
Method | Bias | Std | Bias | Std | Est | Std | Est | Std
$\gamma_{uy}=-0.2$
naive | -1.011 | 0.160 | — | — | 2.500 | 0.053 | 0.500 | 0.022
CC | 0.124 | 0.184 | -0.514 | 0.269 | 2.863 | 0.139 | 0.760 | 0.115
MI | 0.215 | 0.162 | -0.656 | 0.233 | 2.808 | 0.152 | 0.712 | 0.115
WQ-EE | -0.067 | 0.239 | 0.075 | 0.324 | 2.972 | 0.072 | 0.898 | 0.073
$\gamma_{uy}=-0.1$
naive | -1.002 | 0.153 | — | — | 2.500 | 0.053 | 0.499 | 0.022
CC | 0.084 | 0.182 | -0.474 | 0.270 | 2.871 | 0.132 | 0.764 | 0.116
MI | 0.185 | 0.159 | -0.624 | 0.230 | 2.820 | 0.141 | 0.715 | 0.115
WQ-EE | -0.062 | 0.228 | 0.065 | 0.318 | 2.971 | 0.072 | 0.897 | 0.075
$\gamma_{uy}=0.1$
naive | -1.002 | 0.157 | — | — | 2.499 | 0.053 | 0.500 | 0.022
CC | 0.046 | 0.182 | -0.441 | 0.258 | 2.886 | 0.119 | 0.780 | 0.105
MI | 0.147 | 0.163 | -0.595 | 0.231 | 2.832 | 0.143 | 0.734 | 0.107
WQ-EE | -0.021 | 0.233 | 0.022 | 0.307 | 2.974 | 0.074 | 0.898 | 0.077
$\gamma_{uy}=0.2$
naive | -1.010 | 0.160 | — | — | 2.497 | 0.051 | 0.499 | 0.023
CC | 0.046 | 0.188 | -0.444 | 0.263 | 2.883 | 0.124 | 0.772 | 0.112
MI | 0.146 | 0.168 | -0.590 | 0.234 | 2.835 | 0.140 | 0.732 | 0.111
WQ-EE | 0.023 | 0.242 | -0.011 | 0.308 | 2.972 | 0.079 | 0.895 | 0.082

  * •

Note: Value, the value of the estimated regime; Opt%, the percentage of
patients correctly classified to their true optimal treatments in both stages;
Std, the sample standard error.

Overall, our simulation studies demonstrated that the proposed WQ-EE method
performed well under correct model specification and mild model
misspecification, in comparison with the naive, CC, and MI methods.

##  4 Application to the MIMIC-III data

We applied the proposed methods to the MIMIC-III data (Johnson et al., 2016)
for investigating the optimal 2-stage fluid strategy to treat sepsis patients
in ICUs.

Sepsis is a life-threatening syndrome caused by the body’s response to
infection, resulting in damage to its own tissues and organs. Timely and
effective fluid resuscitation is crucial for stabilizing tissue hypoperfusion
or septic shock in sepsis patients. The ‘Surviving Sepsis Campaign’ guideline
strongly recommends administrating at least 30 mL/kg of intravenous (IV) fluid
within the first 3 hours of ICU admission (Rhodes et al., 2017). To address
the lack of robust evidence for the optimal fluid resuscitation strategy in
the early hours of treatment, Speth et al. (2022) utilized MIMIC-III data to
assess the optimal DTR during the 0-3 hours and 3-24 hours following admission
to the medical ICU. Building upon Speth et al. (2022) and in line with the
recommendation of the ‘Surviving Sepsis Campaign’ that additional fluid
administration after initial resuscitation should be guided by reassessing the
patient’s hemodynamic status, we incorporated hemodynamic variables in
determining the optimal fluid resuscitation strategy during the 3-24 hours
period after ICU admission (stage 2).

We followed similar patient selection criteria as outlined by Speth et al.
(2022), focusing on adult septic patients admitted to the medical ICU after
initially presenting to the emergency department. The detailed cohort
eligibility criteria can be found in Web Appendix D. The treatments in both
stages, the baseline covariates, and the final outcome were defined according
to Speth et al. (2022). Specifically, the stage 1 and stage 2 treatments were
categorized as either fluid restrictive ($<$30 ml/kg) or fluid liberal ($\geq$
30 ml/kg) strategies in 0-3 hours and 3-24 hours post-ICU admission,
respectively. The baseline covariates included gender, age, weight, racial
groups, and Elixhauser comorbidity score. The outcome of interest was the
Sequential Organ Failure Assessment (SOFA) score evaluated at 24 hours post-
admission. The SOFA score is a clinical tool used to assess the severity of
organ dysfunction by assigning scores to various organ systems. A higher SOFA
score indicates a more severe impairment.

For the intermediate variables considered before stage 2 treatment, we
included mechanical ventilation and vasopressors within the first 3-hour
period, as well as the patients’ SOFA score evaluated at 3 hours post-
admission. Additionally, we also considered heart rate, blood pressure, SpO2,
respiratory rate, temperature at 3 hours post-admission and urine output
within 0-3 hours after admission. These hemodynamic variables contained
missing values. As discussed previously, these variables were likely to be
MNAR due to informative monitoring of the patients in the MIMIC-III database.
The missing proportion of each hemodynamic variable is provided in Web Table
3.

There were 973 patients in the selected cohort, among them $53.2\%$ were male
and $78.7\%$ were Caucasian. 45.5% of the patients received the fluid liberal
strategy in 0-3 hours post-admission while $53.5\%$ received the fluid liberal
strategy in 3-24 hours post-admission. Overall, 67.1% of patients had fully
observed covariates at stage 2, while the rest of the patients had at least
one covariate with missing value.

As discussed in Section 2.3, it is plausible that Assumption 4 holds, i.e.,
that the missing data mechanism of the hemodynamic variables at 3 hours post-
admission is unrelated to the patient’s SOFA score at 24 hours post-admission,
conditioning on the patient’s medical history and treatments up to stage 2.
Denote the SOFA score at 24 hours post-admission under the optimal stage 2
treatment as $Y_{pse,1}$. It is plausible that weight and gender were
independent of the missingness of hemodynamic variables given baseline
covariates and $Y_{pse,1}$. Since weight and gender were associated with
$Y_{pse,1}$, they were suitable candidates for nonresponse instrumental
variables. Additionally, the Elixhauser comorbidity score can be regarded as a
mismeasured proxy of $Y_{pse,1}$, which may also serve as a candidate for the
nonresponse instrumental variable, according to Miao and Tchetgen (2018). We
investigated the impact of using different sets of candidates for nonresponse
instrumental variables in Table 4. Besides, we applied CC, MI, and the naive
method that ignored hemodynamic variables for comparison. The improvement of
SOFA at 24 hours was then calculated as the difference between the means of
the observed and the expected SOFA scores at 24 hours under the estimated
optimal DTR. We performed 1000 bootstrap resamplings and applied the 5-fold
cross-validation procedure to each bootstrap sample to assess the performance
of the estimated optimal DTRs. More details of the estimation procedure are
given in Web Appendix D.

Table 4: Improvement of SOFA at 24 hours post-admission for sepsis patients in
the selected cohort of the MIMIC-III database under the estimated optimal
fluid strategies with different approaches to handling MNAR covariates.

| Mean | Bootstrap std | Bootstrap 95% CI
---|---|---|---
WQ-EE (ECS) | 1.294 | 0.163 | (0.971, 1.634)
WQ-EE (weight) | 1.439 | 0.156 | (1.141, 1.760)
WQ-EE (ECS, weight) | 1.334 | 0.157 | (1.032, 1.668)
WQ-EE (gender) | 1.428 | 0.181 | (1.113, 1.775)
WQ-EE (gender, ECS) | 1.311 | 0.167 | (1.052, 1.716)
WQ-EE (gender, weight) | 1.462 | 0.165 | (1.167, 1.799)
WQ-EE (gender, weight, ECS) | 1.339 | 0.170 | (1.044, 1.692)
CC | 1.193 | 0.148 | (0.924, 1.476)
MI | 1.093 | 0.110 | (0.884, 1.316)
naive | 0.911 | 0.102 | (0.785, 1.127)

  * •

Note: Mean, the mean of bootstrap point estimates under 5-fold cross-
validation; Bootstrap Std, the sample standard error of 1000 cross-validation
estimates with bootstrap samples; Bootstrap 95% CI, 95% confidence interval
constructed by bootstrap percentiles; ECS, Elixhauser comorbidity score; WQ-EE
($\cdot$), the WQ-EE method based on the corresponding nonresponse
instrumental variables listed in parenthesis.

Table 4 presents the bootstrap mean, standard error, and 95% confidence
interval of the improvement of SOFA score at 24 hours post-admission. Notably,
all three methods that considered hemodynamic variables outperformed the naive
method, which aligns with the recommendation from the Surviving Sepsis
Campaign. Furthermore, the estimated improvement of SOFA at 24 hours post-
admission based on the WQ-EE method were notably larger than those based on
the CC and MI methods. Importantly, while the WQ-EE method exhibited
robustness to the choice of the three candidate nonresponse instrumental
variables, the estimated improvement of SOFA at 24 hours post-admission would
be lower if the Elixhauser comorbidity score was included in the set of
nonrsponse instrumental variables. This could be attributed to the association
between the Elixhauser comorbidity score and the missingness of hemodynamic
variables. In Simulation 2, we have demonstrated that this association can
affect the performance of the WQ-EE method. Therefore, we recommend excluding
the Elixhauser comorbidity score from the set of nonresponse instrumental
variables. The remaining three combinations of nonresponse instrumental
variables showed similar means and standard errors for the estimated
improvement of the final outcome.

However, we cannot rule out that gender and weight could be weakly associated
with the missingness of hemodynamic variables, even after conditioning on
other covariates and the pseudo-outcome. Therefore, we also applied the
proposed WQ-SA method in case none of the aforementioned variables were valid
nonresponse instrument variables. We assumed $\gamma^{\prime}_{1}\leq 0$
because patients with higher SOFA scores were more severely ill, and they may
undergo more intensive monitoring of their hemodynamic status. We considered a
range of values for $\gamma^{\prime}_{1}$ within the interval $[-2,0]$. The
rationale for selecting this range was that when we included gender and weight
as nonresponse instrumental variables, the estimated value of $\gamma_{1}$ in
the working model (1) results in a bootstrap mean of $-0.630$ with a bootstrap
standard error of $0.469$. We believe that the coefficient of $Y_{pse,1}$
would not change considerably if gender and weight were included in the
working model (4) for the conditional missingness probability. Furthermore,
when $\gamma^{\prime}_{1}=-2$, it deviated by almost three times the standard
error from the bootstrap mean, allowing a wide range for the residual
association with $Y_{pse,1}$. Table 5 shows that when $\gamma^{\prime}_{1}=0$,
the estimated improvement in SOFA at 24 hours post-admission was nearly
identical to that of the CC method, which also assumed a future-independent
missingness assumption and ignored the relationship between $Y_{pse,1}$ and
$R_{pse,1}$. The estimated improvement in SOFA at 24 hours post-admission
increased rapidly as $\gamma^{\prime}_{1}$ changed from $0$ to $-1$. When
$\gamma^{\prime}_{1}$ changes from -1 to -2, the estimated improvement in SOFA
at 24 hours post-admission stabilized at around 1.5.

Table 5: The improvement of SOFA at 24 hours post-admission for sepsis
patients in the selected cohort of the MIMIC-III database under the estimated
optimal fluid strategies with different values of the sensitivity parameter
$\gamma_{1}^{\prime}$ in the WQ-SA method

$\gamma^{\prime}_{1}$ | Mean | Bootstrap std | Bootstrap 95% CI
---|---|---|---
0 | 1.185 | 0.151 | (0.894, 1.479)
-0.5 | 1.349 | 0.160 | (1.060, 1.680)
-1 | 1.461 | 0.177 | (1.130, 1.827)
-1.5 | 1.511 | 0.185 | (1.172, 1.891)
-2 | 1.559 | 0.209 | (1.198, 2.027)

  * •

Note: Mean, the mean of bootstrap point estimates under 5-fold cross-
validation; Bootstrap Std, the sample standard error of 1000 cross-validation
estimates with bootstrap samples; Bootstrap 95% CI, 95% confidence interval
constructed by bootstrap percentiles.

While the challenge of optimizing the management of septic patients is
intricate and multifaceted, we employed the proposed methods to investigate
whether incorporating partially missing hemodynamic variables could influence
the fluid administration after the initial resuscitation to enhance the
overall outcome of the septic patients. In line with the best practice
recommendations of the Surviving Sepsis Campaign, our analysis results
indicated that certain hemodynamic variables should be taken into account when
formulating the fluid resuscitation strategy during the 3-24 hours following
admission to the medical ICU because accounting for these MNAR covariates with
the proposed methods increased the estimated improvement of SOFA at 24 hours
post-admission.

##  5 Discussion

In this article, we proposed two weighted Q-learning approaches to estimating
optimal DTRs with MNAR covariates and consequent MNAR pseudo-outcomes. To the
best of our knowledge, this is the first attempt to estimate optimal DTRs
under these challenging MNAR conditions. Specifically, we utilized the future-
independent missingness assumption and developed a semi-parametric estimating
equation-based approach to estimating the conditional missingness probability
of the pseudo-outcomes with valid nonresponse instruments. After obtaining the
estimated confitional missingness probabilities, we applied inverse
probability weighting to estimate the Q-function parameters. When nonresponse
instrumental variables were not available, we specified a sensitivity
parameter in the working model for missingness probability and conducted a
sensitivity analysis to assess the robustness of the estimated optimal DTRs.
In our simulation studies, we focused on evaluating the performance of the WQ-
EE method as the performance of the WQ-SA method would rely on the correct
specification of the sensitivity parameter. Our simulation results showed that
the estimators from the WQ-EE method were unbiased when Assumptions 1-5 hold.
They were also robust to mild model misspecification when Assumption 5 was
violated, while the alternative methods led to sub-optimal DTRs when
covariates were MNAR. In the application to the MIMIC-III data, we showed that
incorporating hemodynamic variables into optimal fluid strategy estimation and
handling MNAR problems with the proposed methods can improve the SOFA score
for sepsis patients in ICUs.

There are several directions for future work. First, this article exclusively
addresses the issue of missing pseudo-outcomes due to MNAR covariates in
optimal DTR estimation with Q-learning. However, in observational studies or
EMR databases, patients may not progress through all stages of clinical
intervention due to experiencing a terminal event or being censored before the
end of the study. In such cases, it is desirable to integrate the terminal
event and longitudinal outcome measurements for estimating optimal DTRs while
accounting for the censoring and MNAR covariates. Second, the dimension of the
covariates is often high in practice, especially in multiple-stage settings.
In such scenarios, the kernel regression method is susceptible to the ‘curse
of dimensionality’ problem. To address this challenge, one may consider
dimension reduction techniques proposed in Tang et al. (2014). Third, the
approach proposed for handling MNAR covariates and pseudo-outcomes can be
extended to other dynamic programming methods for DTRs, such as A-learning.

## References

  * Awad et al. (2017) Awad, A., Bader-El-Den, M., McNicholas, J., and Briggs, J. (2017).  Early hospital mortality prediction of intensive care unit patients using an ensemble learning approach.  International journal of medical informatics 108, 185–195.
  * Chakraborty and Moodie (2013) Chakraborty, B. and Moodie, E. E. (2013).  Statistical methods for dynamic treatment regimes.  Springer-Verlag. doi 10, 4–1.
  * Cho et al. (2023) Cho, H., Holloway, S. T., Couper, D. J., and Kosorok, M. R. (2023).  Multi-stage optimal dynamic treatment regimes for survival outcomes with dependent censoring.  Biometrika 110, 395–410.
  * Ding and Geng (2014) Ding, P. and Geng, Z. (2014).  Identifiability of subgroup causal effects in randomized experiments with nonignorable missing covariates.  Statistics in Medicine 33, 1121–1133.
  * Dong et al. (2020) Dong, L., Laber, E., Goldberg, Y., Song, R., and Yang, S. (2020).  Ascertaining properties of weighting in the estimation of optimal treatment regimes under monotone missingness.  Statistics in Medicine 39, 3503–3520.
  * D’Haultfoeuille (2010) D’Haultfoeuille, X. (2010).  A new instrumental method for dealing with endogenous selection.  Journal of Econometrics 154, 1–15.
  * Goldberg and Kosorok (2012) Goldberg, Y. and Kosorok, M. R. (2012).  Q-learning with censored data.  Annals of statistics 40, 529\.
  * Hansen (1982) Hansen, L. P. (1982).  Large sample properties of generalized method of moments estimators.  Econometrica: Journal of the econometric society pages 1029–1054.
  * Huang and Zhou (2020) Huang, Y. and Zhou, X.-H. (2020).  Identification of the optimal treatment regimen in the presence of missing covariates.  Statistics in medicine 39, 353–368.
  * Johnson et al. (2016) Johnson, A. E., Pollard, T. J., Shen, L., Lehman, L.-w. H., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Anthony Celi, L., and Mark, R. G. (2016).  Mimic-iii, a freely accessible critical care database.  Scientific data 3, 1–9.
  * Kalani et al. (2020) Kalani, C., Venigalla, T., Bailey, J., Udeani, G., and Surani, S. (2020).  Sepsis patients in critical care units with obesity: is obesity protective?  Cureus 12,.
  * Kim and Yu (2011) Kim, J. K. and Yu, C. L. (2011).  A semiparametric estimation of mean functionals with nonignorable missing data.  Journal of the American Statistical Association 106, 157–165.
  * Kosorok and Laber (2019) Kosorok, M. R. and Laber, E. B. (2019).  Precision medicine.  Annual review of statistics and its application 6, 263–286.
  * Little and Rubin (2014) Little, R. J. and Rubin, D. B. (2014).  Statistical Analysis with Missing Data, volume 333.  John Wiley & Sons.
  * Miao et al. (2016) Miao, W., Ding, P., and Geng, Z. (2016).  Identifiability of normal and normal mixture models with nonignorable missing data.  Journal of the American Statistical Association 111, 1673–1683.
  * Miao and Tchetgen (2018) Miao, W. and Tchetgen, E. T. (2018).  Identification and inference with nonignorable missing covariate data.  Statistica Sinica 28, 2049\.
  * Miao and Tchetgen Tchetgen (2016) Miao, W. and Tchetgen Tchetgen, E. J. (2016).  On varieties of doubly robust estimators under missingness not at random with a shadow variable.  Biometrika 103, 475–482.
  * Murphy (2003) Murphy, S. A. (2003).  Optimal dynamic treatment regimes.  Journal of the Royal Statistical Society Series B: Statistical Methodology 65, 331–355.
  * Nahum-Shani et al. (2012) Nahum-Shani, I., Qian, M., Almirall, D., Pelham, W. E., Gnagy, B., Fabiano, G. A., Waxmonsky, J. G., Yu, J., and Murphy, S. A. (2012).  Q-learning: a data analysis method for constructing adaptive interventions.  Psychological methods 17, 478\.
  * Rhodes et al. (2017) Rhodes, A., Evans, L. E., Alhazzani, W., Levy, M. M., Antonelli, M., Ferrer, R., Kumar, A., Sevransky, J. E., Sprung, C. L., Nunnally, M. E., et al. (2017).  Surviving sepsis campaign: international guidelines for management of sepsis and septic shock: 2016.  Intensive care medicine 43, 304–377.
  * Sakr et al. (2013) Sakr, Y., Elia, C., Mascia, L., Barberis, B., Cardellino, S., Livigni, S., Fiore, G., Filippini, C., and Ranieri, V. M. (2013).  The influence of gender on the epidemiology of and outcome from severe sepsis.  Critical Care 17, 1–9.
  * Schulte et al. (2014) Schulte, P. J., Tsiatis, A. A., Laber, E. B., and Davidian, M. (2014).  Q-and a-learning methods for estimating optimal dynamic treatment regimes.  Statistical science: a review journal of the Institute of Mathematical Statistics 29, 640\.
  * Shao and Wang (2016) Shao, J. and Wang, L. (2016).  Semiparametric inverse propensity weighting for nonignorable missing data.  Biometrika 103, 175–187.
  * Shen et al. (2023) Shen, J., Hubbard, R. A., and Linn, K. A. (2023).  Estimation and evaluation of individualized treatment rules following multiple imputation.  Statistics in Medicine 42, 4236–4256.
  * Shi et al. (2018) Shi, C., Fan, A., Song, R., and Lu, W. (2018).  High-dimensional a-learning for optimal dynamic treatment regimes.  Annals of statistics 46, 925\.
  * Shortreed et al. (2014) Shortreed, S. M., Laber, E., Scott Stroup, T., Pineau, J., and Murphy, S. A. (2014).  A multiple imputation strategy for sequential multiple assignment randomized trials.  Statistics in medicine 33, 4202–4214.
  * Simoneau et al. (2020) Simoneau, G., Moodie, E. E., Nijjar, J. S., Platt, R. W., Investigators, S. E. R. A. I. C., et al. (2020).  Estimating optimal dynamic treatment regimes with survival outcomes.  Journal of the American Statistical Association 115, 1531–1539.
  * Speth et al. (2022) Speth, K. A., Elliott, M. R., Marquez, J. L., and Wang, L. (2022).  Penalized spline-involved tree-based (pensit) learning for estimating an optimal dynamic treatment regime using observational data.  Statistical Methods in Medical Research 31, 2338–2351.
  * Tang et al. (2014) Tang, N., Zhao, P., and Zhu, H. (2014).  Empirical likelihood for estimating equations with nonignorably missing data.  Statistica Sinica 24, 723\.
  * VanderWeele and Robins (2007) VanderWeele, T. J. and Robins, J. M. (2007).  Directed acyclic graphs, sufficient causes, and the properties of conditioning on a common effect.  American journal of epidemiology 166, 1096–1104.
  * Wang et al. (2014) Wang, S., Shao, J., and Kim, J. K. (2014).  An instrumental variable approach for identification and estimation with nonignorable nonresponse.  Statistica Sinica pages 1097–1116.
  * Wen and Seaman (2018) Wen, L. and Seaman, S. R. (2018).  Semi-parametric methods of handling missing data in mortal cohorts under non-ignorable missingness.  Biometrics 74, 1427–1437.
  * Xu et al. (2016) Xu, Y., Müller, P., Wahed, A. S., and Thall, P. F. (2016).  Bayesian nonparametric estimation for dynamic treatment regimes with sequential transition times.  Journal of the American Statistical Association 111, 921–950.
  * Xu et al. (2019) Xu, Y., Thall, P. F., Hua, W., and Andersson, B. S. (2019).  Bayesian non-parametric survival regression for optimizing precision dosing of intravenous busulfan in allogeneic stem cell transplantation.  Journal of the Royal Statistical Society Series C: Applied Statistics 68, 809–828.
  * Yang et al. (2019) Yang, S., Wang, L., and Ding, P. (2019).  Causal inference with confounders missing not at random.  Biometrika 106, 875–888.
  * Yu et al. (2019) Yu, C., Dong, Y., Liu, J., and Ren, G. (2019).  Incorporating causal factors into reinforcement learning for dynamic treatment regimes in hiv.  BMC medical informatics and decision making 19, 19–29.
  * Yu and Zhang (2022) Yu, H. and Zhang, Q. (2022).  A deep reinforcement computation model for sepsis treatment.  In 2022 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pages 3217–3222. IEEE.
  * Zhao et al. (2015) Zhao, Y.-Q., Zeng, D., Laber, E. B., and Kosorok, M. R. (2015).  New statistical learning methods for estimating optimal dynamic treatment regimes.  Journal of the American Statistical Association 110, 583–598.