#  Bayesian Semiparametric Model for Sequential Treatment Decisions with
Informative Timing

Arman Oganisian
Department of Biostatistics, Brown University
Kelly D. Getz
Department of Biostatistics, Epidemiology, and Informatics,
University of Pennsylvania
Todd A. Alonzo
Department of Preventive Medicine, University of Southern California
Richard Aplenc
Division of Oncology, Children’s Hospital of Philadelphia
Jason A. Roy
Department of Biostatistics and Epidemiology, Rutgers University

###### Abstract

We develop a Bayesian semi-parametric model for the estimating the impact of
dynamic treatment rules on survival among patients diagnosed with pediatric
acute myeloid leukemia (AML). The data consist of a subset of patients
enrolled in the phase III AAML1031 clinical trial in which patients move
through a sequence of four treatment courses. At each course, they undergo
treatment that may or may not include anthracyclines (ACT). While ACT is known
to be effective at treating AML, it is also cardiotoxic and can lead to early
death for some patients. Our task is to estimate the potential survival
probability under hypothetical dynamic ACT treatment strategies, but there are
several impediments. First, since ACT was not randomized in the trial, its
effect on survival is confounded over time. Second, subjects initiate the next
course depending on when they recover from the previous course, making timing
potentially informative of subsequent treatment and survival. Third, patients
may die or drop out before ever completing the full treatment sequence. We
develop a generative Bayesian semi-parametric model based on Gamma Process
priors to address these complexities. At each treatment course, the model
captures subjects’ transition to subsequent treatment or death in continuous
time under a given rule. A g-computation procedure is used to compute a
posterior over potential survival probability that is adjusted for time-
varying confounding. Using this approach, we conduct posterior inference for
the efficacy of hypothetical treatment rules that dynamically modify ACT based
on evolving cardiac function.

##  1 Introduction

Our application involves a subset of data from the AAML1031 phase III clinical
trial run by the Children’s Oncology Group (COG). Upon trial enrollment,
subjects with pediatric acute myeloid leukemia (AML) were randomized to a
four-course treatment sequence. At each course, patients were given a standard
chemotherapy treatment comprised of multiple agents which may or may not
include anthracycline (ACT). In addition to the standard chemotherapy, some
patients were also randomized to receive bortezomib. Despite its demonstrated
effectiveness in treating AML, ACT agents are known to be cardiotoxic. The
impact of ACT on overall survival is complex - for some it could increase
survival while for others it may cause cardiotoxicity and early death. This
presents a potential risk-benefit tradeoff associated with ACT. Ahead of each
course, physicians conducted an echocardiogram to measure patients’ left
ventricular ejection fraction (EF) - a marker of cardiac function with lower
values indicating worse function. Along with other information, EF is used to
decide whether to administer ACT or withhold ACT. The precise ways in which EF
is used varies, but often involves withholding ACT if decline in EF is below
some relative threshold (e.g. a 10% decline from baseline (Neuendorff et al.,
2020)) or if EF is below some absolute threshold (e.g. any EF value below 50%)
- but the utility of these thresholds and their effects on survival are still
debated. In a previous trial, COG AAML0531, about 12% of patients experienced
cardiotixicity within a year of treatment initiation. In COG AAML1031, the
percentage is 22% which is likely due to more complete monitoring of patients.
Additionally, patients who experienced cardiotoxicity were found to have lower
overall survival in both the AAML0531 and AAML1031 trials (Getz et al., 2019a,
b). Thus, there is strong clinical interest in estimating (and optimizing) the
impact of ACT treatment on survival.

Fortunately, a subset of the AAML1031 patients were treated at centers that
contribute to a separate Pediatric Health Information System (PHIS) database,
allowing us to assess ACT administration via PHIS billing data. While the data
are rich, they pose several impediments to estimation. First, ACT
administration is not randomized, but modified as needed according to the
patient’s clinical status at each treatment course. Thus, its effect on
survival is confounded by time-varying variables such as EF and other
information available just ahead of each treatment decision. Second, while
each subject receives at most four treatment courses, the timing of the
courses vary by subject. This is partially because AML chemotherapy causes
prolonged neutropenia. Subjects initiate the next course after they achieve
hemotalogic recovery and are clinically stable. Since patients with longer
recovery times may have greater toxicities, the waiting times between
treatment courses may be informative of subsequent ACT decision and survival.
Third, patients may die or drop out of the study before completing the full
sequence so that the total number of treatment courses varies across patients
as well. Naive approaches such as restricting the analysis to only subjects
who have survived through the full treatment sequence will generally yield a
biased subset (likely systematically healthier) of subjects. Common approaches
such as discretizing the time period (e.g. into months) do not respect the
continuous-time nature of the decision-making process, often leading to
difficulties in dealing with patients who initiate treatments between discrete
periods. Standard g-computation approaches (Robins, 1986) were designed for
settings where all subjects receive, say, the second treatment at the same
state of their disease process. These are inappropriate in settings where
timing may confound subsequent treatment and survival.

To address these concerns, we frame ACT assignment as a dynamic treatment
regime in which a rule maps a patient’s available history at each course to a
treatment decision. Under certain assumptions about the underlying causal
structure, we identify the potential survival probability under a given
treatment rule in terms of the cause-specific hazards of subsequent treatment
and death as well as the distribution of the time-varying confounders (e.g.
EF) at each course. It is recognized that even when the causal assumptions
hold, inferences about potential survival probability will be biased if
modeling assumptions are violated. This motivates a Bayesian semiparametric
specification of the cause-specific hazards that is more robust to
misspecification. In particular, we specify proportional hazard models in
which the baseline hazards follow a Gamma Process prior. Together, these
semiparametric cause-specific hazards characterize patients’ transition
between successive treatments and death in continuous time. Markov Chain Monte
Carlo (MCMC) samplers for such models can be difficult to tune due to the
large number of parameters present even with just four time periods, so we
outline an adaptive sampling scheme for posterior inference. Finally, a
g-computation procedure is outlined that simulates from the continuous-time
transition process under hypothetical treatment rules and computes the
corresponding potential survival probability. The end product is full
posterior inference for potential survival rates (and their functionals) under
treatment rules of interest which are appropriately adjusted for time-varying
confounding and informative waiting times while offering robustness to
misspecification.

While causal estimation and optimization with dynamic treatment rules have
been discussed for some time (Robins, 1986; Murphy, 2003), Bayesian
nonparametric (BNP) estimation strategies in particular have only recently
grown in popularity. A key draw of BNP is the ability to incorporate posterior
uncertainty about the decision process into a single posterior distribution
for the target estimand while leveraging flexible models that impose only weak
structure on the data. Moreover, simulations from this posterior can be
transformed to allow inference about functionals of interest. We build upon
this literature with several distinct contributions. For instance, Xu et al.
(2016) proposed a BNP model for sequential AML treatments based on Dependent
Dirichlet Processes, however their setting only required adjustment for a
single set of baseline confounders. In our setting, we have important time-
varying confounders whose evolution must be modeled over time. Guan et al.
(2020) and Hua et al. (2021) also develop Bayesian approaches in a setting
with random, potentially informative treatment times. Specifically, Guan et
al. (2020) consider optimizing a value function based on a proportion outcome,
which was not subject to right-censoring. Hua et al. (2021) develop a Bayesian
point-process model of a decision process with a survival outcome where all
decision points of interest are before the observed time on study and a
Weibull proportional hazard submodel is fit for the event time. In contrast,
we have a fixed number of decision points of interest, and not all subjects
will survive through the full sequence. Moreover, we take a semiparametric
approach to modeling the event time via Gamma Process priors. Murray et al.
(2017) propose a design of a randomized clinical trial assessing the efficacy
of a treatment sequence where the timing of a third-line therapy is random. In
our study, ACT administration is not randomized, which necessitates careful
thought about potential outcomes and confounding adjustment.

In the following sections we describe the observed AAML1031 data structure,
formalize and identify our target estimand in terms of potential outcomes, and
discuss Bayesian inference. As our approach is tailored to deal with the
complexities of our data, we will frequently motivate particular features by
referencing the data, but a complete data analysis will be reserved for the
last section.

##  2 Description of AAML1031-PHIS Data

Subject 1$Y_{1}=0$$Y_{2}$$Y_{3}$$Y_{4}$$T$$W_{1}$$W_{2}$$W_{3}$$W_{4}$Subject
2$Y_{1}=0$$Y_{2}$$Y_{3}$$C$$W_{1}$$W_{2}$$W_{3}$Subject
3$Y_{1}=0$$Y_{2}$$T$$W_{1}$$W_{2}$

Figure 1: Some possible patient trajectories in our data. While at most $K=4$
treatments are possible, the timing of the $k^{th}$ treatment, $Y_{k}$,
varies. Patients may die after completing the full sequence (Subject 1) but
may also drop out (Subject 2) or die (Subject 3) before completing the full
sequence.

Consider a scenario where patients move through a sequence of $K$ treatment
courses, indexed by $k=1,2,\dots,K$. In our AML study, $K=4$. At the beginning
of course $k$, physicians monitor a vector of $P$ patient characteristics,
$L_{k}\in\mathcal{L}_{k}$, to guide their decision to include ACT or withhold
ACT in favor of treating with non-anthracyclines (nACT) . We let
$A_{k}\in\mathcal{A}_{k}=\\{0,1\\}$ denote this decision with $A_{k}=1$
indicating ACT inclusion. In our analysis, these include time-constant
features such as sex and genotype information as well as time-varying features
like ejection fraction and presence bloodstream infection. These features
typically affect both treatment decisions and survival, making them important
time-varying confounders which must be adjusted for in order to estimate the
impact of the treatment sequence on survival. Importantly, while the maximum
number of treatment courses is fixed at $K$, patients initiate each course at
different times depending on how quickly they recover from their previous
course. Treatment course $k$ is initiated at time $Y_{k}>0$ with available
information $H_{k}=(\bar{L}_{k},\bar{A}_{k-1},\bar{Y}_{k})\in\mathcal{H}_{k}$,
where the overbar $\bar{X}_{k}=(X_{1},X_{2},\dots,X_{k})$ denotes the history
of $X$ up to and including $X_{k}$. We use underbar notation to denote future
observations, e.g. $\underaccent{\bar}{X}_{k}=(X_{k},X_{k+1},\dots,X_{K})$.
All patients undergo their first treatment course, $A_{1}$, at time
$Y_{1}:=0$. The remaining treatment times $Y_{k}$ for $k>1$ vary across
subjects and are measured relative to $Y_{1}=0$. We adopt the convention that
at time $Y_{k}$, $L_{k}$ is measured before $A_{k}$ is administered. Finally,
let $T>0$ denote the patient’s survival time, also measured from $Y_{1}=0$. In
our data, subjects may die before completing the full sequence of $K$
treatments. Additionally, patients may drop out at any point after the first
treatment, inducing right-censoring. This means that each subject will undergo
a random number of treatment courses, $\kappa\in\\{1,2,\dots,K\\}$. Since all
subjects undergo the first treatment course, $\kappa\neq 0$. Let $W=\min(T,C)$
be the observed time in study, which is the minimum of death time and
censoring time, $C>0$, and let $\delta=I(T<C)$ be the corresponding indicator.
For course $k=1,\dots,K$ we define a sequence of waiting times until the next
event. Specifically, for subjects still in the study at the start of treatment
course $k$ (i.e. $\kappa\geq k$), $W_{k}$ is the time from treatment $k$
initiation, $Y_{k}$, to either death at time $T$, drop out at time $C$, or
subsequent treatment at time $Y_{k+1}$ \- whichever comes first. Thus, the
observed waiting time is $W_{k}=\min(Y_{k+1},T,C)-Y_{k}$ and the corresponding
indicator $\delta_{k}\in\\{1,0,-1\\}$ indicates a transition to death,
censoring, or subsequent treatment respectively. Figure 1 depicts a few
possible patient trajectories. We additionally define cause-specific
transition indicators for death and subsequent treatment,
$\delta_{Tk}=I(\delta_{k}=1)$ and $\delta_{Yk}=I(\delta_{k}=-1)$. At the last
treatment course, there is no subsequent treatment so $W_{K}=\min(T,C)-Y_{K}$,
$\delta_{K}\in\\{1,0\\}$, $\delta_{TK}=\delta_{K}$ and $\delta_{YK}$ is
undefined.

Figure 2:  Left: Distribution of waiting times between treatment courses among
those at-risk. Note the long tail as some patients take especially long to
recover before starting their next course. Right: scatterplot of observed
waiting times between course 2 and 3 ($W_{2}$) and course 4 and 3 ($W_{3}$)
among those who had all four treatment courses visualizing positive
association in waiting times.

In this setup, the observed time for a subject who underwent $\kappa$
treatment courses can be expressed as the sum of the waiting times between
events $W=\sum_{j=1}^{\kappa}W_{j}$. We equivalently write the history
available just before treatment decision $A_{k}$ as
$H_{k}=(\bar{L}_{k},\bar{A}_{k-1},\bar{Y}_{k})=(\bar{L}_{k},\bar{A}_{k-1},\bar{W}_{k-1})$,
since the $\bar{W}_{k}$ are just transformations of $\bar{Y}_{k}$. We observe
data on $i=1,2,\dots,n$ independent subjects. The data for subject $i$
consists of
$\mathcal{D}_{i}=(\bar{L}_{\kappa_{i}},\bar{A}_{\kappa_{i}},\bar{W}_{\kappa_{i}},\bar{\delta}_{Y\kappa_{i}},\bar{\delta}_{T\kappa_{i}})$
and we let $\mathcal{D}=\\{\mathcal{D}_{i}\\}_{i=1}^{n}$ denote the full data.

The left panel of Figure 2 shows appreciable variation in waiting times around
35 days, with a tail of subjects who start their next course later than
average. At any given course, the time since previous course potentially
informs the treatment decision (e.g., subjects who took longer than usual to
recover should perhaps be treated conservatively) and may also be related to
survival time accrued afterward. Indeed, the right panel of Figure 2 shows
that observed waiting times between courses are positively correlated in these
data. These complications, along with time-varying confounding of the
treatment sequence, motivate the formal causal approach to estimating
treatment effects described in the next section.

##  3 Target Estimand and Identification

Here, we define relevant target estimand in terms of potential outcomes had we
assigned treatment via some rule and state the main identification result that
expresses this estimand in terms of objects that can be estimated from the
observed data. Suppose the treatment decision at course $k$ is the output of
decision rule $r_{k}:\mathcal{H}_{k}\rightarrow\mathcal{S}_{k}(h_{k})$ that
maps from the available history to a treatment in the feasible set,
$d_{k}^{r}\in\mathcal{S}_{k}(h_{k})\subset\mathcal{A}_{k}$, and denote the
collection of rules as $r=(r_{1},r_{2},\dots,r_{K})$. It is necessary to
define feasible subsets of $\mathcal{A}_{k}$ since, depending on a subject’s
history $h_{k}$, some options in $\mathcal{A}_{k}$ may not be possible. For
instance, in the AAML1031, all subjects were generally ineligible for ACT at
the third course and so $\mathcal{S}_{3}(h_{3})=\\{0\\}\subset\mathcal{A}_{3}$
for all $h_{3}$. The ACT/nACT treatment sequence is not set at baseline, but
adapts dynamically to each patient’s unique disease progression: at time
$Y_{1}$, we first measure $L_{1}$ then use it to arrive at a decision
$d_{1}^{r}=r_{1}(L_{1},Y_{1})$ and assign $A_{1}=d_{1}^{r}$. At the second
treatment course, at time $Y_{2}$, we measure $L_{2}$ the make decision
$d_{2}^{r}=r_{2}(\bar{L}_{2},A_{1}=d_{1}^{r},\bar{Y}_{2},)=r_{2}(H_{2})$ and
set $A_{2}=d_{2}^{r}$. Similarly, at the $k^{th}$ treatment course, we arrive
at decision $d_{k}^{r}=r_{k}(H_{k})$ and set $A_{k}=d_{k}^{r}$. We let
$\bar{d}_{k}^{r}$ denote the decision history up to and including decision
$k$. For $k=2,\dots,K$, let $\kappa_{k}(\bar{d}_{k-1}^{r})\in\\{0,1\\}$
indicate whether subject would have survived long enough to initiate course
$k$, had previous treatment decision been made via rule $r$. By construction,
$\kappa_{k}(\bar{d}_{k-1}^{r})=1$ implies $\kappa_{j}(\bar{d}_{j-1}^{r})=1$
for all $j<k$ and $\kappa_{k}(\bar{d}_{k-1}^{r})=0$ implies
$\kappa_{j}(\bar{d}_{j-1}^{r})=0$ for all $j>k$. Thus,
$\kappa(r)=\min_{k}\\{k:\kappa_{k}(\bar{d}_{k-1}^{r})=1\\}$ denotes the
potential number of treatments a patient would have received. For
$k=2,\dots,K$ we define potential time of treatment $k$ as
$Y_{k}(r)=Y_{k}(\bar{d}_{k-1}^{r})$ had they been treated according to rule
$r$. Since $\kappa_{1}=1$ and $Y_{1}=0$ for all subjects, they are not
potential outcomes of $r$. We also define potential survival time under rule
$r$, $T(r)=T(\bar{d}^{r}_{\kappa(r)})$, and
$W_{k}(r)=W_{k}(\bar{d}_{k}^{r})=\min(T(r),Y_{k+1}(r))-Y_{k}(r)$ the
corresponding potential waiting time from treatment $k$ to the first of
$Y_{k+1}(r)$ or $T(r)$. Finally, we define the potential confounder value
measured at treatment course $k$, $L_{k}(r)=L_{k}(\bar{d}_{k-1}^{r})$.

Our target estimand is the potential probability of surviving past $t$ under
rule $r$, $\Psi^{r}(t)=P(T(r)>t)$. Note that $\Psi^{r}(t)$ is a function of
potential outcomes, not observed data. As is typical in causal inference
settings, identification assumptions are required to express $\Psi^{r}(t)$ \-
in terms of observed data, $\mathcal{D}$. Throughout, we assume the following
assumptions about the causal mechanism hold:

  1. 1.

Sequential conditional exchangeability: at the $k^{th}$ treatment course,

$\underaccent{\bar}{W}_{k}(\bar{d}_{k}^{r}),\underaccent{\bar}{\kappa}_{k+1}(\bar{d}_{k}^{r}),\underaccent{\bar}{L}_{k+1}(\bar{d}_{k}^{r})\perp
A_{k}\mid\bar{A}_{k-1}=\bar{d}_{k-1}^{r},\bar{L}_{k},\bar{W}_{k-1},\kappa\geq
k$

Among those who underwent course $k$, treatment $A_{k}$ is as good as random
conditional on available history. I.e. assignment is unrelated to future
potential waiting times, confounder trajectories, and number of subsequent
treatments. This can be viewed as emulating a sequentially randomized trial.

  2. 2.

Treatment positivity: $P(A_{k}=a_{k}\mid h_{k},\kappa\geq k)>0$ for all
histories $f(h_{k})>0$. This must hold for each feasible treatment option
$a_{k}\in\mathcal{S}_{k}(h_{k})$. This could be violated structurally if, for
instance, subset of patients with history
$(\bar{a}_{k-1},\bar{l}_{k},\bar{w}_{k})$ would always receive one of the two
treatments - rendering the causal effect in this subset undefined.

  3. 3.

Non-informative censoring: The cause-specific hazard of being censored $w$
time units after $Y_{k}$ is given by $\lambda_{C}(w\mid h_{k},\kappa\geq k)$.
Note that it does not depend on potential outcomes conditional on available
history, $h_{k}$. This ensures that censoring is conditionally random and
rules out situations in which, say, subjects with systematically higher
potential survival times are more likely to be censored. Additionally, we
require $P(C\geq Y_{k}+w\mid W_{k}(\bar{d}_{k}^{r})=w,h_{k},\kappa\geq k)>0$.
At each course there is a positive probability of reaching the next event at
$W_{k}(\bar{d}_{k}^{r})=w$. Else there would be some subgroup in terms of
$h_{1}$ for whom we would never observe the full outcome history.

  4. 4.

SUTVA : For uncensored subjects with observed treatment vector
$\bar{A}_{k}=\bar{a}_{k}$, the potential outcomes under $\bar{a}_{k}$ are
equal to the observed outcomes: $W_{k}(\bar{a}_{k})=W_{k}$,
$L_{k}(\bar{a}_{k})=L_{k}$, and $\kappa_{k}(\bar{a}_{k})=\kappa_{k}$. For

As in standard survival settings, the survival curve, $\Psi^{r}(t)$, is
nonparametrically identifiable only up to the maximum observed death time.
Under these assumptions, the joint distribution of the potential outcomes can
be identified in terms of observed data in settings where both the number of
treatments and timing is random (Tsiatis et al., 2019). Various quantities of
interest can then computed as integrals over this joint. In our case,
$\Psi^{r}(t)=\int_{t}^{\infty}f_{T(r)}(w)dw$. By the total law of probability,
this can be expressed as the joint averaged over the potential number of
treatments
$\Psi^{r}(t)=\int_{t}^{\infty}\sum_{k=1}^{K}f_{T(r),\kappa(r)}(w,\kappa(r)=k)dw$.
Identification of $\Psi^{r}(t)$ is achieved by identifying
$f_{T(r),\kappa(r)}(w,\kappa(r)=k)$ for each $k$. We present the result for
$k=K$ and provide a derivation in the supplement

$\small\begin{split}&f_{T(r),\kappa(r)}(w,\kappa(r)=K)=\int_{\bar{\mathcal{L}}_{K}}\int_{\bar{\mathcal{W}}_{K-1}}\prod_{k=1}^{K-1}\Big{\\{}f_{k}(l_{k}\mid\bar{d}_{k-1}^{r},\bar{l}_{k-1},\bar{w}_{k-1},\kappa\geq
k)\\\ &\ \ \ \times
S_{k}(w_{k}\mid\bar{d}_{k}^{r},\bar{l}_{k},\bar{w}_{k-1},\kappa\geq
k)\lambda_{Yk}(w_{k}\mid\bar{d}_{k}^{r},\bar{l}_{k},\bar{w}_{k-1},\kappa\geq
k)\Big{\\}}\\\ &\ \ \ \times
f_{K}(l_{k}\mid\bar{d}_{K}^{r},\bar{l}_{K-1},\bar{w}_{K-1},\kappa\geq
K)S_{TK}(w-\sum_{k=1}^{K-1}w_{k}\mid\bar{d}_{K}^{r},\bar{l}_{K},\bar{w}_{K-1},\kappa\geq
K)\\\ &\ \ \
\times\lambda_{TK}(w-\sum_{k-1}^{K-1}w_{k}\mid\bar{d}_{K}^{r},\bar{l}_{K},\bar{w}_{K},\kappa\geq
K)\ d\bar{l}_{K}\ d\bar{w}_{K-1}\end{split}$ (1)

This is a version of the g-formula (Robins, 1986) where we integrate over the
joint distribution of the observed confounders and waiting times under
interventions that are dynamically set via rule $r$, $\bar{d}_{K}^{r}$. The
process of computing the integrals involved in (1) is known as g-computation,
which we do via Monte Carlo simulation in practice. Above, $f_{k}$ is the
density function of confounders measured at time point $k$ conditional on
available history. Since there is no previous history at $k=1$,
$\bar{X}_{0}=\emptyset$. So, for instance, in the above we have
$f_{1}(l_{1}\mid\bar{a}_{0},\bar{l}_{0},\bar{w}_{0})=f_{1}(l_{1})$. For
$k=1,2,\dots,K-1$, $\lambda_{Yk}(w\mid-)$ is the cause-specific hazard of
subsequent treatment after treatment $k$

$\begin{split}\lambda_{Yk}(w\mid-)=\lim_{dw\rightarrow 0}dw^{-1}P(w\leq
W_{k}<w+dw,\delta_{k}=-1\mid W_{k}\geq w,-)\end{split}$ (2)

For compactness, the conditioning set is abbreviated as “$-$”. Conditional on
not having died or initiated the subsequent treatment by $w$ time-units after
treatment $k$, $\lambda_{Yk}$ is the instantaneous probability of being
treated again at that time. Similarly, the function $\lambda_{Tk}(w\mid-)$ for
$k=1,2,\dots,K$ is the cause-specific hazard of death after treatment course
$k$ and has the same form as (2) except with $\delta_{k}=1$. This represents
the instantaneous probability of dying $w$ time units after treatment $k$,
conditional on not having died or received treatment $k+1$ by that time. The
function $S_{k}(w\mid-\ )=\exp(-\int_{0}^{w}\\{\lambda_{Tk}(s\mid-\
)+\lambda_{Yk}(s\mid-\ )ds\\})$ is the probability of not having died or
received subsequent treatment within $w$ time-units after treatment $k$ \-
i.e. the probability of being at-risk for subsequent treatment or death at
$w$. Since there are no more treatments after course $K$, $S_{TK}(w\mid-\
)=\exp(-\int_{0}^{w}\lambda_{TK}(s\mid-)ds)$ \- which does not involve the
hazard of subsequent treatment. These functional forms follow from the
standard relationships between hazard, density, and survival functions in the
competing risks literature (Kalbfleisch and Prentice, 2011) \- where after
each treatment, subsequent treatment and death are competing to be the next
event. Alternatively, (1) can be understood in terms of transition
probabilities between states of treatment and death. In (1), the product
$P(\delta_{k}=-1\rightarrow\delta_{k+1}=-1)=S_{k}(w_{k}\mid-)\lambda_{Yk}(w_{k}\mid-)$
is the conditional probability of transitioning to treatment course $k+1$
$w_{k}$ time units after course $k$. Similarly,
$P(\delta_{k}=-1\rightarrow\delta_{k+1}=1)=S_{k}(w_{k}\mid-)\lambda_{Tk}(w_{k}\mid-)$
is the conditional probability of transitioning to death $w_{k}$ time units
after course $k$. Thus these cause-specific hazards characterize how subjects
transition to either subsequent treatment or death states, in continuous time,
under rule $r$. In the following sections, we propose Bayesian semi-parametric
models for the unknowns in (1), outline an MCMC procedure for posterior
sampling, and describe a g-computation procedure for computing the integrals
in (1).

##  4 Generative Bayesian Semiparametric Models

The quality of our inference about $\Psi^{r}(t)$ relies on correctly specified
models for $f_{k}$, $\lambda_{Yk}$, and $\lambda_{Tk}$. Misspecification of
these nuisances may bias estimates of $\Psi^{r}(t)$, even if the causal
identification assumptions hold. To be more robust to misspecification, we
propose proportional hazard models with a flexible baseline hazard
specification. For the hazard of subsequent treatment after treatment $k$, we
set

$\begin{split}\lambda_{Yk}(w\mid
H_{k},\lambda_{0Yk},\beta_{Yk})&=\lambda_{0Yk}(w)\exp\big{(}H_{k}^{\prime}\beta_{Yk}\big{)}\\\
\Lambda_{0Yk}&\sim GP(\alpha_{Yk}\Lambda_{Yk}^{*})\\\ \beta_{Yk}&\sim
f_{\beta_{Yk}}(\beta_{Yk})\end{split}$ (3)

Here, $\Lambda_{0Yk}(w)=\int_{0}^{w}\lambda_{0Yk}(u)du$ is the baseline
cumulative hazard and $GP(\alpha\Lambda_{Yk}^{*})$ denotes the Gamma Process
prior (Kalbfleisch, 1978) over the unknown baseline cumulative hazard.
$f_{\beta_{Yk}}$ is a prior over the multiplicative effects, $\beta_{Yk}$.
Throughout we set this to multivariate normal with density
$f_{\beta_{Yk}}=N(0,cI)$. Note that on a hazard ratio scale, $c=1$ is fairly
informative - putting most prior mass on ratios in $\exp(\pm 3\cdot 1)$. The
Gamma Process (GP) is a stochastic process that generates random positive,
non-decreasing functions. Realizations from the GP in (3) are random
cumulative hazards centered around a specified prior cumulative hazard
$\Lambda_{Yk}^{*}$ with dispersion controlled by $\alpha_{Yk}$. The key
property of GPs is that if $\Lambda_{0Yk}(w)$ follows a GP, then the hazard
rate over any finite interval follows a Gamma distribution that is centered
around the hazard rate over that interval given by the prior
$\lambda_{Yk}^{*}$. This motivates a piecewise constant specification of
$\lambda_{0Yk}(w)$ over a pre-specified partition of time
$u_{Yk1}=0<u_{Yk2}<u_{Yk1}<\dots<u_{J_{Yk}}=\max(W_{k})$:
$\lambda_{0Yk}(w)=\sum_{j=1}^{J_{Yk}-1x}\mathbbm{1}_{[u_{j},u_{j+1})}(w)\lambda_{0Yk}^{j}$.
The Gamma Process induces independent Gamma priors on the $j^{th}$ hazard
rates $(u_{j+1}-u_{j})\lambda_{0Yk}^{j}\sim
Gam\big{(}\alpha_{Yk}(\Lambda_{Yk}^{*}(u_{j+1})-\Lambda_{Yk}^{*}(u_{j})),\alpha_{Yk}\big{)}$.
A finer partition leads to a higher-dimensional, flexible hazard model. On the
other extreme, a single partition completely shrinks to a one-dimensional
constant hazard model. With a fine partition, this provides an alternative to,
say, a $Weibul(a,s)$ that imposes a rigid form for the hazard,
$\lambda_{0Yk}(w)=\frac{a}{s^{a}}w^{a-1}$. We take an empirical approach that
sets $\Lambda_{Yk}^{*}$ to an exponential (linear) cumulative hazard with rate
equal to inverse sample mean of the waiting time. Throughout, we specify a
fine partition with $\alpha_{Yk}$ close to zero. For small $\alpha_{Yk}$,
realizations of the GP are widely dispersed around $\Lambda_{Yk}^{*}$ \-
representing a weak prior. Thus the posterior can deviate from this linear
prior if the data suggests nonlinearity. An analogous model is specified for
the waiting time to death after treatment $k$, $\lambda_{Tk}(w\mid
H_{k},\Lambda_{0Tk},\beta_{Tk})=\lambda_{0Tk}(w)\exp\big{(}H_{k}^{\prime}\beta_{Tk}\big{)}$,
with priors $\Lambda_{0Tk}\sim GP(\alpha_{Tk}\Lambda_{Tk}^{*})$ and
$\beta_{Tk}\sim f_{\beta_{Tk}}(\beta_{Tk})$.

While the specification of confounder distribution is application specific, it
is common to assume conditional independence among the $P$ covariates in
$L_{k}$ at each time point such that
$f_{k}(l_{k}\mid\bar{H}_{k-1})=\prod_{p=1}^{P}f_{kp}(l_{kp}\mid\bar{H}_{k-1})$
and model each component separately. For generality of presentation, we write
$f_{k}(l_{k}\mid\bar{H}_{k},\eta_{k})$, where $\eta_{k}$ is a vector of
parameters that governs the joint distribution of the confounders at time $k$.
These models can be as flexible or smooth as needed depending on sample size
considerations. For instance, continuous $L_{kp}$ can be modeled as Gaussian
with conditional mean dependent on $H_{k}$ through, say, a regression
$E[L_{kp}\mid\bar{H}_{k};\eta_{kp}]=H_{k}^{\prime}\eta_{kp}$, where
$\eta_{kp}$ is a vector of regression coefficients that may include an
intercept if a constant 1 is included in $H_{k}$. If more data is available we
could be more flexible and specify
$E[L_{kp}\mid\bar{H}_{k};\eta_{kp}]=g(h_{k};\eta_{kp})$, where $g$ is some
function of history, and specify a BART (Chipman et al., 2010) prior on the
function $g\sim BART$. Here, $\eta_{kp}$ would consist of the BART tree
parameters. While time-varying covariates must usually be modeled as a
function of history, nonparametric approaches like the Bayesian bootstrap can
be used for the baseline confounder distribution, $f_{1}$. This is essentially
a Bayesian analogue of the empirical distribution. Similarly, appropriate
logistic models or BART models can be specified for binary covariates. Under
the non-informative censoring assumption, all subjects who undergo course $k$
contribute to the likelihood for $\eta_{k}$ and to the hazards of the next
possible events, $\lambda_{Tk}$ and $\lambda_{Yk}$ and time-varying covariate
distributions. Details of the likelihood construction are given in the
Supplement.

In all, this model is semi-parametric in the sense that we have a flexible
specification for the waiting time hazards. It is generative in the sense that
we model the full joint distribution of the data rather than, say, just the
moments of some conditional outcome model. Generative models are useful
because functionals of the joint can be computed easily by simulating from the
model. For instance, in (1), $\Psi^{r}(t)$ is just a functional of the unknown
confounder distributions, $f_{k}$, and the unknown cause specific hazards,
$\lambda_{Yk}$ and $\lambda_{Tk}$, which are completely characterized by the
parameters $\omega=\\{\omega_{k}\\}_{k=1}^{K}$, where
$\omega_{k}=\Big{\\{}\lambda_{0Tk},\beta_{Tk},\lambda_{0Yk},\beta_{Yk},\eta_{k}\Big{\\}}$.
Note that for the last course, $K$, $\omega_{K}$ does not contain
$\\{\lambda_{0YK},\beta_{YK}\\}$ since there is no subsequent course. Thus,
from a Bayesian perspective, a posterior over $\omega$ induces a posterior
over relevant functionals such as $\Psi^{r}(t)$, provided that we compute the
requisite integrals.

##  5 Adaptive Algorithm for Posterior Computation

Inference for the required models follows from the joint posterior over the
unknown cause-specific hazards and the parameters governing the sequence of
confounder distributions,
$f_{\Omega|\mathcal{D}}(\omega_{1},\omega_{2},\dots,\omega_{K}\mid\mathcal{D})$.
In the following we will use bracket notation to denote generic posteriors.
E.g. the joint posterior above will be denoted as
$[\omega_{1},\omega_{2},\dots,\omega_{K}]$, where conditioning on
$\mathcal{D}$ is suppressed for compactness. Since this posterior is not
available in closed form, we implement a blocked Metropolis-in-Gibbs sampler.
Specifically, let
$\omega_{k}^{(m)}=\Big{\\{}\lambda_{0Tk}^{(m)},\beta_{Tk}^{(m)},\lambda_{0Yk}^{(m)},\beta_{Yk}^{(m)},\eta_{k}^{(m)}\Big{\\}}$
denote the $m^{th}$ posterior draw. We sample from the joint one
block/variable at a time, conditional on the others, in the following way.
First, we initialize values for all parameters to $\omega_{k}^{(0)}$ for all
$k$ and. Then, at each iteration $m=1,2,\dots,M$, update the model parameters
for treatment courses $j=1,2,\dots,K$ sequentially. Using the usual bracket
notation to denotes posteriors, we begin by updating from conditional
posterior of parameters at $j=1$,
$\omega_{1}^{(m)}\sim[\omega_{1}\mid\omega_{2}^{(m-1)},\dots,\omega_{K}^{(m-1)}]$.
To obtain
$\omega_{1}^{(m)}=\\{\lambda_{0T1}^{(m)},\beta_{T1}^{(m)},\lambda_{0Y1}^{(m)},\beta_{Y1}^{(m)},\eta_{1}^{(m)}\\}$,
we update one component at a time conditional on the others. E.g. first,
update
$\lambda_{0T1}^{(m)}\sim[\lambda_{0T1}\mid\beta_{T1}^{(m-1)},\lambda_{0Y1}^{(m-1)},\beta_{Y1}^{(m-1)},\eta_{1}^{(m-1)}]$,
then update
$\beta_{T1}^{(m)}\sim[\beta_{T1}\mid\lambda_{0T1}^{(m)},\lambda_{0Y1}^{(m-1)},\beta_{Y1}^{(m-1)},\eta_{1}^{(m-1)}]$,
and so on. Having obtained $\omega_{1}^{(m)}$, we move to update
$\omega_{2}^{(m)}\sim[\omega_{2}\mid\omega_{1}^{(m)},\omega_{2}^{(m-1)},\dots,\omega_{K}^{(m-1)}]$
in the same way. Doing this for $j=1,2,\dots,K$ yields the $m^{th}$ posterior
draw
$\omega^{(m)}=\\{\omega_{1}^{(m)},\omega_{2}^{(m)},\dots,\omega_{K}^{(m)}\\}\sim[\omega_{1},\omega_{2},\dots,\omega_{K}]$.
Typically we run the sampler for some amount of “burn-in” iterations $M^{*}<M$
and base inferences only on draws $M^{*}<m\leq M$.

Since conjugate priors for the coefficients of the proportional hazard model
in (3) do not exist, we use a Metropolis step to update from their
conditionals. The efficiency of the Metropolis algorithm greatly depends on
the choice of proposal covariance, which can be difficult to tune manually in
our setting given the number of parameters and models. Instead, we use an
adaptive procedure (Haario et al., 2001) that tunes the proposal covariance in
the burn-in period. The sampler is initialized with a symmetric, diagonal
proposal covariance and run for $M^{*}$ iterations. At iteration $M^{*}$, we
compute the estimated covariance of the $M^{*}$ draws. The remaining draws
$M^{*}<m\leq M$ are then obtained via a proposal distribution with covariance
equal to a scale factor of this estimated covariance. Details and forms for
these conditionals are provided in the Supplement.

###  5.1 Posterior G-computation and Maximization

In this section, we describe posterior g-computation and optimization for
$\Psi^{r}(t)$ and utility functions more generally. In the Bayesian framework,
we conduct g-computation for each posterior draw of the parameters
$\\{\omega_{k}^{(m)}\\}_{k=1}^{K}$, for $m>M^{*}$. These parameters govern the
continuous-time transition process between treatment and death as well as the
time-varying confounder evolution over time. Given each draw of these
parameters, the idea is to simulate from the process under hypothetical
treatment rule $r$, then compute $\Psi^{r}(t)$ by averaging over those
simulations. Here, we describe this process for $K=2$ for simplicity.
Specifically, for the $m^{th}$ posterior draw, we do the following. For
$b=1,2,\dots,B$

  1. 1.

Simulate confounder at first treatment course, $L_{1}^{(b)}\sim
f_{1}(L_{1}\mid\eta_{1}^{(m)})$ and make a treatment decision
$d_{1}^{r,(b)}=r_{1}(L_{1}^{(b)})$

  2. 2.

Simulate time from first to second treatment:
$W_{Y1}^{(b)}\sim\lambda_{Y1}(w\mid
L_{1}^{(b)},A_{1}=d_{1}^{r,(b)}\lambda_{0Y1}^{(m)},\beta_{Y1}^{(m)})$

  3. 3.

Simulate time from first treatment to death
$W_{T1}^{(b)}\sim\lambda_{T1}(w\mid
L_{1}^{(b)},A_{1}=d_{1}^{r,(b)},\lambda_{0T1}^{(m)},\beta_{T1}^{(m)})$

  4. 4.

If $W_{T1}^{(b)}<W_{Y1}^{(b)}$, then set $T^{(b)}=W_{T1}^{(b)}$,
$\kappa^{(b)}=1$, and move to iteration $b+1$. Otherwise, store
$W_{Y1}^{(b)}$, set $\kappa^{(b)}=2$, and proceed to treatment course $k=2$:

  5. 5.

Simulate second course confounder $L_{2}^{(b)}\sim f_{2}(L_{2}\mid
L_{1}^{(b)},A_{1}=d_{1}^{r,(b)},W_{Y1}^{(b)},\eta_{2}^{(m)})$ and make a
decision
$d_{2}^{r,(b)}=r_{2}(\bar{L}_{2}^{(b)},A_{1}=d_{1}^{r,(b)},W_{Y1}^{(b)})$

  6. 6.

Simulate
$W_{T2}^{(b)}\sim\lambda_{T2}(w\mid\bar{L}_{2}^{(b)},\bar{A}_{2}=\bar{d}_{2}^{r,(b)},W_{Y1}^{(b)},\lambda_{0T2}^{(m)},\beta_{T2}^{(m)})$.
Set $T^{(b)}=W_{Y1}^{(b)}+W_{T2}^{(b)}$. That is, total time is the waiting
time between course one and two added with the waiting time from course two to
death.

After running for $B$ iterations, the $m^{th}$ posterior draw of the survival
curve can be computed $\Psi^{r}(t)^{(m)}=B^{-1}\sum_{b=1}^{B}I(T^{(b)}>t)$.
The sampling of waiting times in steps 2,3, and 6 can be done efficiently via
the inverse-CDF method under the Gamma Process models. For each $t$, this
yields a set of $M-M^{*}$ posterior draws for $\Psi^{r}(t)$. We can form a
point estimate by taking the mean across draws. Pointwise $100(1-\alpha)\%$
credible interval can be formed by using the $\alpha/2$ and $(1-\alpha/2)$
percentiles of draws as the endpoints.

Comparing survival under two different rules $r$ and $r^{\prime}$ does not
require re-running the MCMC sampler. They are functionals of the same
parameters, so we need only post-process the posterior draws,
$\\{\omega_{k}^{(m)}\\}_{k=1}^{K}$, under another rule $r^{\prime}$, as above,
to produce draws for $\Psi^{r^{\prime}}(t)$. In this way, we can also obtain
posterior draws of contrasts such as $\Psi^{r}(t)/\Psi^{r^{\prime}}(t)$ by
dividing the posterior draws $\Psi^{r}(t)^{(m)}/\Psi^{r^{\prime}}(t)^{(m)}$.
Absolute contrasts can be done by subtraction of draws. This is one of the
main benefits of our generative approach. Once posterior draws of the nuisance
parameters are obtained, they can be used to obtain posterior draws for
complex transformations without re-running the sampler.

In addition to effect estimation, optimization within a class of rules
$\mathcal{R}$ and can also be considered. For interpretability, we often
restrict to simple rules that are characterized by a low-dimensional set of
decision parameters $\tau\in\mathcal{T}$, the class of rules is then given by
$\mathcal{R}=\\{r_{\tau}:\tau\in\mathcal{T}\\}$. In this case, we denote
$r_{\tau}$ simply by $\tau$, since it completely determines the rule. In our
setting, this is suitable since rules of interest are governed by a few
thresholds for EF. The task is then to search $\mathcal{T}$ for the rule
parameters that maximize the probability of surviving past some time, had we
treated according to a rule with those parameters, $\Psi^{\tau}(t)$. Note that
$\Psi^{\tau}(t)=\Psi^{\tau}(t;\omega)$ is a function of the unknown model
parameters, $\omega$, written explicity as such for emphasis. Thus a posterior
on these model parameters induces a posterior on the optimal rule parameters,
$\tau^{*}(\omega)=\operatorname*{argmax}_{\tau\in\mathcal{\mathcal{T}}}\Psi^{\tau}(t;\omega)$.
Computationally, maximization must be done for each posterior draw
$\omega^{(m)}$, which yields a posterior draw of the optimal parameters
$\tau^{*}(\omega^{(m)})\sim f(\tau^{*}(\omega)|\mathcal{D})$. In our
application, $\tau$ is a pair of thresholds - one for relative change in EF
and another for absolute EF - while $\mathcal{T}$ is specified as a two-
dimensional grid of possible thresholds. This relatively small search space
permits an exhaustive search over the grid. For each $\omega^{(m)}$,
g-computation is done for each $\tau$ in the grid we choose
$\tau^{*}(\omega^{(m)})$ to be the threshold with the highest survival
probability. This consititutes a draw from the bivariate posterior mass
function (pmf), $f(\tau^{*}(\omega)|\mathcal{D})$. The posterior mode across
the $m$ draws can be taken as a point estimate of $\tau^{*}$. A
$100(1-\alpha)\%$ credible set of threshold pairs,
$C(1-\alpha)\subset\mathcal{T}$, can be constructed by lowering a flat plane
onto the pmf from above until $100(1-\alpha)\%$ of the mass lies above the
plane. $C(1-\alpha)$ contains the pairs with mass above this plane. This is a
bivariate analogue of the usual highest posterior density interval (HDI).

While our interest is primarily in overall survival, our generative approach
allows for more general utility function outcomes (Murray et al., 2017) that
directly capture potential risk-benefit tradeoffs between survival and
cardiotoxicity. For instance, consider $U^{r}(s,t)=\Psi^{r}(t)-\Phi^{r}(s)$,
where $\Phi^{r}(s)=P(\min_{k>2}L_{k}(r)<s)$. If $L$ represents EF, then
$\Phi^{r}(s)$ is the probability that EF drops below threshold $s$ after
initial treatment. In the clinical literature, $s=.50$ is one common
definition of an adverse cardiac event that we would ideally want to avoid.
This function deducts utility based on the cardiac risk of the rule and can be
computed just as above except at the end we compute
$U^{r}(s,t)^{(m)}=\frac{1}{B}\sum_{b=1}^{B}I(T^{(b)}>t)-I(L_{2}^{(b)}<s)I(\kappa^{(b)}=2)$.
We could consider more general rules that weigh the risk and benefit
components differently. The idea is that rules can then be contrasted and
optimized on the basis this utility. By construction, both the benefit
($\Psi^{r}(t)$) and the risk ($\Phi^{r}(s)$) are on the same (probability)
scale, which eases interpretation. For instance, a rule where the risk exceeds
benefit ($\Phi^{r}(s)>\Psi^{r}(t)$) would have negative utility. If
$r^{\prime}$ has the same survival probability as $r$, but comes with lower
risk, then it will have higher utility and would be preferred.

##  6 Simulation Experiments

Since this is a complex scenario with intricate treatment and time-varying
confounding mechanisms, we ran extensive simulation studies evaluating the
proposed method’s ability to capture the survival time distribution in
situations that mirror our application. Thus, the settings were carefully
chosen to match our data: First, as in our data, there are four treatment
courses. In a typical simulated data set, about $60\%$ of subjects complete
the full sequence (the rest are either censored or die before completing) and
about 40% of the patients die while the remaining are censored. In the data,
the corresponding proportions are about $60\%$ and $40\%$ as well. For each
subject, we simulate two time-varying confounders - one continuous and one
binary - as well as four time-constant confounders. This matches our data
application which has a continuous and binary time-varying confounder (EF and
presence of bloodstream infection, respectively) as well as several time-
constant factors. At each course, we simulate treatment conditional on time-
constant confounders, current and previous values of time-varying confounders,
and previous treatment. Then, we simulate the waiting time until the next
course, death, or censoring from separate Weibull proportional hazard models.
Current and previous treatment, current and previous time-varying confounders,
time-constant confounders, and the time since the previous course initiation
all affect the waiting times until the next event by multiplying the Weibull
baseline hazards. Longer times since last treatment are set to increase the
time to next treatment, while decreasing time to death. This induces the
positive correlation in Figure 2 and encodes the idea that longer waiting
times between treatments indicate worse survival prospects. We simulate 1,000
such data sets in a large sample setting (n=1000) and small sample setting
(n=300). For each simulated data set, we estimate posterior mean and 95%
credible intervals for $\Psi^{r}(t)$ at pre-specified $t$, with
$r(L_{k1})=I(L_{k1}<0)$. Here $L_{k1}$ is the continuous time-varying
confounder with sufficient support around $0$, and $r$ assigns treatment
$A_{k}=1$ if $L_{k1}<0$. The proposed approach is compared with a parametric
version where the waiting times are modeled via a correctly specified Bayesian
Weibull hazard model. The results are presented in Table 1.

Table 1: Simulation Results: Absolute bias of posterior mean of $\Psi^{r}(t)$,
for various $t$ (as percentage of the truth). Coverage of 95% credible
interval for $\Psi^{r}(t)$. Average width of 95% credible interval.

|  | Gamma Process | Weibull
---|---|---|---
Setting | $t$ | Bias (%) | Coverage | Int. Width | Bias (%) | Coverage | Int. Width
N=1000 | 5 | 0.07 | 95.4 | 0.04 | 0.14 | 94.5 | 0.04
10 | 0.03 | 95.2 | 0.07 | 0.22 | 94.8 | 0.07
15 | 2.29 | 90.4 | 0.08 | 0.24 | 95.8 | 0.08
20 | 185.00 | 18.8 | 0.13 | 6.58 | 94.5 | 0.06
N=300 | 5 | 0.34 | 94.9 | 0.07 | 0.53 | 95.8 | 0.07
10 | 0.13 | 94.7 | 0.12 | 0.98 | 95.2 | 0.12
15 | 4.52 | 89.3 | 0.14 | 0.51 | 94.9 | 0.14
20 | 227.00 | 14.7 | 0.18 | 25.00 | 90.3 | 0.11

Across simulations, the proportion of patients still at risk for subsequent
treatment or death at time points 5, 10, 15, and 20 are about 95%, 85%, 60%,
and 2%, respectively. Accordingly, at early time points, $t\in\\{5,10,15\\}$
the bias is quite low: .07%, .03%, and 2.29%, respectively, in the $N=1000$
setting. Similarly, the coverage of the 95% credible intervals are close to
nominal. At $t=20$, only about 2% of subjects are left in the sample in a
typical simulation run. Thus, the quality of the semiparametric estimate is
degraded. With such low proportion at risk, all flexible estimators will
perform poorly at the tails. The Weibull model is correctly specified in this
setting, so it does better in terms of bias and coverage even at the tail. In
practical applications, however, a Weibull model is unlikely to be correct in
applications and we would avoid extrapolating so far into the tail. These
results demonstrate that GP-based posterior estimates have comparable
frequentist properties to the correctly specified model in a setting with
complex time-varying confounding and informative treatment timing. Despite
being unspecified, the GP posterior picks up the Weibull baseline hazard form.
The small sample setting with $N=300$ is more consistent with our data
application and exhibits largely the same performance. As expected, interval
widths tend to be larger for both models relative to the larger sample case.
Similarly, bias tends to be higher relative to the large sample setting as
well. Additional details on model specification and other settings can be
found in the Supplement. Because the GP models require a partitioning of the
time interval, the appendix includes a sensitivity for the $n=1000$ setting
where the number of partitions is halved, yielding a coarser partition. The
results are generally consistent with only a slight increase in bias and
coverage (1% and 96%, respectively)

##  7 Evaluating Ejection-Fraction Based Treatment Rules

We analyze data from the AAML1031-PHIS cohort with n=292 subjects who had a
median followup time of 2.5 years and underwent a maximum of four ACT
treatment courses. While a small sample, there are only about 500 new
diagnoses of pediatric AML annually in the United States. Table 2 displays
some sample statistics for our data. About 36% of the patients either died or
were censored before completing the four-course treatment protocol. Death is
observed for about 40% of patients. As illustrated in Figure 2, the
distributions of the waiting times between treatments tend to be centered
around 35 days, with some right-skewness due to subjects who take longer than
usual to recover from their previous course. This motivates the need for
flexible models for these waiting times.

Table 2: Summary Statistics: Median and interquartile range presented for continuous features. Counts and proportions presented for discrete features.  Features | $N=292$
---|---
Follow-up (years) | 2.6 (1.3-3.7 )
Death | 114 (40%)
Num. Treatment Courses, ($\kappa$) |
1 | 22 (8%)
2 | 36 (12%)
3 | 46 (16%)
4 | 188 (63%)
AML Risk Classification (high) | 64 (23%)
WBC (cells/uL) | 20 (6.8-65)
Male | 152 (52%)

We also observe several important patient-level confounders. Time-varying
confounders include EF (recorded ahead of each treatment decision) and the
(binary) presence/absence of a bloodstream infection since the start of the
previous treatment course, which we denote as $V_{k}$. The treatment options
of interest consist of chemotherapy treatments that include ($A_{k}=1$) or
exclude ACT ($A_{k}=0$), which is recorded for each visit. Each of the two
possible treatments are feasible in courses 1,2, and 4 for all histories. That
is $\mathcal{S}_{k}(h_{k})=\\{0,1\\}$ for all $h_{k}$ at $k=1,2,4$. However,
in general subjects were not eligible for ACT at course three, meaning
$\mathcal{S}_{3}(h_{3})=\\{0\\}$, $\forall h_{3}$. Several potential baseline
confounders including age, sex, baseline white blood cell (WBC) count, and AML
risk classification are also recorded. Baseline white blood cell count
indicates amount of leukemia at the start of the treatment protocol and may
inform subsequent ACT decisions and survival. AML risk classification (RS) is
a binary covariate (high versus low) based on favorable/unfavorable molecular
and cytogenetic features linked with prognosis of worse survival. We assume
sequential ignorability and non-informative censoring (as well as the other
identification assumptions discussed) all hold conditional on these factors.
While these assumptions cannot be verified in the data, we know that these
covariates (EF in particular) play an influential role in treatment decisions.

Using the notation of our model, $L_{k}=(EF_{k},V_{k},age,sex,WBC,RS)$. For a
subject reaching course $k$, the available history is
$H_{k}=(\bar{A}_{k-1},\bar{W}_{k-1},\bar{EF}_{k},\bar{V}_{k},age,sex,WBC,RS)$.
We fit the models described in Section 4, adjusting for information in current
and previous treatment course in all models. Since $EF_{k}\in(0,1)$ is a
proportion, at each course we model its distribution as a Beta with
conditional mean being a function of history from course $k-1$ through a logit
link. Similarly, we model $V_{k}$ with a sequence of logistic regressions to
respect the binary domain. At each course, the hazard models for subsequent
death or treatment conditions on current EF, $V$, and ACT decisions as well as
their values from the previous treatment course in addition to all baseline
confounder. To avoid linearity assumptions, the time since previous treatment,
$W_{k-1}$, is included as a categorical variable with four levels with cutoffs
at the 25th, 50th, and 75th percentiles of the observed values. Detailed
specifications of these models are provided in the Supplement. We use the
outlined MCMC approach to obtain draws from the posterior of all the unknowns.

Figure 3:  Left: Posterior inference for $\Psi^{r}(t)$ under GP versus Weibull
models for $r$. Red ticks indicate median time of treatment sequence 2, 3, and
4 initiation in order. Middle: Posterior inference for $\Psi^{r^{\prime}}(t)$
under GP versus Weibull models for $r$. Right: risk ratios comparing rule $r$
with $r^{\prime}$.

Current criteria for ACT withholding varies between trials and thresholds have
not been rigorously validated. For instance, one rule that has been
recommended in the clinical literature is to withhold ACT ($A_{k}=0$) if
current EF has declined $10\%$ from baseline to below $50\%$. I.e.
$A_{k}=r_{k}(EF_{k},EF_{1})$, where
$r_{k}(EF_{k},EF_{1})=1-I(EF_{k}/EF_{1}-1<-.1)I(EF_{k}<.5)$ (Neuendorff et
al., 2020). To assess the impact of this rule, the left panel of Figure 3
visualizes posterior inference for survival, $\Psi^{r}(t)$, under rule
$r(\bar{E}F_{K})=(A_{1}=1-I(EF_{1}<.5),A_{2}=r_{2}(EF_{2},EF_{1}),A_{3}=0,A_{4}=r_{4}(EF_{4},EF_{1}))$.
At $k=3$, the rule always withholds ACT since it is not feasible at the third
course. Figure 3 presents results under both the outlined GP model and a
parametric Weibull model. Note that the first treatment is assigned only using
baseline EF since percentage change from itself is zero. It can be seen from
Figure 3 that the parametric model smooths over intricacies of the survival
function, whereas the GP captures important features - for instance, the steep
drop-off in survival around month 4. In the middle panel, we present the
posterior survival curve for rule
$r^{\prime}(\bar{E}F_{K})=(A_{1}=0,A_{2}=1-I(EF_{2}/EF_{1}-1<-.1),A_{3}=0,A_{4}=1-I(EF_{4}/EF_{1}-1<-.1))$.
This starts everyone on nACT, then withholds ACT only if relative change in EF
decreases more than 10% (regardless of absolute EF level). This rule relaxes
the criterion for withholding ACT begins by initially withholding ACT. The
right panel presents the ratio in survival probabilities for rule $r$ versus
$r^{\prime}$ and shows relative survival benefit for $r$, which has stricter
criteria for withholding ACT. In summary, posterior mean 3-year survival
probability under $r$ is $.64$ $[95\%\text{CrI}:.58-.69]$. Under rule
$r^{\prime}$ posterior mean 3-year survival probability is $.60$
$[95\%\text{CrI}:.46-.7]$. The right panel of Figure 3 shows the relative
survival probabilities under each rule at various time points. There is higher
survival under rule $r$ across time.

The actual thresholds of $.5$ and $-.1$ in $r$ are still debated in the
literature. To inform this debate, we consider finding thresholds that
optimize 3-year survival, an important endpoint in pediatric AML.
Specifically, we define a class of rules
$r(\bar{E}F_{K};\tau)=(A_{1}=1-I(EF_{1}<\tau_{2}),A_{2}=r_{2}(EF_{2},EF_{1};\tau),A_{3}=0,A_{4}=r_{k}(EF_{4},EF_{1};\tau))$.
Here,
$r_{k}(EF_{k},EF_{1};\tau)=1-I(EF_{k}/EF_{1}-1<\tau_{1})I(EF_{k}<\tau_{2})$
and $\tau=(\tau_{1},\tau_{2})$ parameterizes the rule. For
$\tau\in\mathcal{\mathcal{T}}=\\{0,-.1,-.2,-.3,-.4,-.5\\}\times\\{.4,.5,.6,.7,.8,.9\\}$,
this defines a class of rules, where the range of threshold values are
selected with clinical guidance and positivity constraints in mind. We
optimize this rule by finding a posterior over
$\tau^{*}=\operatorname*{argmax}_{\tau\in\mathcal{\mathcal{T}}}U^{\tau}(36,.5)$,
where $U^{\tau}(36,.5)$ is the 3-year utility had we treated according to rule
$r(\bar{E}F_{K};\tau)$. As discussed in Section 5.1, this utility is the risk
of post-baseline EF dipping below .5 subtracted from the 3-year survival
probability. We do this as outlined in Section 5.1 and display the posterior
draws of $\tau^{*}$ as points in the right panel over Figure 4. The points are
jittered and transparent to create a kind of heat map with more posterior
draws leading to darker clouds at that $\tau$ pair. The posterior mode
threshold pair is $(\tau_{1}=0,\tau_{2}=.7)$. This favors a rule that
withholds ACT if there is any decline from baseline, as long as EF falls below
.7. If a patient experiences an increase in EF from baseline, then ACT is
assigned. This makes clinical sense as ACT is cardiotoxic and we may wish to
assign it if heart function is improved. The red points in the figure indicate
all threshold pairs that are in the 90% credible set,
$C(.90)\subset\mathcal{T}$, as discussed in Section 5.1. Note that the
previously discussed rule of thumb with $(\tau_{1}=-.1,\tau_{2}=.5)$ is in
this 90% credible set and, in that sense, is justified by this data. The left
panel of Figure 4 are draws from the posterior of the 3-year survival
probability under posterior mode rule, $(\tau_{1}=0,\tau_{2}=.7)$. The
posterior distribution of the utility is shifted left since there is about a
10% posterior mean risk of EF dipping below .5 under this rule.

This approach characterizes the optimal rule as an unknown about which we have
uncertainty. There are many ways this can be translated to inform clinical
decision-making in practice. One possibility is outputting all rules in the
credible set $C(.90)$ and leaving the final decision up to the clinician. The
credible set gives some guidance about how well-supported various rules are by
the data. Another option is to provide the entire list of thresholds, sorted
by the posterior probability of each threshold pair being the optimal. The
clinician may be wary of choosing pairs at the bottom of the list and more
comfortable choosing ones higher on the list. In either case, we view this
output as a decision support tool that allows clinicians to gauge how strongly
the data support various treatment options. The ultimate decision is made by
the clinician and patient together.

Figure 4:  Left: Posterior distribution of 3-year survival and utility under
posterior mode optimal rule, $(\tau_{1}=0,\tau_{2}=.7)$. Under this rule,
posterior mean 3-year survival is about 60% with about 10% risk of EF dipping
below .5 during the treatment course. Right: Posterior draws of the optimal
rule $\tau^{*}$. Red points indicate threshold pairs in the 90% credible set,
$C(1-\alpha)\subset\mathcal{T}$ for $\alpha=.10$

##  8 Discussion

Approaches for estimating causal effects of dynamic treatment rules exist, but
our data presented several unique challenges: a fixed number of treatment
courses that occur at random, informative times and occurrence of death and
censoring that prevent subjects from completing the full sequence of treatment
courses. We use results from causal inference to identify a causal survival
probability under a particular dynamic treatment rule and specify a generative
Bayesian semiparametric model for how patients transition from one treatment
to subsequent treatment (or death) conditional on treatment and confounder
history. An appealing aspect of this approach is full posterior inference for
all functions of the transition process - including causal survival
probabilities, contrasts of survival probabilities under alternate rules, and
even posterior inference over optimal rules. In simulations we assess the
frequentist performance of this Bayesian estimator and report nominal coverage
and low bias for the 95% credible interval and posterior mean survival
probability, respectively. In our analysis of data from the AAML trial, we
show that our semiparametric model captures complexities in the survival curve
that the parametric Weibull model smooths over. We also illustrate that
current rules of thumb for ACT administration fall within a 90% credible set
of optimal rules - indicating that current practice has some support in the
data.

A central difficulty with time-varying treatments is the difficulty in
estimating the effect of all possible treatment sequences in small samples.
Suppose there are four treatment courses with two possible treatments in each
course leading to $2^{4}$ possible treatment sequences. In small samples there
may be very few patients who follow some of these sequences, essentially
leading to random violations of positivity. An appealing aspect of this
Bayesian approach is that the credible intervals on the estimated survival
probabilities under such sequences are quite wide - reflecting our increased
uncertainty due to sparsity of data. For instance, when we try to estimate the
effect of a rule that may assign ACT in course three, the 95% credible
interval around 3-year survival probability covers almost the entire range
$[0,1]$. Due to the fact that ACT is not feasible (only two subjects were
assigned ACT in the third course as exceptions to protocol), our credible
intervals widen to reflect the increased uncertainty that stems from
extrapolating beyond the range of the observed data.

### Acknowledgements and Funding

This work was partially funded by: 1) the Patient-Centered Outcomes Research
Institute (PCORI) contract ME-2021C3-24942; 2) National Institutes of Health K
Award K01HL143153; 3); National Institutes of Health R01 Award R01HL164925. 4)
Brown University Office of the Vice President for Research (OVPR Salomon
Faculty Research Award awarded.

## References

  * Chipman et al. (2010) Chipman, H. A., E. I. George, and R. E. McCulloch (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics 4(1), 266–298.
  * Getz et al. (2019a) Getz, K. D., L. Sung, B. Ky, R. B. Gerbing, K. J. Leger, A. B. Leahy, L. Sack, W. G. Woods, T. Alonzo, A. Gamis, and R. Aplenc (2019a).  Occurrence of treatment-related cardiotoxicity and its impact on outcomes among children treated in the aaml0531 clinical trial: A report from the children’s oncology group.  Journal of Clinical Oncology 37(1), 12–21.  PMID: 30379624.
  * Getz et al. (2019b) Getz, K. D., L. Sung, B. Ky, R. B. Gerbing, K. J. Leger, A. B. Leahy, L. Sack, W. G. Woods, T. Alonzo, A. Gamis, and R. Aplenc (2019b).  Occurrence of treatment-related cardiotoxicity and its impact on outcomes among children treated in the aaml0531 clinical trial: A report from the children’s oncology group.  Journal of Clinical Oncology 37(1), 12–21.  PMID: 30379624.
  * Guan et al. (2020) Guan, Q., B. J. Reich, E. B. Laber, and D. Bandyopadhyay (2020).  Bayesian nonparametric policy search with application to periodontal recall intervals.  Journal of the American Statistical Association 115(531), 1066–1078.  PMID: 33012901.
  * Haario et al. (2001) Haario, H., E. Saksman, and J. Tamminen (2001).  An adaptive Metropolis algorithm.  Bernoulli 7(2), 223 – 242.
  * Hua et al. (2021) Hua, W., H. Mei, S. Zohar, M. Giral, and Y. Xu (2021).  Personalized Dynamic Treatment Regimes in Continuous Time: A Bayesian Approach for Optimizing Clinical Decisions with Timing.  Bayesian Analysis, 1 – 30.
  * Kalbfleisch (1978) Kalbfleisch, J. D. (1978).  Non-parametric bayesian analysis of survival time data.  Journal of the Royal Statistical Society. Series B (Methodological) 40(2), 214–221.
  * Kalbfleisch and Prentice (2011) Kalbfleisch, J. D. and R. L. Prentice (2011).  The statistical analysis of failure time data.  John Wiley & Sons.
  * Murphy (2003) Murphy, S. A. (2003).  Optimal dynamic treatment regimes.  Journal of the Royal Statistical Society. Series B (Statistical Methodology) 65(2), 331–366.
  * Murray et al. (2017) Murray, T. A., P. F. Thall, Y. Yuan, S. McAvoy, and D. R. Gomez (2017).  Robust treatment comparison based on utilities of semi-competing risks in non-small-cell lung cancer.  Journal of the American Statistical Association 112(517), 11–23.  PMID: 28943681.
  * Neuendorff et al. (2020) Neuendorff, N. R., K. P. Loh, A. S. Mims, K. Christofyllakis, W.-K. Soo, B. Bölükbasi, C. Oñoro-Algar, W. G. Hundley, and H. D. Klepin (2020, 02).  Anthracycline-related cardiotoxicity in older patients with acute myeloid leukemia: a Young SIOG review paper.  Blood Advances 4(4), 762–775.
  * 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(9), 1393–1512.
  * Tsiatis et al. (2019) Tsiatis, A. A., M. Davidian, S. T. Holloway, and E. B. Laber (2019).  Dynamic treatment regimes: Statistical methods for precision medicine.  Chapman and Hall/CRC.
  * Xu et al. (2016) Xu, Y., P. Müller, A. S. Wahed, and P. F. Thall (2016).  Bayesian nonparametric estimation for dynamic treatment regimes with sequential transition times.  Journal of the American Statistical Association 111(515), 921–950.  PMID: 28018015.


# Universal solution for optimal protocols of finite-time and weak processes

Pierre Nazé<EMAIL_ADDRESS>Departamento de Física, Instituto de
Geociências e Ciências Exatas, Universidade Estadual Paulista “Júlio de
Mesquita Filho”, 13506-900, Rio Claro, SP, Brazil

###### Abstract

The analytical expression for the optimal protocol of the thermodynamic work
and its variance for finite-time, isothermal and weak processes is presented.
The method that solves the Euler-Lagrange integral equation is quite general
and depends only on the properties of time-reversal symmetry of the optimal
protocol. The solution is proven to be physically consistent and many examples
are solved to illustrate the method. To overcome the hypothesis consistency
problem of the singular part of the solution, an interpretation of the
appearance of the delta peaks and their derivatives is presented.

##  I Introduction

In times when computational endeavor has achieved its maximum power, a
detriment in the search for simple analytical methods is a severe symptom.
Scientific analyses restrain themselves in case-to-case examples and
generalizations, which without a doubt is an important part of every
scientific investigation, become a herculean task. That is the case for
instance of problems that depend on methods to solve integral equations, which
rely basically on the analysis of types of kernels and equations [1].

In such a scenario, it was quite disappointing to attack the problem of
finding the optimal protocols of the thermodynamic work and its variance in
finite-time and weak processes [2], since its Euler-Lagrange integral equation
was of the worst type possible – Fredholm integral equation of first kind and
constant limits of integration –, and highly depend on the kernel proposed
[3]. Alternative methods, like genetic programming, although they present a
result with some easiness, predicted in some cases optimal protocols outside
the range of linear-response theory, which would defy its utility for this
case.

This work tries to solve all the complications described above by solving the
Euler-Lagrange integral equation of the referred problem in a general way. The
method is quite general and depends only on the symmetric property of the
optimal protocol, which was already found in previous work [3]. I find that
the optimal protocol, for any reasonable physical system in the regime
proposed, is universal, being composed of a straight line with jumps at the
beginning and end of the process, and a series of Dirac deltas and their
derivatives. Also, it depends exclusively on the ratio between the switching
time of the process and the relaxation timescale of the system. To illustrate
the method, I present many examples, from Brownian motion to the Ising chain,
whose results corroborate the approach. Finally, to overcome the problem if
the singular part of the equation is in accord with the hypothesis of linear
response theory, an interpretation of the appearance of the Dirac deltas and
their derivatives is presented.

##  II Preliminaries

I start defining notations and developing the main concepts to be used in this
work.

Consider a classical system with a Hamiltonian
$\mathcal{H}(\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t)),\lambda(t))$,
where $\mbox{\boldmath$z$}(\mbox{\boldmath$z_{0}$},t)$ is a point in the phase
space $\Gamma$ evolved from the initial point $z_{0}$ until time $t$, with
$\lambda(t)$ being a time-dependent external parameter. During a switching
time $\tau$, the external parameter is changed from $\lambda_{0}$ to
$\lambda_{0}+\delta\lambda$, with the system being in contact with a heat bath
of temperature $\beta\equiv{(k_{B}T)}^{-1}$, where $k_{B}$ is Boltzmann’s
constant. The average work performed on the system during this interval of
time is

$\overline{W}\equiv\int_{0}^{\tau}\left\langle\overline{\partial_{\lambda}\mathcal{H}}(t)\right\rangle_{0}\dot{\lambda}(t)dt,$
(1)

where $\partial_{\lambda}$ is the partial derivative in respect to $\lambda$
and the superscripted dot the total time derivative. The generalized force
$\left\langle\overline{\partial_{\lambda}\mathcal{H}}\right\rangle_{0}$ is
calculated using the averaging $\overline{\cdot}$ over the stochastic path and
the averaging $\langle\cdot\rangle_{0}$ over the initial canonical ensemble.
The external parameter can be expressed as

$\lambda(t)=\lambda_{0}+g(t)\delta\lambda,$ (2)

where, to satisfy the initial conditions of the external parameter, the
protocol $g(t)$ must satisfy the following boundary conditions

$g(0)=0,\quad g(\tau)=1.$ (3)

We consider as well that $g(t)\equiv g(t/\tau)$, which means that the
intervals of time are measured according to the switching time unit.

Linear-response theory aims to express average quantities until the first-
order of some perturbation parameter considering how this perturbation affects
the observable to be averaged and the process of average [4]. In our case, we
consider that the parameter does not considerably changes during the process,
$|g(t)\delta\lambda/\lambda_{0}|\ll 1$, for all $t\in[0,\tau]$. In that
manner, using such framework, the generalized force can be approximated until
first-order as

$\begin{split}\left\langle\overline{\partial_{\lambda}\mathcal{H}}(t)\right\rangle_{0}=&\,\left\langle\partial_{\lambda}\mathcal{H}\right\rangle_{0}+\delta\lambda\left\langle\partial_{\lambda\lambda}^{2}\mathcal{H}\right\rangle_{0}g(t)\\\
&-\delta\lambda\int_{0}^{t}\phi_{0}(t-t^{\prime})g(t^{\prime})dt^{\prime}.\end{split}$
(4)

The quantity $\phi_{0}(t)$ is the so-called response function [4], which can
be conveniently expressed as the derivative of the relaxation function
$\Psi_{0}(t)$ [4]

$\phi_{0}(t)=-\frac{d\Psi_{0}}{dt}.$ (5)

In our particular case, the relaxation function is calculated as

$\Psi_{0}(t)=\beta\left\langle\partial_{\lambda}\mathcal{H}(0)\overline{\partial_{\lambda}\mathcal{H}}(t)\right\rangle_{0}-\mathcal{C},$
(6)

where the constant $\mathcal{C}$ is calculated to vanish the relaxation
function for long times [4]. We define as the relaxation timescale of the
system the quantity

$\tau_{R}=\int_{0}^{\infty}\frac{\Psi(t)}{\Psi(0)}dt.$ (7)

The generalized force, written in terms of the relaxation function, can be
expressed as

$\begin{split}\left\langle\overline{\partial_{\lambda}\mathcal{H}}(t)\right\rangle_{0}=&\,\left\langle\partial_{\lambda}\mathcal{H}\right\rangle_{0}-\delta\lambda\widetilde{\Psi}_{0}g(t)\\\
&+\delta\lambda\int_{0}^{t}\Psi_{0}(t-t^{\prime})\dot{g}(t^{\prime})dt^{\prime},\end{split}$
(8)

where
$\widetilde{\Psi}_{0}\equiv\Psi_{0}(0)-\left\langle\partial_{\lambda\lambda}^{2}\mathcal{H}\right\rangle_{0}$.
Finally, combining Eqs. (1) and (8), the average work performed at the linear
response of the generalized force is

$\begin{split}\overline{W}=&\,\delta\lambda\left\langle\partial_{\lambda}\mathcal{H}\right\rangle_{0}-\frac{\delta\lambda^{2}}{2}\widetilde{\Psi}_{0}\\\
&+\delta\lambda^{2}\int_{0}^{\tau}\int_{0}^{t}\Psi_{0}(t-t^{\prime})\dot{g}(t^{\prime})\dot{g}(t)dt^{\prime}dt.\end{split}$
(9)

We observe that the double integral on Eq. (9) vanishes for long switching
times [5]. Therefore the other terms are part of the contribution of the
difference of free energy, since this quantity is exactly the average work
performed for quasistatic processes in isothermal drivings. Thus, we can split
the average work into the difference of free energy $\Delta F$ and
irreversible work $W_{\rm irr}$

$\Delta
F=\delta\lambda\left\langle\partial_{\lambda}\mathcal{H}\right\rangle_{0}-\frac{\delta\lambda^{2}}{2}\widetilde{\Psi}_{0},$
(10)
$\begin{split}W_{\text{irr}}=\delta\lambda^{2}\int_{0}^{\tau}\int_{0}^{t}\Psi_{0}(t-t^{\prime})\dot{\lambda}(t^{\prime})\dot{\lambda}(t)dt^{\prime}dt.\end{split}$
(11)

In particular, the irreversible work can be rewritten using the symmetric
property of the relaxation function [4]

$\begin{split}W_{\rm
irr}=\frac{\delta\lambda^{2}}{2}\int_{0}^{\tau}\int_{0}^{\tau}\Psi_{0}(t-t^{\prime})\dot{\lambda}(t^{\prime})\dot{\lambda}(t)dt^{\prime}dt.\end{split}$
(12)

I establish at this point the regimes where linear-response theory is able to
describe thermodynamic processes. Those regimes are determined by the relative
strength of the driving with respect to the initial value of the protocol,
$\delta\lambda/\lambda_{0}$, and the rate by which the process occurs with
respect to the relaxation time of the system, $\tau_{R}/\tau$. See Fig. 1 for
a diagram depicting the regimes. In region 1, the so-called slowly-varying
processes, the ratio $\delta\lambda/\lambda_{0}$ is arbitrary, while
$\tau_{R}/\tau\ll 1$. By contrast, in region 2, the so-called finite-time and
weak processes, the ratio $\delta\lambda/\lambda_{0}\ll 1$, while
$\tau_{R}/\tau$ is arbitrary. In region 3, the so-called arbitrarily far-from-
equilibrium processes, both ratios are arbitrary. Linear-response theory can
only describe regions 1 and 2 [5]. In this work, we are going to focus on
region 2 only.

Figure 1: (Color online) Diagram of nonequilibrium regions. Region 1: slowly-
varying processes, Region 2: finite-time but weak processes and Region 3:
arbitrarily far-from-equilibrium processes. Linear response theorem can
describe regions 1 and 2.

Consider the irreversible work rewritten in terms of the protocols $g(t)$
instead of its derivative

$\displaystyle W_{\rm irr}=$
$\displaystyle\frac{\delta\lambda^{2}}{2}\Psi(0)+\delta\lambda^{2}\int_{0}^{\tau}\dot{\Psi}_{0}(\tau-t)g(t)dt$
(13)
$\displaystyle-\frac{\delta\lambda^{2}}{2}\int_{0}^{\tau}\int_{0}^{\tau}\ddot{\Psi}(t-t^{\prime})g(t)g(t^{\prime})dtdt^{\prime}.$
(14)

Using calculus of variations, we can derive the Euler-Lagrange equation that
furnishes the optimal protocol $g^{*}(t)$ of the system that will minimize the
irreversible work [3]

$\int_{0}^{\tau}\ddot{\Psi}_{0}(t-t^{\prime})g^{*}(t^{\prime})dt^{\prime}=\dot{\Psi}_{0}(\tau-t).$
(15)

In particular, the optimal irreversible work will be [3]

$W_{\rm
irr}^{*}=\frac{\delta\lambda^{2}}{2}\Psi(0)+\frac{\delta\lambda^{2}}{2}\int_{0}^{\tau}\dot{\Psi}_{0}(\tau-t)g^{*}(t)dt.$
(16)

Also, the Euler-Lagrange equation (15) furnishes also the optimal protocol
that minimizes the variance of the work [2]. In this case, the optimal
variance of work is

$\sigma_{\rm
W}^{2^{*}}=\frac{\beta\delta\lambda^{2}}{4}\Psi(0)+\frac{\beta\delta\lambda^{2}}{4}\int_{0}^{\tau}\dot{\Psi}_{0}(\tau-t)g^{*}(t)dt.$
(17)

The objective of this work is to solve the Euler-Lagrange equation (15) for
any kind of relaxation function of isothermal processes. Mathematically
speaking, I will present a new method to solve Fredholm integral equations of
the first type for symmetric kernels. Such a method relies only on the time-
reversal symmetry property of the solution.

##  III Universal solution

To derive the universal solution of Eq. (15), I will basically use the
symmetric property of the optimal protocol

$g^{*}(t)=1-g^{*}(\tau-t).$ (18)

First, I open the right-hand side of Eq. (15) in the appropriate integrals

$\int_{0}^{t}\ddot{\Psi}(t-t^{\prime})g^{*}(t^{\prime})dt^{\prime}+\int_{t}^{\tau}\ddot{\Psi}(t-t^{\prime})g^{*}(t^{\prime})dt^{\prime}=\dot{\Psi}(\tau-t)$
(19)

Using the symmetric property (18) in the second term of the right-hand side of
Eq. (19), one can show that

$\int_{0}^{t}\ddot{\Psi}(t-t^{\prime})g^{*}(t^{\prime})dt^{\prime}=\int_{0}^{\tau-t}\ddot{\Psi}(\tau-
t-t^{\prime})g^{*}(t^{\prime})dt^{\prime}.$ (20)

Therefore, Eq. (20) must be equal to a symmetric function $h(t)$, such that,
$h(t)=h(\tau-t)$. Using the solution of the Euler-Lagrange equation, the above
result is equal to

$\int_{0}^{t}\ddot{\Psi}(t-t^{\prime})g^{*}(t^{\prime})dt^{\prime}=\int_{0}^{t}\ddot{\Psi}(t-t^{\prime})(1-g^{*}(\tau-t^{\prime}))dt^{\prime},$
(21)

which is another way to express the symmetry of the optimal protocol.
Consider, by simplicity, a function $g_{0}(t)$, such that

$\int_{0}^{t}\ddot{\Psi}(t-t^{\prime})g_{0}(t^{\prime})dt^{\prime}=C_{0},$
(22)

where $C_{0}$ is a constant in time. Applying the convolution theorem we have

$g_{0}(t)=\mathcal{L}^{-1}_{s}\left\\{\frac{C_{0}}{s\mathcal{L}_{t}\\{\ddot{\Psi}(t)\\}(s)}\right\\}(t),$
(23)

where $\mathcal{L}_{t}\\{\cdot\\}$ and $\mathcal{L}_{s}^{-1}\\{\cdot\\}$ are
respectively the Laplace and inverse Laplace transform. Assuming that
$\mathcal{L}_{t}\\{\Psi(t)\\}(s)$ can be expressed as a Taylor series in
$s=0$, we have by applying Horner-Ruffini method

$\frac{1}{\mathcal{L}_{t}\\{\ddot{\Psi}(t)\\}(s)}=\frac{1}{s^{2}\mathcal{L}_{t}\\{\Psi(t)\\}(s)-s\Psi(0)}=\sum_{n=-2}^{\infty}a_{n}s^{n}.$
(24)

We have then as solution

$g_{0}(t)=C_{0}\left(a_{-2}t+a_{-1}+\sum_{n=0}^{\infty}a_{n}\delta^{(n)}(t)\right).$
(25)

We demand that the constant $C_{0}$ must be equal to a number where it holds
the time-reversal symmetry

$g_{0}(t)-C_{0}\sum_{n=0}^{\infty}a_{n}\delta^{(n)}(\tau-t)=1-g_{0}(\tau-t)+C_{0}\sum_{n=0}^{\infty}a_{n}\delta^{(n)}(t).$
(26)

Since the Dirac deltas and their derivatives will cancel out, the constant
will be

$C_{0}=\frac{1}{a_{-2}\tau+2a_{-1}}.$ (27)

In this manner, the optimal protocol, by construction, will be

$\displaystyle g^{*}(t)=\frac{1}{2}\left(g_{0}(t)+1-g_{0}(\tau-t)+\right.$
(28)
$\displaystyle\left.C_{0}\sum_{n=0}^{\infty}a_{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))\right),$
(29)

in which, by substituting Eq. (25) and (27), we will have

$g^{*}(t)=-\frac{t-a_{-1}}{a_{-2}\tau+2a_{-1}}-\sum_{n=0}^{\infty}\frac{a_{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))}{a_{-2}\tau+2a_{-1}},$
(30)

which is the main result of this work. As expected, such a solution shows, by
the definitions of the terms $a_{n}$ that compose it, that the optimal
protocol is very related to the relaxation function of the system.

##  IV Examples

To illustrate the consistency of the method, I find the optimal protocols of
many examples of physical systems performing an isothermal process by using
the procedure above.

###  IV.1 Overdamped Brownian motion

We consider in this example a white noise overdamped Brownian motion subjected
to a time-dependent harmonic potential, with the mass of the system equal to
one, $\gamma$ as a damping coefficient and $\omega_{0}$ as the natural
frequency of the potential. The relaxation function for both moving laser and
stiffening traps [3] are given by

$\Psi(t)=\Psi_{0}\exp{\left(-\frac{|t|}{\tau_{R}}\right)},$ (31)

where $\tau_{R}$ is the relaxation timescale of each case. Applying the
method, the terms $a_{n}$ will be

$a_{-2}=-1,\quad a_{-1}=-\tau_{R},\quad a_{n}=0,\,n\geq 0.$ (32)

Therefore

$g^{*}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}},$ (33)

which was first calculated by Schmiedl and Seifert [6] for the full dynamics,
but which is identical to the linear-response regime [3].

###  IV.2 Underdamped Brownian motion

We consider in this example a white noise overdamped Brownian motion subjected
to a time-dependent harmonic potential, with $m$ as the mass of the particle,
$\gamma$ as a damping coefficient, and $\omega_{0}$ as the natural frequency
of the potential. The relaxation function for moving laser trap [3] is given
by

$\Psi(t)=\Psi_{0}\exp{\left(-\frac{\gamma}{2}|t|\right)}\left(\cos{\omega
t}+\frac{\gamma}{2\omega}\sin{\omega|t|}\right),$ (34)

where $\omega=\sqrt{\omega_{0}^{2}-\gamma^{2}/4}$. Applying the method, the
terms $a_{n}$ will be

$a_{-2}=-1,\quad a_{-1}=-\tau_{R},\quad a_{0}=1/\omega_{0}^{2},\quad
a_{n}=0,\,n\geq 1.$ (35)

Therefore, the optimal protocol will be

$g^{*}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}}+\frac{\tau_{R}(\delta(t)-\delta(\tau-t))}{\gamma(\tau+2\tau_{R})},$
(36)

where $\tau_{R}$ is the relaxation timescale of the system. This result was
first calculated by Gomez-Marin and co-authors for the full dynamics [7], but
which is identical to the linear-response regime.

For the stiffening trap case, the relaxation function is [5]

$\Psi(t)=\Psi_{0}\exp{\left(-\gamma|t|\right)}\left[\frac{2\omega_{0}^{2}}{\omega^{2}}\right.\\\
\left.+\left(\frac{\omega^{2}-2\omega_{0}^{2}}{\omega^{2}}\right)\cos{\omega
t}+\frac{\gamma}{\omega}\sin{\omega|t|}\right],$ (37)

where $\omega=\sqrt{4\omega_{0}^{2}-\gamma^{2}}$. The coefficients $a_{n}$
will be

$a_{-2}=-1,\quad a_{-1}=-\tau_{R},\quad$ (38)
$a_{n-2}=-\frac{1}{n!}\frac{\partial^{n}}{\partial
s^{n}}\frac{(s+\gamma)(s^{2}+2\gamma
s+4\omega_{0}^{2})}{2\gamma(2\gamma+s)}\Big{|}_{s=0},$ (39)

for $n\geq 2$, where $\tau_{R}$ is the relaxation timescale of the system. The
optimal protocol will be

$g^{*}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}}+\sum_{n=0}^{\infty}\frac{a_{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))}{\tau+2\tau_{R}}.$
(40)

###  IV.3 Sinc relaxation function

In Ref. [8], when we apply the method of time average in a thermally isolated
system performing an adiabatic process, we produce a new one performing an
isothermal one with a typical relaxation time. In particular, for thermally
isolated systems that have a relaxation function equal to

$\Psi(t)=\Psi_{0}\cos{(\omega t)},$ (41)

will have for time-averaged relaxation function

$\overline{\Psi}(t)=\Psi_{0}\,\text{sinc}\left(\frac{\pi}{2}\frac{t}{\tau_{R}}\right),$
(42)

where $\tau_{R}$ is the relaxation timescale of the system. Applying the
method, the terms $a_{n}$ will be

$a_{-2}=-1,\quad a_{-1}=-\tau_{R},$ (43)
$a_{n-2}=\frac{1}{n!}\frac{\partial^{n}}{\partial
s^{n}}\frac{1}{(-1+2\tau_{R}s\,\text{acot}(2\tau_{R}s/\pi)/\pi)}\Big{|}_{s=0},\,n\geq
2.$ (44)

The optimal protocol will be

$g^{*}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}}+\sum_{n=0}^{\infty}\frac{a_{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))}{\tau+2\tau_{R}}.$
(45)

###  IV.4 Gaussian relaxation function

A relaxation function that satisfies the criteria of compatibility with the
Second Law of Thermodynamics [5] is the Gaussian relaxation function

$\Psi(t)=\Psi_{0}\exp\left(-\frac{\pi}{4}\left(\frac{t}{\tau_{R}}\right)^{2}\right),$
(46)

where $\tau_{R}$ is the relaxation timescale of the system. Applying the
method, the terms $a_{n}$ will be

$a_{-2}=-1,\quad a_{-1}=-\tau_{R},$ (47)

and

$a_{n-2}=\frac{1}{n!}\frac{\partial^{n}}{\partial
s^{n}}\frac{1}{(-1+\tau_{R}s\exp(\tau_{R}^{2}s^{2}/\pi)\text{Erfc}(\tau_{R}s/\sqrt{\pi})}\Big{|}_{s=0},$
(48)

for $n\geq 2$. Here $\text{Erfc}(x)$ is the complementary error function. The
optimal protocol will be

$g^{*}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}}+\sum_{n=0}^{\infty}\frac{a_{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))}{\tau+2\tau_{R}}.$
(49)

###  IV.5 Bessel relaxation function

The Bessel relaxation function is given by

$\Psi(t)=\Psi_{0}J_{0}\left(\frac{t}{\tau_{R}}\right),$ (50)

where $J_{0}$ is the Bessel function of the first kind with $\nu=0$ and
$\tau_{R}$ is its relaxation timescale. It satisfies the criteria for
compatibility with the Second Law of Thermodynamics. Such relaxation function
can model the Ising chain subjected to a time-dependent magnetic field and
evolving in time at equilibrium accordingly to Glauber-Ising dynamics [9].
Applying the method, the terms $a_{n}$ will be

$a_{-2}=-1,\quad a_{-1}=-\tau_{R},$ (51)

and

$a_{n-2}=-\frac{1}{n!}\frac{\partial^{n}}{\partial
s^{n}}\frac{1}{(1+\tau_{R}^{2}s^{2}+\tau_{R}s\sqrt{1+\tau_{R}^{2}s^{2}})}\Big{|}_{s=0},$
(52)

with $n\geq 2$. The optimal protocol will be

$g^{*}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}}+\sum_{n=0}^{\infty}\frac{a_{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))}{\tau+2\tau_{R}}.$
(53)

##  V Discussion

###  V.1 Continuous part

For all examples treated here, the continuous part of the optimal protocol
$g^{*}_{C}(t)$ was given by

$g^{*}_{C}(t)=\frac{t+\tau_{R}}{\tau+2\tau_{R}}.$ (54)

Such a solution is consistent with previous results for extreme cases where
$\tau\gg\tau_{R}$ and $\tau\ll\tau_{R}$ [3], where

$\lim_{\tau\gg\tau_{R}}g^{*}_{C}(t)=\frac{t}{\tau},\quad\lim_{\tau\ll\tau_{R}}g^{*}_{C}(t)=\frac{1}{2}.$
(55)

Also, this continuous part is restrained

$0\leq g^{*}_{C}(t)\leq 1,$ (56)

for all $t\in[0,\tau]$, $0\leq\tau/\tau_{R}\leq\infty$, meaning that it can be
used in linear-response theory. Therefore, the solution is physically
consistent.

###  V.2 Singular part

For all examples treated here, the singular part of the optimal protocol
$g^{*}_{C}(t)$ was given by

$g^{*}_{S}(t)=\sum_{n=0}^{\infty}\frac{b_{n}\tau_{R}^{n}(\delta^{(n)}(t)-\delta^{(n)}(\tau-t))}{\tau+2\tau_{R}}.$
(57)

where $b_{n}$ is a number independent of $\tau/\tau_{R}$. For extreme cases
where $\tau\gg\tau_{R}$ and $\tau\ll\tau_{R}$, we have

$\lim_{\tau\gg\tau_{R}}g^{*}_{S}(t)=0,\quad\lim_{\tau\ll\tau_{R}}g^{*}_{C}(t)=0.$
(58)

Another important point is the behavior of the value of the Dirac deltas and
their derivatives calculated at the extreme points $t=0$ and $t=\tau$. First,
I consider the Dirac delta defined as

$\frac{d\Theta}{dt}=\delta(t),$ (59)

where $\Theta(t)$ is the Heaviside theta function, which makes its value and
the values of its derivatives calculated at $t=0$ free of choice. In order to
satisfy the time-reversal symmetry (26), one must choose

$\delta^{(n)}(0)=\delta^{(n)}(\tau),$ (60)

for any $n\in\mathcal{N}$. Since $\tau$ is a positive number, then
$\delta^{(n)}(0)=\delta^{(n)}(\tau)=0$. In particular, it is expected too that
the integral of the optimal protocol must be time-reversal symmetric.
Therefore, one must have

$\Theta(0)=\Theta(\tau)=1.$ (61)

Also, will the Dirac deltas be outside the hypothesis of linear response
theory? To answer this question is necessary to understand how can one
implement such delta peaks in practice in the laboratory. In my opinion,
operating protocols considering “peaks” to infinity and considering “peaks” in
the derivatives does not make any sense. How can one guarantee the optimality
of the irreversible work then? I suggest just including the constant term
calculated with the singular part in the definition of the irreversible work,
whose primary definition was made considering only jumps at the beginning and
final point of the process. The new optimal irreversible work will become

$\displaystyle W_{\rm irr}^{*}$ $\displaystyle=\
\frac{\delta\lambda^{2}}{2}\Psi(0)+\frac{\delta\lambda^{2}}{2}\int_{0}^{\tau}\dot{\Psi}_{0}(\tau-t)g^{*}(t)dt$
(62)
$\displaystyle+\frac{\delta\lambda^{2}}{2}\sum_{n=0}^{\infty}\frac{a_{n}(-1)^{n}[\dot{\Psi}^{(n)}(\tau)+\dot{\Psi}^{(n)}(0)]}{\tau+2\tau_{R}},$
(63)

where the continuous part of the optimal protocol is now the unique part of
the optimal protocol. So, in practice, the Dirac deltas are nothing more than
artifices that the method of solving the integral equation presents to taking
into account the terms that have been not defined in the treatment of the
irreversible work, where the constants terms that appear from the Dirac deltas
and its derivatives were not considered. In this manner, the Dirac deltas and
derivatives are in accord to the hypothesis of linear response theory.

##  VI Final remarks

In this work, I presented a method to solve the Euler-Lagrange integral
equation that furnishes the optimal protocol for the irreversible work and its
variance. It relies only on the time-reversal symmetric property of the
optimal protocol. Many examples are solved to illustrate the consistency of
the method. The solution is proven to be physically consistent. To overcome
the hypothesis consistency problem with the linear response theory of the
singular part, an interpretation of the appearance of delta peaks is
presented, where the constant terms of such singular solutions should be
included at irreversible work at the beginning of the optimization problem.

## References

  * Polyanin and Manzhirov [2008] A. D. Polyanin and A. V. Manzhirov,  _Handbook of integral equations_ (CRC press, 2008).
  * Nazé [2023] P. Nazé, arXiv preprint arXiv:2304.11965 (2023).
  * Nazé  _et al._ [2022] P. Nazé, S. Deffner, and M. V. Bonança, Journal of Physics Communications 6, 083001 (2022).
  * Kubo  _et al._ [2012] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical physics II: nonequilibrium statistical mechanics_, Vol. 31 (Springer Science & Business Media, 2012).
  * Nazé and Bonança [2020] P. Nazé and M. V. S. Bonança, Journal of Statistical Mechanics: Theory and Experiment 2020, 013206 (2020).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Physical review letters 98, 108301 (2007).
  * Gomez-Marin  _et al._ [2008] A. Gomez-Marin, T. Schmiedl, and U. Seifert, The Journal of chemical physics 129, 024114 (2008).
  * Nazé [2022] P. Nazé, arXiv preprint arXiv:2210.16116 (2022).
  * Glauber [1963] R. J. Glauber, Journal of mathematical physics 4, 294 (1963).