# Dynamic treatment regime characterization via a value function surrogate
with an application to partial compliance

Nikki L. B. Freeman1,∗, Sydney E. Browder2, Katharine L. McGinigle3,
and Michael R. Kosorok1
1Department of Biostatistics, University of North Carolina, Chapel Hill, NC
2Department of Epidemiology, University of North Carolina, Chapel Hill, NC
3Department of Surgery, University of North Carolina, Chapel Hill, NC

###### Abstract

Precision medicine is a promising framework for generating evidence to improve
health and health care. Yet, a gap persists between the ever-growing number of
statistical precision medicine strategies for evidence generation and
implementation in real world clinical settings, and the strategies for closing
this gap will likely be context dependent. In this paper, we consider the
specific context of partial compliance to wound management among patients with
peripheral artery disease. Through the use of a Gaussian process surrogate for
the value function, we expand beyond the common precision medicine task of
learning an optimal dynamic treatment regime to characterization of classes of
dynamic treatment regimes and how those findings can be translated into
clinical contexts.

##  1 Introduction

In precision medicine, dynamic treatment regimes (DTRs), sequences of
treatment rules over time that recommend treatments to patients based on
observable, possibly time-evolving characteristics, are often the object of
primary statistical interest. Common precision medicine analytic goals include
inference for the value of a fixed regime, inference for comparisons between
two or more regimes, and estimation of optimal regimes. Given the promise of
optimal DTRs to improve treatment decisions and patient outcomes when
employed, the focus on DTRs is well warranted. In clinical practice, however,
DTRs and applications of precision medicine analytic treatment recommendation
algorithms have not been widely adopted.

Bridging the gap between precision medicine evidence and clinical decision
support is important, context dependent, and may require expanding the focus
of precision medicine analyses. For example, in the treatment of patients with
lower limb wounds related to peripheral artery disease (PAD), a chronic
vascular condition characterized by arterial blockages and poor limb
perfusion, vascular surgeons attempt to save patients’ limbs by recommending a
combination of wound care and surgery to improve blood flow. Often outpatient
wound management in which patients complete regular wound clinic appointments
over weeks to months is sufficient to heal the wound, and surgery to
revascularize the limb can be avoided. With the goal of healing wounds and
preserving limbs, determining which patients would benefit most from surgery
and which patients can achieve equal outcomes with wound management alone is
of great clinical interest. In the parlance of precision medicine, learning an
optimal DTR for this single stage decision is of great importance.

However, given that the optimal treatment recommendation may depend on
adherence to a rigorous outpatient schedule for wound management, embedded in
this question is the need to understand if compliance matters in the
revascularization decision-making calculus. If so, it is also important to
know how much compliance is needed to heal wounds without surgery and how it
relates to equivalent or near-equivalent clinical outcomes. These questions
fall outside the objectives that traditionally have been pursued in precision
medicine, yet they represent natural questions related to understanding and
actually using optimal DTRs in a real world context.

Strategies for dealing some of these issues in treatment effect estimation,
such as imperfect compliance to treatment assignment, have been proposed
including the use of instrumental variables (IVs) (Angrist et al., 1996) and,
more generally, principal stratification (Frangakis and Rubin, 2002). In
precision medicine, (Cui and Tchetgen Tchetgen, 2021) offered an application
of IVs to identify the value of a particular regime and to learn and optimal
DTR under binary noncompliance. Artman et al. (2020) took a principal
stratification approach to estimate the value of regimes within strata in the
setting of a Sequential Multiple Assignment Randomization Trial (SMART)
(Lavori and Dawson, 2000, 2004; Murphy and Bingham, 2009). Still, some of the
questions raised in the PAD setting have not been the focus of precision
medicine yet are important to the translation of precision medicine into the
clinical setting.

In this paper, we bring a new perspective to the optimal DTR learning problem
by focusing on tools for understanding DTRs and how to make them useful in
practice which is particularly useful when considering real world complexities
like partial compliance. Through the use of a Gaussian process surrogate for
the value function, we (1) demonstrate the feasibility of employing Bayesian
optimization for optimal policy search over a policy class with finite
parameters, (2) exploit surrogate response surface methodology to characterize
the approximate value of all policies in a finite-parameter policy space, and
(3) describe how these findings can be translated to clinical contexts. The
paper proceeds as follows: In Section 2 we set notation and formalize the
setting. In Section 3, we describe the general strategy for our approach. In
Section 4 we present simulation results for using Bayesian optimization for
finding optimal DTRs, and in Section 5 we present an application of our
approach to the non-healing wounds in PAD setting that motivates our work. We
conclude in Section 6 with a discussion of our strategy, limitations, and
future work.

##  2 Notation and setting

###  2.1 Notation for the observed data

We consider three time points $k=0,1,2$ and two possible treatment options
that are feasible to all patients $\mathcal{A}=\\{0,1\\}$. We let $A_{0}$ and
$A_{1}$ denote the treatment assignment at $k=0$ and $k=1$, respectively, and
we further assume that all patients are initially assigned to treatment $0$,
i.e. $A_{0}\equiv 0$ for all patients. In the non-healing PAD wound setting,
this is equivalent to giving patients an initial trial of non-surgical wound
management instead of immediately opting for surgical intervention. Let
$H_{0}\equiv X_{0}\in R^{p}$ denote the baseline covariates observed at the
time of treatment initiation, $C_{0}\in[0,1]$ denote compliance with the
initial treatment assignment, and $X_{1}\in R^{q}$ denote the updated
covariates that arise in the interim after the initial treatment is assigned.
At $k=1$, the accumulated patient history is denoted by
$H_{1}=(X_{0},A_{0},C_{0},X_{1})$, where $\mathcal{H}_{1}$ denotes the set of
all such histories, and is available for the treatment assignment decision
between at $k=1$. We let $C_{1}\in[0,1]$ denote the measured compliance with
$A_{1}$, and let $Y$ denote the outcome of interest that is measured at $k=2$.
We assume our data contains i.i.d. sequences
$(X_{0,i},A_{0,i},C_{0,i},X_{1,i},A_{1,i},C_{1,i},Y_{i})$ for patients
$i=1,\ldots,n$.

###  2.2 Potential outcomes and the precision medicine framework

Our first goal is to learn an optimal dynamic treatment regime for the single
decision of whether to continue with treatment $0$ or to switch to the
treatment $1$ at the first decision point where treatment $0$ may be a lower
intensity treatment with the potential for partial patient compliance and
treatment $1$ is a higher intensity treatment with perfect compliance, i.e.
one-sided partial compliance. In the setting we have thus far described, a
dynamic treatment regime (DTR) $d$ is a map from the accumulated observed
patient history at $k=1$ to a treatment option,
$d:\mathcal{H}_{1}\longrightarrow\mathcal{A}_{1}$. Our interest is in regimes
that are indexed by a finite parameter $\boldsymbol{\theta}$ which we denote
by
$\mathcal{D}=\\{d(\boldsymbol{\theta}):\boldsymbol{\theta}\in\mathbf{\Theta}\subset\mathbb{R}^{m},m<\infty\\}$.
Using the potential outcomes framework (Splawa-Neyman et al., 1990; Rubin,
1974), let $C_{1}(a)$ and $Y(a,C_{1}(a))$ denote the compliance level and
outcome we would observe, respectively, after assigning treatment level $a$ at
$k=1$. Similarly, let $C_{1}(d)$ and $Y(d,C_{1}(d))$ denote the compliance
level and outcome, respectively, we would observe after assignment treatment
according to $d$ at $k=1$. We define the value of a DTR $d$ as

$\mathcal{V}(d)=\mathbb{E}[Y(d,C_{1}(d))]$

where the expectation is with respect to the trajectory generated under policy
$d$. An optimal dynamic treatment regime $d^{opt}$ is a regime such that
$\mathcal{V}(d^{opt})\geq\mathcal{V}(d)$ for all $d\in\mathcal{D}$.

###  2.3 Identification under one-sided partial compliance

To identify the value of a DTR, which is defined in terms of potential
outcomes, from the observed data, we must take care given that compliance to
the assigned treatment is an intermediate variable on the causal pathway
between treatment assignment and the outcome. In causal inference, one
approach to handling post-treatment intermediate outcome variables is
principal stratification, an approach that classifies individuals based on
sets of the values of their potential intermediate variables (Frangakis and
Rubin, 2002). The canonical example in causal effect estimation is the case of
two-sided noncompliance with binary treatments where strata are defined by the
two potential compliance values for each individual. Letting $0$ and $1$
represent the two treatments and $C(0)$ and $C(1)$ denote the potential
compliance levels for assignment to $0$ and $1$, respectively, the principal
strata are $\mathcal{G}_{NC}=\\{(0,0),(0,1),(1,0),(1,1)\\}$. Sometimes these
principal strata are referred to as “never takers” $(0,0)$, “compliers”
$(0,1)$, “defiers” $(1,0)$, and “always takers” $(1,1)$. A benefit of
principal strata is that they are unaffected by treatment assignment and thus
can be viewed as a pretreatment variable, however principal strata present
their own difficulties since they are unobserved. Treatment effects are often
estimated conditional on a principal strata and are referred to as principal
causal effects, e.g., $\mathbb{E}[Y(1)-Y(0)|G]$; identification of principal
causal effects relies on structural or modeling assumptions, e.g. monotonicity
or the exclusion restriction among others, with the plausibility of such
assumptions depending on the specific problem context (Frangakis and Rubin,
2002; Ding and Lu, 2017; Imbens and Rubin, 2015).

In our setting, the treatment is binary and we assume one-sided partial
compliance (OSPC), that is $C_{1}(1)=1$ and $C_{1}(0)\in[0,1]$. Principal
strata in this setting are thus
$\mathcal{G}=\\{(C_{1}(0),1):C_{1}(0)\in[0,1]\\}$ which we can further simply
as $\\{C_{1}(0):C_{1}(0)\in[0,1]\\}$. This differs from the canonical
noncompliance setting with binary treatment in the structural assumption that
compliance is one-sided and that compliance may be partial rather than all-or-
nothing. Further, in causal effect estimation, principal causal effects are
often the target of estimation and provide an estimate of the treatment
effect, for example, among those who comply. In contrast, to learn rules for
optimal treatment recommendations, we are interested in rules that work across
the target population, not just a targeted subgroup such as the complier
subgroup. For that reason, using the principal stratification framework, we
write the value of a DTR as

$\displaystyle\mathcal{V}(d)=\int_{0}^{1}\mathbb{E}[Y(d)|C_{1}(0)=c_{1}]dP_{C_{1}(0)}(c_{1}).$
(1)

To identify (1) from the observed data, we assume that the Stable Unit
Treatment Value Assumption (SUTVA) (Rubin, 1980) holds. We further assume that
conditional on $H_{1}$, treatment assignment is independent of the potential
outcomes $\\{Y(0),Y(1),C_{1}(0)\\}$ which we refer to as the no unmeasured
confounders (NUC) assumption, and we assume positivity for each level of
treatment $P(A_{1}=a|H_{1}=h_{1})>0\textrm{ for all
}a\in\mathcal{A},h_{1}\in\mathcal{H}_{1}$. We further assume that $Y(1)\perp
C_{1}(0)|H_{1}$ which is often referred to as principal ignorability (PI) (Jo
and Stuart, 2009). In causal effect estimation, the PI assumption has
motivated the use of principal scores for estimation of principal causal
effects.

###### Theorem 1

For the binary treatment setting with OSPC, if the SUTVA, NUC, positivity and
PI assumptions hold and if the class
$\mathcal{P}_{\mathcal{H}_{1}}=\\{P(C_{1}(0)|H_{1}=h_{1}):h_{1}\in\mathcal{H}_{1}\\}$
is complete, then (1) is identifiable.

We note that in cases when the PI assumption is not plausible, it may be
possible to identify (1) using an auxillary variable $W$ such that $Y(a)\perp
W|C_{1}(0),H_{1}$ for $a=0,1$. We refer to this assumption as auxillary
independence (AI) (Jiang and Ding, 2021). In our setting, a candidate
auxillary variable might be compliance to initial treatment to wound
management $C_{0}$.

##  3 Approach

In the setting described and given a parametric policy class $\mathcal{D}$
with parameters $\boldsymbol{\theta}$, our objective is to optimize
$\mathcal{V}$ over $\boldsymbol{\theta}$ and to characterize the estimated
values of the policies in $\mathcal{D}$. In this section we propose a simple
strategy for estimating the value of a fixed DTR that incorporates partial
compliance information. We follow by proposing the use of a Gaussian process
surrogate model for the value function to aid in both the optimization goal
and the characterization goal.

###  3.1 Value estimation with partial compliance

We use a simple imputation-based estimator for (1). $C_{1}(0)$ is missing for
those assigned to $A=1$. To impute this missing value, we propose estimating
$C_{1}(0)\sim H_{1}$ from the class of complete distributions among those with
$A_{1}=0$. Under the assumption of principal ignorability, we can then impute
$C_{1}(0)$ from the posterior predictive distribution of $C_{1}(0)|H_{1}$. Now
that $C_{1}(0)$ is complete, we can estimate a model $Y|C_{1}(0),H_{1},A$. To
estimate the value of a fixed regime, we take an approach similar to
(Chakraborty and Moodie (2013), Chapter 6.2): For a fixed, finite-parameter
regime $d(\boldsymbol{\theta})$

  1. 1.

Draw a large number of observations from the sample with replacement.

  2. 2.

For each observation, if $C_{1}(0)$ is unobserved, draw $C_{1}(0)$ from the
posterior predictive distribution $C_{1}(0)|H_{1}$.

  3. 3.

Fix the treatment vector such that $A$ is consistent with the fixed regime of
interest $d$.

  4. 4.

For each observation, draw from the posterior predictive distribution of
$Y|C_{1}(0),H_{1},A$.

To get a draws from the distribution of mean counterfactual outcomes under
regime $d$, we (5) average the draws in step $(4)$ and repeat steps (2)-(5) a
large number of times. One then may use the mean or median of the mean
counterfactual outcomes under regime $d$ as an estimate for (1).

###  3.2 General approach for DTR class characterization

Because one of our goals is to characterize policies in a policy class by
their value, an important aspect of our approach is estimation of the value
function for a fixed policy, i.e. policy evaluation. A number of strategies
for policy evaluation have been proposed including inverse probability
weighting (Robins et al., 2000), augmented inverse probability weighting Zhang
et al. (2012), and G-computation (Robins, 1986). While policy evaluation is
important for some precision medicine questions, it is not immediately
practical for finding optimal policies among a large policy class nor for
characterizing the value of an entire policy class. For example, one could
achieve these goals using grid search over the policy class which is
impractical except for the simplest of policies and settings.

Our approach starts with noting that for a finite parameter policy class the
value function can be viewed as a function mapping from the policy parameter
space to a real number
$\mathcal{V}:\boldsymbol{\theta}\longrightarrow\mathbb{R}$. While the value
function is generally unknown, we conceptualize policy evaluation as a
strategy for sampling the value function at a given policy parameter value. We
then model $\mathcal{V}$ using a Gaussian process regression (GPR) model, that
is, we place a Gaussian process (GP) prior on $\mathcal{V}$ which we update by
sampling $\mathcal{V}$ to yield a posterior predictive distribution that can
be used to approximate the value at unsampled policy parameters.

A Gaussian process $W$ indexed by a set $T$ is a stochastic process such that
finite-dimensional distributions of $(W(t_{1}),\ldots,W(t_{m}))$, for all
$t_{1},\ldots,t_{m}\in T$, $m\in\mathbb{N}$ are jointly Gaussian distributed
and completely specified by mean function $\mu:T\longrightarrow\mathbb{R}$ and
covariance kernel $K:T\times T\longrightarrow\mathbb{R}$ with
$\mu(t)=\mathbb{E}[W(t)]$ and $K(s,t)=cov(W(s),W(t))$ Rasmussen and Williams
(2005); Ghosal and Van der Vaart (2017). As random functions, Gaussian
processes can serve as priors on function spaces and from this view can be
considered as a flexible method for nonparametric Bayesian inference directly
on a function space and is often referred to as Gaussian process regression
(GPR). By Mercer’s theorem, GPR can be viewed as Bayesian linear regression
with possibly infinitely many basis functions Rasmussen and Williams (2005);
Ghosal and Van der Vaart (2017).

We write the GPR model for $\mathcal{V}$ as

$\displaystyle\mathcal{V}(\boldsymbol{\theta})\sim\mathcal{GP}(\boldsymbol{0},K(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}))$

where $\mathcal{GP}$ denotes Gaussian process prior and we have taken the mean
function to be zero. Assuming that we have a strategy for evaluating policies
in $\mathcal{D}$, e.g. inverse-probability weighting, augmented inverse-
probability weighting, G-computation, we can select a small number of policies
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{n}$, where the selection
of points are by design, such as a space-filling design, to evaluate and yield
corresponding estimates of the value
$\widehat{\mathcal{V}}^{(1)},\ldots,\widehat{\mathcal{V}}^{(n)}$. The $n$
evaluation pairs
$\\{(\boldsymbol{\theta},\widehat{\mathcal{V}}^{(1)}),\ldots,(\boldsymbol{\theta},\widehat{\mathcal{V}}^{(n)})\\}$
can be used to update the prior and yield the predictive distribution for
predicting the value for unevaluated policies, i.e. policies that correspond
to
$\boldsymbol{\theta}\notin\\{\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{n}\\}$.
In particular, if we let $\widehat{\mathcal{V}}^{(\cdot)}$ denote
$(\widehat{\mathcal{V}}^{(1)},\ldots,\widehat{\mathcal{V}}^{(n)})^{\top}$,
$\bar{\boldsymbol{\theta}}$ denote
$(\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{n})^{\top}$,
$\boldsymbol{\theta}_{n+1}$ denote the parameter indexing a policy whose value
we wish to approximate based on our model, and
$\widehat{\mathcal{V}}^{(n+1)}_{\ast}$ as the value corresponding to
$\boldsymbol{\theta}^{(n+1)}$, then we can write the joint distribution of the
evaluation pairs and the value function at $\boldsymbol{\theta}^{(n+1)}$ as

$\displaystyle\begin{bmatrix}\widehat{\mathcal{V}}^{(\cdot)}\\\
\widehat{\mathcal{V}}^{(n+1)}_{\ast}\end{bmatrix}\sim\mathcal{N}\begin{pmatrix}\boldsymbol{0},\begin{bmatrix}K(\bar{\boldsymbol{\theta}},\bar{\boldsymbol{\theta}})+\widehat{\sigma}^{2}_{n}I&K(\bar{\boldsymbol{\theta}},\boldsymbol{\theta}^{(n+1)})\\\
K(\boldsymbol{\theta}^{(n+1)},\bar{\boldsymbol{\theta}})&K(\boldsymbol{\theta}^{(n+1)},\boldsymbol{\theta}^{(n+1)})\end{bmatrix}\end{pmatrix}$

where $\widehat{\sigma}_{n}^{2}$ is the estimated variance for the value
estimates. The conditional distribution of
$\widehat{\mathcal{V}}^{(n+1)}_{\ast}$ given $\bar{\boldsymbol{\theta}}$,
$\widehat{\mathcal{V}}^{(\cdot)}$, and $\boldsymbol{\theta}^{(n+1)}$ follows
from basic distribution theory for the multivariate Gaussian distribution and
is given by

$\displaystyle\widehat{\mathcal{V}}^{(n+1)}_{\ast}|\bar{\boldsymbol{\theta}},\widehat{\mathcal{V}}^{(\cdot)},\boldsymbol{\theta}^{(n+1)}$
$\displaystyle\sim\mathcal{N}(K(\boldsymbol{\theta}^{(n+1)},\bar{\boldsymbol{\theta}})[K(\bar{\boldsymbol{\theta}},\bar{\boldsymbol{\theta}})+\widehat{\sigma}^{2}_{n}I]^{-1}\widehat{\mathcal{V}}^{(\cdot)},\text{cov}(\widehat{\mathcal{V}}^{(n+1)}_{\ast})),\text{
where}$ (2) $\displaystyle\text{cov}(\widehat{\mathcal{V}}^{(n+1)}_{\ast})$
$\displaystyle=K(\boldsymbol{\theta}^{(n+1)},\boldsymbol{\theta}^{(n+1)})-K(\boldsymbol{\theta}^{(n+1)},\bar{\boldsymbol{\theta}})[K(\bar{\boldsymbol{\theta}},\bar{\boldsymbol{\theta}})+\widehat{\sigma}^{2}_{n}I]^{-1}K(\bar{\boldsymbol{\theta}},\boldsymbol{\theta}^{(n+1)}).$

We use the estimated GPR model in two ways, first to help us find the optimal
DTR and second, to help us characterize the policy class at large. We describe
these two functions for the estimated GPR model in the next two subsections.

###  3.3 Optimal DTR learning

Bayesian optimization is a machine learning approach to solving problems of
the form $\underset{x\in\mathcal{X}}{\max}f(x)$ where
$\mathcal{X}\in\mathbb{R}^{d}$ is the feasible set, typically a simple set
like a hyper-rectangle or $d-$dimensional simplex, and $f$ is the objective
function. Typically, Bayesian optimization methods employ a surrogate for the
expensive-to-evaluate function $f$ and define an acquisition function from the
surrogate to guide which points to evaluate next. Different surrogates and
different acquisition functions correspond to different Bayesian optimization
strategies, and the fundamental idea is that the acquisition function enables
computationally efficient search for optima over $\mathcal{X}$ (Frazier,
2018).

In the optimal DTR learning context, the objective function we want to
optimize is the value function, the GPR model serves as our surrogate for the
value function, and an acquisition function is used to guide which value of
$\boldsymbol{\theta}$ to evaluate next. Multiple acquisition functions and
approaches to Bayesian optimization have been proposed, each with their own
rationale and varying theoretical guarantees. One popular approach is the
Expectation Improvement (EI) algorithm (Jones et al., 1998) whose acquisition
function can be viewed as balancing exploration and exploitation. EI may be
implemented to learn an optimal DTR as follows:

  1. 1.

Select $b$ initial values of $\boldsymbol{\theta}\in\mathbf{\Theta}$ using a
space-filling design; label them as
$\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(b)}$ and the
corresponding policies they index as $d^{(1)},\ldots,d^{(b)}$. Evaluate the
value function under each of the policies $d^{(1)},\ldots,d^{(b)}$, and label
them as $\mathcal{V}^{d^{(1)}},\ldots,\mathcal{V}^{d^{(b)}}$.

  2. 2.

Fit a Gaussian process to the counterfactual draws of
$-\mathcal{V}^{d^{(1)}},\ldots,-\mathcal{V}^{d^{(b)}}$ over the policy
parameter space $\boldsymbol{\Theta}$

$\displaystyle-\mathcal{V}\sim\mathcal{GP}(m(\boldsymbol{\theta}),k(\boldsymbol{\theta},\boldsymbol{\theta})).$
(3)

Here, we negate the draws from the value because EI is traditionally written
as a minimization routine.

  3. 3.

Compute the expected improvement $EI$

$EI(\boldsymbol{\theta})=(f^{(b)}_{min}-\mu_{(b)}(\boldsymbol{\theta}))\Phi\left(\frac{f^{(n)}_{min}-\mu_{b}(\boldsymbol{\theta})}{\sigma_{b}(\boldsymbol{\theta})}\right)+\sigma_{b}(\boldsymbol{\theta})\phi\left(\frac{f^{(b)}_{min}-\mu_{b}(\boldsymbol{\theta})}{\sigma_{b}(\boldsymbol{\theta})}\right)$
(4)

where $\Phi(\cdot)$ is the cumulative distribution function for the standard
normal distribution, $\phi(\cdot)$ is the probability density function for the
standard normal distribution, $\mu_{(b)}(\cdot)$ is the predictive mean
function of the fitted Gaussian process, and $\sigma_{b}(\cdot)$ is the
predictive variance function of the fitted Gaussian process.

  4. 4.

Find the value of $\boldsymbol{\theta}$ that maximizes EI and evaluate the
value function for the associated policy.

  5. 5.

Let $b\leftarrow b+1$ and repeat Step 3 and Step 4 until maximum EI small or
the sampling budget is met.

The estimated optimal DTR is the policy indexed by
$\boldsymbol{\theta}^{\ast}=\max\\{\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(b)}\\}$
so that $d^{opt}=d(H_{1};\boldsymbol{\theta}^{\ast})$.

####  3.3.1 Policy class characterization

While optimization of the value function may yield one or more points in the
policy parameter space that are optimal, solely optimizing the value function
alone cannot tell us how quickly the value changes with a change in one or
more of the parameters. Moreover, in the case of value functions where the
optima occur in flat, non-peaky regions of the value function, it may be the
case that an entire set of decision rules yield nearly equivalent estimated
values and equivalent benefit in terms of clinical significance. For example,
using the predictive mean function (2) of the fitted Gaussian process, we can
quickly approximate the value of policies in the entire policy class and
thereby identify multiple optima and the points in policy parameter space
where benefits begin to drop off. With the ability to sample from the
posterior predictive distribution, a wide range of numerical and visual
summaries can be computed which may further aid in characterizing DTRs.

##  4 Simulation

We conducted simulation studies to examine the performance of our approach for
DTR class characterization, both learning an optimal DTR and for approximating
the value of DTRs in a given policy class which are the novel aspects of our
approach. For Simulation Setting 1, we drew $U\sim Uniform(0,1)$ and set
$X=wU$ where $w$ is a positive scale parameter. We drew $A\sim Bern(0.5)$, and
set $Y=\gamma_{0}+\gamma_{1}x+A[I(x_{1}<0.75\text{ or
}x_{1}>0.25)+0.5*x*(1-I(x_{1}<0.75\text{ or }x_{1}>0.25))]$. For Simulation
Setting 2, we drew $U$, $A$, and $X$ in the same way as Simulation Setting 1
and set $Y=\gamma_{0}+\gamma_{1}x_{1}+A*cos(2\pi x_{1})$. Simulation Setting 3
is the same as Simulation Settings 1 and 2 except for setting
$Y=\gamma_{0}+\gamma_{1}x+A*\cos(4x\pi)$. In each simulation setting, we
searched over the policy class
$\mathcal{D}=\\{d(x;\beta):d(x;\beta)=\mathbbm{1}(x<\beta_{1}\text{ or
}x>\beta_{2}),\text{ for }\beta_{1},\beta_{2}\in[0,1]\\}$.

For each simulation setting, we considered $w\in\\{0.75,1,1.25\\}$. The true
value of policies in $\mathcal{D}$ for each setting and each value of $w$ is
visualized by the contour plots of the value as a function of policy
parameters $\beta_{0}$ and $\beta_{1}$ in Panel A of Figures 1, 2, and 3. In
Simulation Setting 1, there are many optimal policies for each value of $w$.
In Simulation Setting 2, the optimal policies are in a narrower range of
$\beta_{0}$ and $\beta_{1}$ and the optimal policy region is more complex than
that of Simulation Setting 1. Simulation Setting 3 is characterized by
multiple peaks.

Each simulation study was characterized by the setting and policy class under
examination, the sample size $N=200,500,1000$, $w=0.75,1,1.25$, and the policy
evaluation method employed: IPW, stabilized IPW, G-computation (gcomp), or
AIPWE. For the G-computation and AIPWE value estimators, we misspecified the
conditional mean models; for IPW, sIPW, and AIPWE we used the true sampling
weights. For each study we computed 1000 simulation runs. For each run we
initially evaluated 50 random policies to fit a GPR model for the value
function and then used EI to evaluate 50 more policies. The estimated optimal
policy is the policy, among those evaluated, with the highest value. To
characterize the policy class, we fixed a fine grid over the policy parameter
space and used the posterior predictive mean function of the GPR model (Matern
kernel + white noise) to approximate the value of policies on the grid. To
assess the performance of the GPR and EI approach to learn an optimal policy,
we computed the mean square error (MSE) between the estimated value and the
oracle value over the 1000 simulation runs; full results are given in the
appendix and visualized in Panels B, C and D of Figures 1, 2, and 3. To assess
the performance of the GPR to approximate the policy class, we computed the
$L^{1}$ norm and $L^{2}$ norm between the posterior predictive mean of the
value function GP surrogate and the true value at each point on the
aforementioned fine grid; full results are included in the appendix.

In Panels B, C, and D of Figures 1, 2, 3 the background contours correspond to
the true value of the policy indexed by $\beta_{0}$ and $\beta_{1}$ and each
open point represents the estimated value of the optimal DTR from a single
simulation run. By visual inspection, in Simulation Setting 1, each value
estimation method yielded optimal policies in the region of the policy
parameter space that corresponds to the highest value, meaning that optimal
policies were found (MSE for each case $\leq 0.007$). In Simulation Settings 2
and 3, the grossly misspecified conditional mean model for the G-computation
value estimator yielded sub-optimal policies as well as in the case of AIPWE
in Simulation Setting 2 for $w=0.75$. This is not surprising since our
strategy hinges on having a good estimator for the value. IPW and stabilized
IPW perform well across settings, with estimated optimal policies
concentrating on the regions of the policy parameter space with the highest
value.

Figure 1: Simulation Setting 1 Figure 2: Simulation Setting 2 Figure 3:
Simulation Setting 3

One particularly promising aspect of using a GPR surrogate for the value
function is the ability to characterize the values for policies in an entire
policy class, thereby finding cases where there are multiple optimal or near-
optimal policies that would yield similarly clinically significant results. We
demonstrate this with the first simulation run for Simulation Setting 3 with
$w=1$ and $N=200$, a setting with 3 distinct sub-classes of policies that
yield the highest values. In Figure 4, we plot as contours the predicted mean
values across the policy parameter space after training a GPR model (Matern
kernel + white noise) with the evaluation points used for optimization. The
red star represents the estimated optimal policy. In this case, the estimated
GPR surrogate identifies the regions of the policy parameter space that yields
policies with the highest values where the distance between the true value
surface and surrogate differs by 0.04 in the $L^{2}$ sense.

Figure 4: Example GPR surface from Simulation Setting 3 (sIPW with $N=200$,
$w=1$)

##  5 Dynamic treatment regimes for partially compliant ambulatory wound care
patients

###  5.1 Setting and data

PAD is narrowing of the arteries outside of the heart and brain. Because
diseased arteries result in reduced blood flow to the lower limbs, many
patients develop wounds that are slow to heal, painful, and at risk of
becoming infected or necrotic. PAD patients are complex with multiple
comorbidities such as diabetes, coronary artery disease, cerebrovascular
disease, and hypertension, and limb loss rates (amputation) and mortality
rates are high. For patients with PAD with non-healing wounds, treatment may
consist of outpatient wound management, resvascularization surgery such as a
bypass or an endovascular procedure, or a combination of both. Given the
multimorbidity of this patient population, non-surgical treatment of non-
healing wounds via wound management may be preferable to surgery when it is
possible to heal wounds with wound care alone (Browder et al., 2022), however,
wound management requires patients to attend scheduled appointments to receive
treatment until the wound is healed. Thus, it is useful to know (1) given a
trial of wound management, which patients would benefit from more intensive
treatment (surgery) and which patients would likely benefit from continued
wound management, and (2) how compliance to a trial of wound management should
drive the decision to switch to a more intensive strategy to heal patients’
wounds.

We analyzed a single-institution cohort of patients with PAD that were
referred to outpatient wound care. We extracted patient characteristics,
including age, sex, race, comorbidities, number of wounds, wound location, and
initial wound size from the electronic health records of patients with PAD
that were seen at UNC hospitals between 2013 and 2018. We linked the EHR data
to administrative scheduling data from a wound management clinic in the UNC
system. Patients were considered compliant with wound management if they
completed their scheduled appointments and non-compliant if they cancelled or
missed their scheduled appointments. We considered the first month following
the first wound management appointment as a trial of wound management and
compliance during that period as baseline compliance. For those with wounds
that remained unhealed at 1 month after initiating wound management and with
limbs that had not been revascularized during the first month of treatment, we
sought to learn an optimal dynamic treatment of the form “Revascularize along
with continued wound care if baseline compliance is less than or equal to
$\theta_{1}$ or if the initial wound size of the largest wound a patient has
is greater than $\theta_{2}$” that optimizes the outcome of being healed,
alive, and with an intact limb (not amputated) at 6 months after the
initiation of wound management. If we let $C_{0}$ denote the proportion of
completed wound care visits in the first month of wound management and $W_{0}$
denote initial wound size, then we can write the class of policies under
examination as
$\mathcal{D}=\\{d(C_{0},W_{0};\theta_{1},\theta_{2})=\mathbbm{1}(C_{0}<\theta_{1}\text{
or }W_{0}>\theta_{1}):\theta_{1}\in[0,1],\theta_{2}\in(0,100]\\}$ where
$d(C_{0}=c_{0},W_{0}=w_{0};\theta_{1},\theta_{2})=1$ recommends
revascularization and wound care for a particular baseline compliance $c_{0}$
and intial wound size $w_{0}$ and
$d(C_{0}=c_{0},W_{0}=w_{0};\theta_{1},\theta_{2})=0$ recommends continuing
with wound care alone. We define compliance after the treatment decision as in
the trial of wound care period, that is, compliance relates only to wound
management and determined from the administrative wound clinic data.

###  5.2 Learning an optimal dynamic treatment regime

Using the notation and under the assumptions in Section 2, compliance strata
$C_{1}(0)$ was modeled as

$\displaystyle C_{1}|H_{0},A=0$ $\displaystyle\sim\textit{Truncated
Normal}(\mu_{C_{1}},\sigma_{C_{1}},0,1;\boldsymbol{\beta}_{C_{1}})$
$\displaystyle\boldsymbol{\beta}_{C_{1}}\sim MVN(\boldsymbol{0},3I),$
$\displaystyle\quad\sigma_{C_{1}}\sim HalfCauchy(5)$

where $\mu_{C_{1}}$ is a linear predictor with parameters
$\boldsymbol{\beta}_{C_{1}}$ and covariates are baseline patient and wound
characteristics and compliance level during the first month of wound care, $I$
denotes the identity matrix, $\sigma_{C_{1}}$ is the standard deviation, and
$0$ and $1$ are lower and upper boundaries, respectively, of the support of
the truncated distribution. The binary outcome of being healed, alive, and
with the limb intact at 6 months after initiating wound care, $Y$, was modeled
as

$\displaystyle P(Y=1)\sim
Bernoulli(g^{-1}(\eta(H_{1},\boldsymbol{\beta}_{Y}))),\quad\boldsymbol{\beta}_{Y}\sim
MVN(\boldsymbol{0},\Sigma_{Y})$

where $g$ is the canonical link function, $\eta$ is the linear predictor with
parameters $\boldsymbol{\beta}_{Y}$ and covariates $H_{1}$ that include
baseline patient and wound features, baseline compliance, treatment, and
compliance strata. $\Sigma_{Y}$ is a diagonal matrix and a hyperparameter for
the prior for $\boldsymbol{\beta}_{Y}$. Both models were fit using Hamilton
Monte Carlo and the No U-Turn Sampler in Python version 3.9.7 with modules
PyMC3 version 3.11.4, ArviZ version 0.11.2, and Bambi version 0.7.1. Trace
plots, sampler statistics, and model summary statistics were analyzed to
assess chain convergence and model fit.

To learn an optimal dynamic treatment regime, we used the expectation
improvement algorithm with 25 initial random starts and a budget of 25
evaluations to minimize the negative of the value function over the policy
parameter space $[0,1]\times(0,100]$. To evaluate the value function, we first
sampled counterfactuals by constructing a hypothetical population of 10,000
samples with replacement from the observed data and drew $C_{1}|H_{1},C_{0}$
and $Y|H_{2}$ from their respective posterior predictive distributions where
treatments were fixed as according to the rule implied by $\theta_{1}$ and
$\theta_{2}$. Repeating this 30 times and averaging over the counterfactual
outcomes yielded a draw from distribution of the value function under the
given rule. The Bayesian optimization step was conducted in Python version
3.9.7 with the module scikit-optimize version 0.9.0.

###  5.3 Results and characterization of dynamic treatment regimes

Our analysis yielded an optimal DTR:
$\widehat{d}^{\text{opt}}(C_{0},W_{0};\widehat{\theta}_{1},\widehat{\theta}_{2})=\mathbbm{1}(C_{0}<0.23\text{
or }W_{0}>89.8)$. This means that for PAD patients with wounds who after a
1-month trial of wound management have not healed their wound, surgical
revascularization in addition to continued wound care should be recommended
for those whose completed less than 1/4 of their wound management appointments
during the first month of outpatient wound care and for those whose initial
wound size is larger than 89.8 cm2. In the observed data, 34.2% of patients
were healed, alive, and with their limb intact at 6 months after initializing
wound care. The estimated value of $\widehat{d}^{\text{opt}}$ is 46.9%, a
12.2% improvement.

In addition to learning an optimal dynamic treatment regime, the goal of our
analysis included characterizing the class of regimes under investigation to
aid in translation of the information encoded in and proximal to optimal
dynamic treatment regimes estimation to those making clinical decisions,
clinicians with their patients. That is, we sought to exploit information
learned through the optimal dynamic treatment regime learning process to
supplement and enrich our understanding of the estimated optimal DTR, the
policy class more generally, and the role and relevance of partial compliance
in treatment decision-making. To do this, we first visualized the policies
evaluated during the expectation improvement step to understand how different
policies relate in value and across the parameters using R version 4.0.4 and
plotly version 4.10.0 (Figure 5). Each point represents a policy in the policy
parameter space. The axes denoted as “Baseline compliance (prop.)” and
“Initial wound size ($\text{cm}^{2}$)” correspond to the policy threshold
parameters $\theta_{1}$ and $\theta_{2}$, and the vertical axis denotes the
estimated value for the corresponding policy. The size of the points
corresponds to the precision of the estimate, the inverse of the standard
deviation of draws from the posterior distribution of the value function. From
inspection, it can be seen that there are a number of explored policies with
similarly optimal values (nearly the same height in the z-axis) and that
optimality of the the rule is closely tied to the policy parameter value for
baseline compliance. We scaled the points in the figure to be proportional to
the precision of the value estimate of the policy to understand which policies
we were more sure of their value and which policies we were less certain of
their value.

Figure 5: Two views of the values of policies evaluated during optimization

In addition analyzing the policies explored during the Bayesian optimization
step, we also examined the mean of a Gaussian process surrogate for the value
function over the parameter space. One could use directly the Gaussian process
used for expectation improvement, although we, for ease of exposition, chose
to fit a new Gaussian process with a Matern kernel with noise using the saved
counterfactual draws from the Bayesian G-computation step. Across a $100\times
100$ grid of the parameter space, we evaluated the mean of the fitted Gaussian
process using Python version 3.9.7, scikit-learn version 0.24.2. We plotted
the results (using R version 4.0.4, ggplot2 version 3.3.5), shown in Figure 6,
as level sets represented by different colors where each level represents a 5%
range in the approximated proportion of individuals that would be healed,
alive, and with their limb intact at 6 months if they were to follow the
policy implied by the parameters for baseline compliance and initial wound
size. We chose levels sets, as opposed to a continuous gradient, to identify
policies with nearly equivalent clinically significant meaning. Each point
represents a policy in the policy parameter space that was evaluated during
the expectation improvement step. The estimated optimal policy is denoted by a
white star. Strikingly, a swath of policies, those that recommend
revascularization along with continued wound care when baseline compliance is
less than about 35% or initial wound size is greater than 12 cm2, yield nearly
equivalent approximated values suggesting that following a number of policies
would lead to outcomes similar to that of the optimal policy. This means that
clinicians and their patients may have some “wiggle room” in optimal decision
making; it also shows under which policies near-optimality drops off. This
latter point may be of particular importance for patients who may face
external barriers (e.g. lack of transportation, child care challenges, work
obligations, etc.) to completing their wound care appointments. Finally we
note that after hyperparameter tuning, the fitted Gaussian process has
lengthscale 1.65 for the baseline compliance parameter and 103 for the initial
wound size parameter ($\sigma^{2}=0.244$, $\nu=1.5$, noise = 0.000274)
providing some evidence that the class of policies of interest, the parameter
for baseline compliance plays a bigger role in determining the value of the
policy than the parameter for initial wound size.

Figure 6: Gaussian process surrogate mean for the value over the policy
parameter space

##  6 Discussion

In this paper, we describe a strategy for computationally efficient learning
of optimal DTRs using Bayesian optimization methodology. Perhaps even more
importantly, we attempt to narrow the gap between precision medicine
statistical analyses and translatable clinical evidence by expanding the
question set being asked and answered in precision medicine. Instead of
focusing solely on finding optimal DTRs, we also show and argue for
identifying policies that yield equally optimal results in terms of clinical
significance and exploring how the value of policies change across policy
parameters which may give insight into what features are driving the benefits
gained from precision medicine decision support and which features are less
important for prescribing treatments. By asking new questions and providing
new types of evidence, we get us closer to realizing the promise of precision
medicine, matching the right treatment to the right patient at the right time,
in the real world.

There are limitations of our work. For example, our approach hinges on having
a good estimator for the value across the policy parameter space.
Additionally, it is not clear how large the exploration budget should be for
Bayesian optimization and guarantees of finding a good optimal policy depend
on the underlying value function at least being sufficiently smooth. Finally,
while simple policy classes are easy to understand and can be visualized well,
picking the right form of policies and the size of the policy class remains a
challenge.

A wide range of future work is motivated by the ideas we have proposed. For
example, extending and analyzing the approach in the multi-stage setting or
with more complex parametric decision rules such as trees and lists. Moreover,
our work is a first pass at expanding the targets of a precision medicine
analysis to generate evidence to directly support translation to the clinic.
Future work for the precision medicine field at large may include better
understanding of how clinicians and patients use the type of evidence we have
suggested generating and the questions and answers it inspires.

## Acknowledgements

The first and fourth authors were supported by a grant from the National
Institute on Drug Abuse (NIDA) (No. R01 DA048764). The project described was
supported by the National Center for Advancing Translational Sciences (NCATS),
National Institutes of Health, through Grant Award Number UL1TR002489. The
content is solely the responsibility of the authors and does not necessarily
represent the official views of the NIH.

## References

  * Angrist et al. (1996) Angrist, J. D., Imbens, G. W., and Rubin, D. B. (1996), “Identification of Causal Effects Using Instrumental Variables,” Journal of the American Statistical Association, 91, 444–455.
  * Artman et al. (2020) Artman, W. J., Ertefaie, A., Lynch, K. G., McKay, J. R., and Johnson, B. A. (2020), “Adjusting for Partial Compliance in SMARTs: a Bayesian Semiparametric Approach,” arXiv:2005.10307.
  * Browder et al. (2022) Browder, S. E., Ngeve, S. M., Hamrick, M. E., Wood, J. E., Parodi, F. E., Pascarella, L. E., Farber, M. A., Marston, W. A., and McGinigle, K. L. (2022), “Analysis of wound healing time and wound-free period in patients with chronic limb-threatening ischemia treated with and without revascularization,” J Vasc Surg.
  * Chakraborty and Moodie (2013) Chakraborty, B. and Moodie, E. E. M. (2013), Statistical methods for dynamic treatment regimes: reinforcement learning, causal inference, and personalized medicine, New York, NY: Springer.
  * Cui and Tchetgen Tchetgen (2021) Cui, Y. and Tchetgen Tchetgen, E. (2021), “A Semiparametric Instrumental Variable Approach to Optimal Treatment Regimes Under Endogeneity,” Journal of the American Statistical Association, 116, 162–173, doi: 10.1080/01621459.2020.1783272.
  * Ding and Lu (2017) Ding, P. and Lu, J. (2017), “Principal stratification analysis using principal scores,” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 79, 757–777.
  * Frangakis and Rubin (2002) Frangakis, C. E. and Rubin, D. B. (2002), “Principal stratification in causal inference,” Biometrics, 58, 21–9.
  * Frazier (2018) Frazier, P. I. (2018), “A tutorial on Bayesian optimization,” arXiv preprint arXiv:1807.02811.
  * Ghosal and Van der Vaart (2017) Ghosal, S. and Van der Vaart, A. (2017), Fundamentals of nonparametric Bayesian inference, vol. 44, Cambridge University Press.
  * Imbens and Rubin (2015) Imbens, G. and Rubin, D. B. (2015), Causal inference for statistics, social, and biomedical sciences : an introduction, New York, NY: New York, NY : Cambridge University Press, 2015.
  * Jiang and Ding (2021) Jiang, Z. and Ding, P. (2021), “Identification of causal effects within principal strata using auxiliary variables,” Statistical Science, 36, 493–508.
  * Jo and Stuart (2009) Jo, B. and Stuart, E. A. (2009), “On the use of propensity scores in principal causal effect estimation,” Statistics in Medicine, 28, 2857–2875.
  * Jones et al. (1998) Jones, D. R., Schonlau, M., and Welch, W. J. (1998), “Efficient global optimization of expensive black-box functions,” Journal of Global optimization, 13, 455–492.
  * Lavori and Dawson (2000) Lavori, P. W. and Dawson, R. (2000), “A Design for Testing Clinical Strategies: Biased Adaptive within-Subject Randomization,” Journal of the Royal Statistical Society. Series A (Statistics in Society), 163, 29–38.
  * Lavori and Dawson (2004) — (2004), “Dynamic treatment regimes: practical design considerations,” Clinical Trials, 1, 9–20.
  * Murphy and Bingham (2009) Murphy, S. A. and Bingham, D. (2009), “Screening Experiments for Developing Dynamic Treatment Regimes,” Journal of the American Statistical Association, 104, 391–408.
  * Rasmussen and Williams (2005) Rasmussen, C. E. and Williams, C. K. I. (2005), Gaussian Processes for Machine Learning, The MIT Press.
  * Robins (1986) Robins, J. (1986), “A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect,” Mathematical Modelling, 7, 1393–1512.
  * Robins et al. (2000) Robins, J. M., Hernan, M. A., and Brumback, B. (2000), “Marginal structural models and causal inference in epidemiology,” Epidemiology, 11, 550–60.
  * Rubin (1974) Rubin, D. B. (1974), “Estimating causal effects of treatments in randomized and nonrandomized studies,” Journal of educational psychology, 66, 688–701.
  * Rubin (1980) — (1980), “Randomization Analysis of Experimental Data: The Fisher Randomization Test Comment,” Journal of the American Statistical Association, 75, 591–593.
  * Splawa-Neyman et al. (1990) Splawa-Neyman, J., Dabrowska, D. M., and Speed, T. P. (1990), “On the Application of Probability Theory to Agricultural Experiments. Essay on Principles. Section 9,” Statistical Science, 5, 465–472, 8.
  * Zhang et al. (2012) Zhang, B., Tsiatis, A. A., Laber, E. B., and Davidian, M. (2012), “A robust method for estimating optimal treatment regimes,” Biometrics, 68, 1010–8.

##  Appendix A Identification of $\mathcal{V}(d)$

In this section, we show how $\mathcal{V}(d)$ is identified from the observed
data. Using a similar argument as in (Jiang and Ding, 2021), observe that
$\mathbb{E}[Y(a)|C_{1}(0)]$ for $a=0,1$ is identified from the observed data
given that the SUTVA, NUC, PI, and OSPC assumptions hold and that the family
of distributions $\\{P(C_{1}(0)|X=x)\\}$ is complete. In particular,

$\displaystyle\mathbb{E}[Y|A_{1}=a,X=x]=$ $\displaystyle\mathbb{E}[Y(a)|X=x]$
(SUTVA, NUC) $\displaystyle=$ $\displaystyle\mathbb{E}[Y(a)|C_{1}(0),X=x]$
(PI, OSPC) $\displaystyle=$
$\displaystyle\int_{0}^{1}\mathbb{E}[Y(a)|C_{1}(0)=c_{1}]dP(C_{1}(0)|X=x)$

is an integral equation that uniquely identifies
$\mathbb{E}[Y(a)|C_{1}(0)=c_{1}]$ from the observed data. Identification of
$\mathcal{V}(d)$ follows by noting that

$\displaystyle\mathcal{V}(d(x))=$
$\displaystyle\mathbb{E}\\{\mathbb{E}[\mathbbm{1}(d(x)=0)Y(0)+\mathbbm{1}(d(x)=1)Y(1)|X=x]\\}.$

##  Appendix B Additional simulation results

Table 1: Mean squared error between value of optimal DTR and estimated optimal DTR. |  | IPW | sIPW | G-comp | AIPWE
---|---|---|---|---|---
N |  CV or $w$ | MSE | MC error | MSE | MC error | MSE | MC error | MSE | MC error
Simulation Setting 1
200 | 0.75 | 0.0043 | 0.0002 | 0.0005 | 0.0000 | 0.0002 | 0.0000 | 0.0002 | 0.0000
1.0 | 0.0057 | 0.0002 | 0.0008 | 0.0000 | 0.0004 | 0.0000 | 0.0004 | 0.0000
1.25 | 0.0073 | 0.0003 | 0.0014 | 0.0001 | 0.0007 | 0.0000 | 0.0007 | 0.0000
500 | 0.75 | 0.0017 | 0.0001 | 0.0002 | 0.0000 | 0.0001 | 0.0000 | 0.0001 | 0.0000
1.00 | 0.0022 | 0.0001 | 0.0004 | 0.0000 | 0.0002 | 0.0000 | 0.0002 | 0.0000
1.25 | 0.0032 | 0.0001 | 0.0006 | 0.0000 | 0.0003 | 0.0000 | 0.0003 | 0.0000
1000 | 0.75 | 0.0009 | 0.0000 | 0.0001 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000
1.00 | 0.0011 | 0.0000 | 0.0002 | 0.0000 | 0.0001 | 0.0000 | 0.0001 | 0.0000
1.25 | 0.0015 | 0.0001 | 0.0003 | 0.0000 | 0.0001 | 0.0000 | 0.0001 | 0.0000
Simulation Setting 2
200 | 0.75 | 0.0008 | 0.0000 | 0.0007 | 0.0000 | 0.0048 | 0.0002 | 0.0050 | 0.0002
1.00 | 0.0029 | 0.0001 | 0.0028 | 0.0001 | 0.0722 | 0.0009 | 0.0029 | 0.0001
1.25 | 0.0052 | 0.0002 | 0.0038 | 0.0002 | 0.0833 | 0.0009 | 0.0061 | 0.0002
500 | 0.75 | 0.0004 | 0.0000 | 0.0004 | 0.0000 | 0.0024 | 0.0001 | 0.0025 | 0.0001
1.00 | 0.0014 | 0.0001 | 0.0014 | 0.0001 | 0.0793 | 0.0006 | 0.0013 | 0.0001
1.25 | 0.0029 | 0.0001 | 0.0016 | 0.0001 | 0.0835 | 0.0006 | 0.0042 | 0.0001
1000 | 0.75 | 0.0018 | 0.0002 | 0.0004 | 0.0000 | 0.0016 | 0.0001 | 0.0016 | 0.0001
1.00 | 0.0010 | 0.0001 | 0.0009 | 0.0001 | 0.0846 | 0.0005 | 0.0009 | 0.0001
1.25 | 0.0014 | 0.0001 | 0.0010 | 0.0000 | 0.0855 | 0.0004 | 0.0040 | 0.0001
Simulation Setting 3
200 | 0.75 | 0.0045 | 0.0002 | 0.0047 | 0.0002 | 0.0254 | 0.0006 | 0.0047 | 0.0002
1.00 | 0.0031 | 0.0002 | 0.0032 | 0.0002 | 0.0151 | 0.0004 | 0.0032 | 0.0002
1.25 | 0.0040 | 0.0002 | 0.0042 | 0.0002 | 0.0118 | 0.0003 | 0.0045 | 0.0002
500 | 0.75 | 0.0039 | 0.0001 | 0.0040 | 0.0001 | 0.0254 | 0.0005 | 0.0027 | 0.0001
1.00 | 0.0013 | 0.0001 | 0.0013 | 0.0001 | 0.0176 | 0.0003 | 0.0013 | 0.0001
1.25 | 0.0015 | 0.0001 | 0.0015 | 0.0001 | 0.0105 | 0.0003 | 0.0019 | 0.0001
1000 | 0.75 | 0.0038 | 0.0001 | 0.0039 | 0.0001 | 0.0246 | 0.0004 | 0.0020 | 0.0001
1.00 | 0.0008 | 0.0000 | 0.0007 | 0.0000 | 0.0193 | 0.0002 | 0.0007 | 0.0000
1.25 | 0.0009 | 0.0000 | 0.0008 | 0.0000 | 0.0112 | 0.0002 | 0.0010 | 0.0001
Table 2: Difference between the true value function and the surrogate value function. |  | IPW | sIPW | G-computation | AIPWE
---|---|---|---|---|---
N | $\sigma$ | $\|\cdot\|_{1}$ | $\|\cdot\|_{2}$ | $\|\cdot\|_{1}$ | $\|\cdot\|_{2}$ | $\|\cdot\|_{1}$ | $\|\cdot\|_{2}$ | $\|\cdot\|_{1}$ | $\|\cdot\|_{2}$
Simulation Setting 1
200 | 0.75 | 0.0499 | 0.0569 | 0.0272 | 0.0351 | 0.0210 | 0.0278 | 0.0210 | 0.0277
1.00 | 0.0594 | 0.0660 | 0.0326 | 0.0404 | 0.0243 | 0.0303 | 0.0243 | 0.0303
1.25 | 0.0693 | 0.0752 | 0.0391 | 0.0466 | 0.0285 | 0.0338 | 0.0285 | 0.0338
500 | 0.75 | 0.0332 | 0.0390 | 0.0196 | 0.0262 | 0.0161 | 0.0224 | 0.0161 | 0.0224
1.00 | 0.0376 | 0.0427 | 0.0224 | 0.0283 | 0.0173 | 0.0222 | 0.0173 | 0.0222
1.25 | 0.0459 | 0.0501 | 0.0262 | 0.0315 | 0.0196 | 0.0236 | 0.0196 | 0.0236
1000 | 0.75 | 0.0252 | 0.0306 | 0.0160 | 0.0223 | 0.0138 | 0.0201 | 0.0138 | 0.0201
1.00 | 0.0280 | 0.0323 | 0.0171 | 0.0221 | 0.0138 | 0.0185 | 0.0138 | 0.0185
1.25 | 0.0323 | 0.0357 | 0.0191 | 0.0235 | 0.0150 | 0.0187 | 0.0150 | 0.0187
Simulation Setting 2
200 | 0.75 | 0.0431 | 0.0492 | 0.0402 | 0.0461 | 0.3759 | 0.3983 | 0.3631 | 0.3878
1.00 | 0.0567 | 0.0624 | 0.0568 | 0.0628 | 0.1051 | 0.1311 | 0.0654 | 0.0724
1.25 | 0.0698 | 0.0742 | 0.0656 | 0.0716 | 0.1713 | 0.1856 | 0.1296 | 0.1385
500 | 0.75 | 0.0282 | 0.0339 | 0.0273 | 0.0329 | 0.3758 | 0.3977 | 0.3626 | 0.3867
1.00 | 0.0382 | 0.0432 | 0.0384 | 0.0436 | 0.0889 | 0.1185 | 0.0445 | 0.0502
1.25 | 0.0455 | 0.0521 | 0.0416 | 0.0463 | 0.1647 | 0.1777 | 0.1182 | 0.1264
1000 | 0.75 | 0.0311 | 0.0426 | 0.0208 | 0.0266 | 0.3766 | 0.3983 | 0.3612 | 0.3852
1.00 | 0.0284 | 0.0333 | 0.0284 | 0.0336 | 0.0813 | 0.1144 | 0.0319 | 0.0371
1.25 | 0.0317 | 0.0354 | 0.0302 | 0.0344 | 0.1653 | 0.1772 | 0.1188 | 0.1270
Simulation Setting 3
200 | 0.75 | 0.0578 | 0.0659 | 0.0576 | 0.0660 | 0.0903 | 0.1041 | 0.0785 | 0.0879
1.00 | 0.0590 | 0.0661 | 0.0575 | 0.0646 | 0.0779 | 0.0904 | 0.0677 | 0.0752
1.25 | 0.0619 | 0.0679 | 0.0639 | 0.0698 | 0.0809 | 0.0909 | 0.0737 | 0.0798
500 | 0.75 | 0.0401 | 0.0483 | 0.0404 | 0.0486 | 0.0766 | 0.0902 | 0.0611 | 0.0694
1.00 | 0.0406 | 0.0471 | 0.0403 | 0.0470 | 0.0597 | 0.0726 | 0.0452 | 0.0518
1.25 | 0.0414 | 0.0467 | 0.0417 | 0.0470 | 0.0608 | 0.0708 | 0.0515 | 0.0568
1000 | 0.75 | 0.0321 | 0.0406 | 0.0320 | 0.0405 | 0.0714 | 0.0848 | 0.0548 | 0.0625
1.00 | 0.0308 | 0.0376 | 0.0302 | 0.0371 | 0.0500 | 0.0645 | 0.0339 | 0.0404
1.25 | 0.0310 | 0.0363 | 0.0310 | 0.0364 | 0.0500 | 0.0601 | 0.0398 | 0.0449