Supplementary_materials

Camille Nevoret

# Debiasing the estimate of treatment effect on the treated with time-varying
counfounders

Camille Nevoret1,2,31,2,3affiliationmark:  Sandrine
Katsahian3,43,4affiliationmark:  and Agathe Guilloux11affiliationmark:
11affiliationmark:  Université Paris-Saclay, CNRS, Univ Evry, Laboratoire de
Mathématiques et Modélisation d’Evry, 91037, Evry-Courcouronnes, France.
22affiliationmark:  CRC, Université de Paris, INRIA, EPI HEKA, INSERM UMR
1138, Sorbonne Université, France.
33affiliationmark:  CEMKA, 43 Bd Maréchal Joffre, 92340, Bourg La Reine,
France.
44affiliationmark:  Inserm, Centre d’Investigation Clinique 1418 (CIC1418)
Epidémiologie Clinique, Paris, France.

###### Abstract

With the increased availability of large health databases comes the
opportunity of evaluating treatment effect on new data sources. Through these
databases time-dependent outcomes can be analysed as events that can be
measured using counting processes. Estimating average treatment effect on the
treated (ATT) requires modelling of time-varying covariate and time-dependent
treatment and outcome. Gran et al. proposed an easy-to-implement method based
on additive intensity regression to estimate ATT. We introduce a debiased
estimate of the ATT based on a generalization of the Gran’s model for a
potentially repeated outcome and in the presence of multiple time-dependent
covariates and baseline covariates. Simulation analyses show that our
corrected estimator outperforms Gran’s uncorrected estimator. Our method is
applied to intensive care real-life data from MIMIC-III databases to estimate
vasoppressors effect on patients with sepsis.

##  1 Introduction

Health data warehouse are increasingly available and could be an alternative
to the prospective cohort studies traditionally designed for epidemiological
research because data are immediately available and projects are thus less
time consuming and expensive to run. In most data sources such as health
registries or electronic health records, patients are followed over time,
resulting in longitudinal data. Statistical analyses are therefore performed
in the presence of time-dependent covariates and time-varying treatment. The
development of longitudinal databases makes the study of recurring events
easier such as repeated events or events that can be measured using counting
processes. For example, they are useful when studying rehospitalizations,
relapses in cancer, asthma attacks, or multiple sclerosis relapses.

Time-dependent confounding occurs when a covariate is both affected by past
treatment and associated with future treatment choice and outcome. Many
longitudinal studies aim to estimate the causal effect of time-varying
exposure on the outcome. The main issue of these analyses is to handle such
confounding. Specific methods address this issue by taking into account time-
confounding and thus avoid biased estimates [1, 2, 3]. The most common method
for dealing with time-dependent confunding is the inverse probability weighted
marginal structural Cox model [4]. Instrumental variable methods [5] or
methods derived from g-methods can also be used.

Two types of treatment effect estimator, the average total treatment effect
(ATE) and the average treatment effect on the treated (ATT), can be
identified. The most commonly used, the ATE, aims identifying the treatment
effect in the whole population, eligible for treatment. It allows to answer
the question "What would happen if all the patients had received the evaluated
treatment, at every timepoint as compared with they had received its
comparator?". The ATT aims at identifying the average treatment effect only in
patients who actually received the treatment of interest. With this method,
the effect of treatment is not assessed in patients who never received the
treatment. The choice of estimator depends on causal research question of
interest and on its availability for the method used. For example, it may be
more interesting to evaluate the treatment effect on the population that
received that treatment rather than on the whole population. Indeed, in
clinical practice, the prescription of a treatment is motivated by a set of
indicators, as patients clinical characteristics or the likelihood that the
treatment will be effective for that patient. In addition, depending on the
data available, only information on treated patients may be available.

The ATT can be estimated mainly by g-estimation of structural nested models
[6]. However, this could be complex and prevents from its common use in
practice [7]. Gran et al. proposed an easy-to-implement method that allows to
estimate the treatment effect in treated patients and to study intermediate
time-dependent covariate trajectories [8]. This model was presented for end-
point outcome and an unique time-varying covariate. It consists of two steps.
The first step is to model time-varying covariates on the time interval where
patients have not yet been treated. Then, the causal treatment effect is
estimated by modelling the untreated evolution of time-varying covariates.
Finally, treatment effect is estimated using intensity regression models. The
proposed method is developed for additive intensity regression models [9].
When studying causality in censored data, this method allows a more explicit
approach to processes [10] and explicit derivations that are not available for
the Cox model [11, 12, 13, 14]. This model seems to be a good alternative to
the Cox model [15, 16].

We propose in this paper to generalize the Gran’s model for a potentially
repeated outcome and in the presence of multiple time-dependent covariates and
baseline covariates. Moreover, we propose a debiased estimate of the ATT based
on the error produced during the modelling of the covariates where patients
are untreated. Simulation analyses show that the corrected estimator that we
propose outperforms Gran’s uncorrected estimator.

The paper is organized as follows. Section section:model_notations introduces
the different notations and assumptions used in the rest of the paper. Section
treatment_estimation presents the additive intensity regression model and the
proposed correction for a general time-dependent covariates modelling method.
Section Section_VAR specifies the correction in case of vector autoregressive
model. Section simulation summarises results from simulations. Section
Application presents the results of the models on intensive care real-life
data from MIMIC-III data base [17].

##  2 Model and notations

###  2.1 Model

Consider an individual who suffers from a specific disease with a follow-up
over time. At time $S$, he/she initiates his/her treatment. We assume that
once started the treatment is continued, so that the treatment process $D$ is
null before $S$ and equals $1$ after $S$. Together with the treatment, we
observe $d_{Z}$-dimensional baseline covariates
$Z=\big{(}Z^{1},\ldots,Z^{d_{Z}}\big{)}$ and time-dependent
$d_{X}$-dimensional covariates $t\mapsto X(t)=(X^{1}(t),\ldots,X^{d_{X}}(t))$.
We focus on pathologies for which the outcome of interest can be measured via
a counting process $N$ ($N(t)$ being the number of events recorded in
$[0,t]$), see examples in Intro and Application sections. We assume also that
$N$ verifies the Aalen additive intensity model [9].

######  Assumption 1.

With respect $\big{(}\mathcal{F}^{\star}_{t}\big{)}_{t\geq 0}$, the historical
filtration spanned by $Z$, $X$, $D$, and $N$, we assume the process $X$ is
predictable and that $N$ has the following intensity:

$\mu^{*}(t)=\alpha_{0}^{*}(t)+\alpha_{Z}^{*}(t)Z+\alpha_{X}^{*}(t)X(t)+\alpha_{D}^{*}(t)D(t).$

As in [8], we will assume further for simplicity in what follows that
independent censoring can occur. Censoring can be assumed to be independent
when analyzing data from intensive care unit or administrative databases but
not in all pathologies database. Methods allowing dependent censoring such as
adjustment using inverse probability of censoring weighting should be used
[18].

###  2.2 Counterfactual quantities, assumptions

We now consider that we can intervene on the treatment decision at time $S$.
Given the definitions introduced above, we have with probability $1$

$\displaystyle D(S)=1,\;D(s)=0\text{ for }s<S\text{ and }\;D(t)=1\text{ for
}t\geq S.$

However assuming that the positivity assumption, see e.g. [19] for example, is
fulfilled implies that

$0<\mathbb{P}\big{(}D(t)=1|D(t-)=0,\mathcal{F}_{t-}\big{)}<1.$

At any time $t$ before treatment initiation, it can be started or not with a
positive probability. This excludes certain pathologies. In practice the
counterfactual reasoning requires this positiving assumption which has to be
checked on the data.

In the hypothetical situation where the treatment is not initiated at $S$, we
denote by $X^{0|S}$ the covariate process after $S$. This counterfactual
process corresponding to what would have been observed after $S$. It is
important to remark that it cannot be observed, so that its values will have
to be estimated (see Section Section_VAR). The potential covariate process in
the situation where the treatment is initiated at $S$ is denoted by $X^{1|S}$.

We also introduce the two potential counting processes $N^{1|S}$ recording the
event history after $S$ when the treatment is actually initiated at $S$ and
$N^{0|S}$ when the treatment is not initiated. Their respective intensities
are denoted by $\mu^{1|S}$ and $\mu^{0|S}$.

The consistency assumption in this case, see [19], implies that for all $t\geq
S=s$

$X^{1|S=s}(t)=X(t)\text{ and }N(t)=N^{1|S=s}(t).$

We also assume conditional exchangeability conditionally to $S=s$, for all
$a=0,1$, all $t\geq s$ and conditionally to the past $\mathcal{F}_{s-}$:

$\displaystyle\mathbb{E}(dN^{a|S=s}(t)|D(s),\mathcal{F}_{s-})=\mathbb{E}(dN^{a|S=s}(t)|\mathcal{F}_{s-})\text{
and }$
$\displaystyle\mathbb{E}(X^{a|S=s}(t)|D(s),\mathcal{F}_{s-})=\mathbb{E}(X^{a|S=s}(t)|\mathcal{F}_{s-}).$

This means that at time $S=s$ conditionally to the past $\mathcal{F}_{s-}$, if
two perfectly similar people (with exactly the same history of treatment and
exposure) one receiving the treatment, the other one not receiving the
treatment, have a counterfactual outcome which is the other one actual
observed outcome. Under Assumption (1), this can be rewritten, for $t\geq
S=s$, as

$\displaystyle\mu^{1|S=s}(t)$
$\displaystyle=\alpha_{0}^{*}(t)+\alpha_{D}^{*}(t)+\alpha_{Z}^{*}Z+\alpha_{X}^{*}(t)X^{1|S=s}(t)\text{
and }$ $\displaystyle\mu^{0|S=s}(t)$
$\displaystyle=\alpha_{0}^{*}(t)+\alpha_{Z}^{*}Z+\alpha_{X}^{*}(t)X^{0|S=s}(t).$

The exchangeability assumption can hold only if no unmeasured covariate
influences our model. This is a strong assumption in practice but is common in
causality reasoning, see e.g. [8, 19] among many others. Following [19], the
clinician opinion expertise is key in proposing a list of covariates
associated with treatment indication to be included in the baseline covariates
$Z$ and the time dependant covariates process $X$ so that this assumption is
fulfilled at least approximately conditional on all variables recorded.

A stronger assumption can be found in the literature [19, 16]: conditionally
to $S=s$, for all $a=0,1$ and conditionally to the past $\mathcal{F}_{s-}$

$X^{a|S}(t)\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}D(s)\text{ and
}N^{a|S}(t)\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}D(s),\text{ for all }t\geq s.$

###  2.3 ATT definition and causal estimate

We assume an additive intensity regression model for the outcome processes.
Under the model and the assumptions presented in previous paragraph, the time
varying causal intensity difference can be then defined as

$\displaystyle d^{*}(t,s)$ $\displaystyle=\mu^{1|S=s}(t)-\mu^{0|S=s}(t)$
$\displaystyle=\alpha_{D}^{*}(t)D(t)+\alpha_{X}^{*}(t)(X^{1|S=s}(t)-X^{0|S=s}(t)).$

With these definitions, we can rewrite the intensity of $N$ as

$\mu^{*}(t)=\alpha^{*}_{0}(t)+\alpha^{*}_{Z}Z+\alpha^{*}_{X}(t)X^{0|S=s}(t)+d^{*}_{i}(t,s)D(t),$

and, for $t\geq S$, the time varying ATT as

$\text{ATT}(t)=\mathbb{E}[d^{*}(t,S)|t\geq S]=\alpha_{D}^{*}(t)+\\\
\alpha_{X}^{*}(t)\mathbb{E}\big{[}X^{1|S}(t)-X^{0|S}(t)\big{|}t\geq S].$

This has been described in [8] and strongly relies on the Aalen additive
intensity model of Assumption (1). We refer the reader to [15] for more
discussions on that matter.

##  3 Additive intensity model for the treatment effect estimation

Our data consists in the observation for $n$ independent individuals
$i=1,\ldots n$ of $Z_{i}$, $X_{i}$, $S_{i}$, $D_{i}$ and $N_{i}$. Each
individual is followed between $t=0$ and $t=\tau_{i}$ which can be the end
time of the study. The intensity of $N_{i}$ can be written as:

$\displaystyle\mu^{*}_{i}(t)$
$\displaystyle=\alpha^{*}_{0}(t)+\alpha^{*}_{Z}Z_{i}+\alpha^{*}_{X}(t)X_{i}^{0|S=S_{i}}(t)+d^{*}_{i}(t,S_{i})D_{i}(t)$
$\displaystyle=W_{i}(t)^{\top}A^{*}(t).$

We can recognize the writing of standard regression where treatment effect is
adjusting on treatment, the baseline covariates and the counterfactual
covariates, $X^{0|S=S_{i}}_{i}(t)$. We remark that in that case the ATT is the
treatment regression coefficient.

###  3.1 Identification of the treatment effect estimate

The aim of these section is to define an estimator of $A^{*}$ in the
hypothetical situation, false in practice, where the trajectories of the
counterfactual covariate processes $X_{i}^{0|S=S_{i}}$ for $i=1,\ldots,n$ and
$t\geq S_{i}$, or equivalently the processes $W_{i}$, would be observed. We
simply show that the traditional empirical least-squares risk (see e.g. [20])
for the additive intensity model would lead to a correct estimator of $A^{*}$.
Let

$A(t)=(\alpha_{0}(t),\alpha_{Z},\alpha_{X}(t),d_{i}(t,S_{i}))^{\top}\text{,
}\forall t\geq 0,$

be a candidate for the estimation of $A^{\star}$. The squared risk of
$A\in\mathbb{R}^{d_{X}+d_{Z}+2}$ is defined as, see [20]:

$\displaystyle r_{n}(A)$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}A(t)^{\top}W_{i}(t)W_{i}(t)^{\top}A(t)dt$
(1)
$\displaystyle-\frac{2}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}A(t)^{\top}W_{i}(t)dN_{i}(t).$

The fact that

$\widehat{A}=\text{arg}\min\limits_{A}r_{n}(A)$ (2)

is a good candidate for the estimation of $A^{*}$ can be checked by computing
the expectation of $r_{n}(A)$. Calculations presented in Supplementary
Materials, section 1.1, allow us to deduce that $A^{*}$ can be expressed as:

$A^{*}=\text{arg}\min\limits_{A}\mathbb{E}[r_{n}(A)],$

which shows that $\widehat{A}$ is an estimator of $A^{*}$.

###  3.2 Estimation of counterfactual covariates

In our situation, the counterfactual trajectories $X_{i}^{0|S=S_{i}}$ cannot
be observed, as in the data, patient $i$ is treated as time $S_{i}$. To
estimate $A^{\star}$, we need to estimated these trajectories. We give an
exemple of such estimations in Section Section_VAR. We consider for now that
we can model the trajectories and estimate the values of $X^{0|S=S_{i}}(t)$
for $t\geq S_{i}$ by $\widetilde{X}^{0|S=S_{i}}_{i}(t)$. Being an estimation,
$\widetilde{X}^{0|S=S_{i}}_{i}(t)$ is not equal to $X^{0|S=S_{i}}(t)$ but we
assume that we can write:

$\displaystyle\widetilde{X}^{0|S=S_{i}}_{i}(t)$
$\displaystyle=X^{0|S=S_{i}}_{i}(t)+\epsilon_{i}(t),\text{ for all }t\geq
S_{i},$ (3)

where $\epsilon_{i}(t)$ is centered and has finite variance. We finally denote
by $\widetilde{W}_{i}(t)$ the vector $W_{i}(t)+\xi_{i}(t)$ with
$\xi_{i}(t)=(0,0_{\mathbb{R}^{d_{z}}},\epsilon_{i}(t),0)$ for $t\geq S_{i}$
and $\xi_{i}(t)=(0,\ldots,0)$ for $t<S_{i}$.

Mimicking the computation in the above paragraph, let us introduce
$\widetilde{R_{n}}(A)$ the squared risk of $A$ using counterfactual covariates
estimate $\widetilde{W}_{i}(t)$:

$\displaystyle\widetilde{R_{n}}(A)=\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}A(t)^{\top}\widetilde{W}_{i}(t)\widetilde{W}_{i}(t)^{\top}A(t)dt$
$\displaystyle-\frac{2}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}A(t)^{\top}\widetilde{W}_{i}(t)dN_{i}(t).$

Computations detailed in Equation (1) of Supplementary materials allow to
write $\widetilde{R_{n}}(A)$ in function of $r_{n}(A)$ using Equation (3) as

$\displaystyle r_{n}(A)$ $\displaystyle=\widetilde{R}_{n}(A)$ (4)
$\displaystyle+\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}\alpha_{X}(t)^{\top}\epsilon_{i}(t)\epsilon_{i}(t)^{\top}\alpha_{X}(t)dt.$

Using Equation (4), the expectation of $\widetilde{R_{n}}$ is given by the
following expression:

$\mathbb{E}[\widetilde{R_{n}}(A)]=\mathbb{E}[r_{n}(A)]\\\
-\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}\alpha_{X}(t)^{\top}\mathbb{E}[\epsilon_{i}(t)\epsilon_{i}(t)^{\top}]\alpha_{X}(t)dt.$
(5)

Hence $\text{arg}\min\limits_{A}\mathbb{E}[\widetilde{R}_{n}(A)]\neq A^{*}$.

In that sens, the proposal of Gran et al. [8] is biased. We will introduce in
the following section a novel estimator of $A^{*}$ (and the ATT) taking into
account the counterfactuals modelling error .

###  3.3 Debiased treatment effect estimate

We showed that minimizing $\widehat{R_{n}}$ would had to a biased estimated.
We propose in this paragraph to correct this bias. Equation (5) implies that

$\widehat{R_{n}}(A)=\widetilde{R_{n}}(A)+\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}\alpha_{X}(t)^{\top}\widehat{\mathbb{E}[\epsilon_{i}(t)\epsilon_{i}(t)^{\top}]}\alpha_{X}(t)dt$

is the correct loss to estimated $A^{\star}$. However, as the unknown
expectations $\mathbb{E}[\epsilon_{i}(t)\epsilon_{i}(t)^{\top}]$ appear, it
cannot be directly used. For some models for the $X_{i}^{0|S=S_{i}}$
trajectories, these expectations can however be estimated. This is the case
with the vector autoregressive model presented in the next section. Now
considering that we have access to such an estimation
$\widehat{\mathbb{E}[\epsilon_{i}(t)\epsilon_{i}(t)^{\top}]}$ for all $i$ and
$t\geq S_{i}$, we can defined our loss as:

$\widehat{R_{n}}(A)=\widetilde{R_{n}}(A)+\widehat{\text{bias}}(A),$

with

$\begin{split}\widehat{\text{bias}}(A)=\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}\alpha_{X}(t)^{\top}\widehat{\mathbb{E}[\epsilon_{i}(t)\epsilon_{i}(t)^{\top}]}\alpha_{X}(t)dt.\end{split}$
(6)

Finally, the debiased estimator of $A^{*}$ is given by:

$\displaystyle\widehat{A}_{\text{debiased}}$
$\displaystyle=\text{arg}\min\limits_{A}\widehat{R_{n}}(A)$
$\displaystyle=\text{arg}\min\limits_{A}\frac{1}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}A(t)\left[\widetilde{W}_{i}(t)\widetilde{W}_{i}(t)^{\top}\right.$
$\displaystyle\quad-\left.\begin{pmatrix}0&\cdots&0\\\
\vdots&\widehat{\mathbb{E}[\epsilon_{i}(t)\epsilon_{i}(t)^{\top}]}&\vdots\\\
0&\cdots&0\end{pmatrix}\right]A(t)dt$
$\displaystyle-\frac{2}{n}\sum^{n}_{i=1}\int_{0}^{\tau_{i}}A(t)\widetilde{W}_{i}(t)dN_{i}(t).$
(7)

We now propose a simple model for the continuous time-varying
$d_{x}$-dimensional covariate. $X_{i}$ in which the estimation
$\widehat{\text{bias}}(A)$ is easily computable.

##  4 Example with a vector autoregressive model VAR(1)

###  4.1 Explicit writing of the bias

In the following, the observation period from $0$ to $\max\limits_{i}\tau_{i}$
is split into observed time intervals defined by following times
$t_{0}=0,\ldots,t_{k},\ldots t_{K}=\max\limits_{i}\tau_{i}$. We assume that
the time varying covariates are constant over these intervals such that for
all $i$, $k$ and $t\in[t_{k},t_{k+1})$, $X_{i}(t)=X_{i}(t_{k})$, and by
extension $W_{i}(t)=W_{i}(t_{k})$.

Considering the case of continuous time-varying covariates, the simplest model
is the vector auto-regressive model [21]. We explicit the computation in the
case of the VAR(1) model:

$X_{i}(t_{k})=\beta_{0}^{i\star}+\Pi^{\star}X_{i}(t_{k-1})+\omega_{i}(t_{k}),$

for all $i$ and $t$, where $\Pi^{\star}$ is the coefficient matrix and
$\omega_{i}$ is an unobservable zero mean white noise vector process with time
invariant covariance matrix $\Sigma^{\star}$.

Following the literature on VAR model, the prediction of $X(t_{k+l})$ using
the past up $t_{k}$ is given iteratively by:

$\displaystyle\widetilde{X}_{i}(t_{k+l})$
$\displaystyle=\beta_{0}^{i\star}+\Pi^{\star}\widetilde{X}_{i}(t_{k+l-1})$
$\displaystyle=\sum_{j=0}^{l-1}{\Pi^{\star}}^{j}\beta_{0}^{i\star}+{\Pi^{\star}}^{l}X_{i}(t_{k}).$

The modelling error can be expressed as follow:

$X_{i}(t_{k+l})-\widetilde{X}_{i}(t_{k+l})=\sum_{j=0}^{l-1}{\Pi^{\star}}^{j}\omega_{i}(t_{k+l-j}),$
(8)

with $\Pi^{0}=\mathbf{I}_{d_{X}}$ the identity matrix of size $d_{X}$. The
mean squared error (MSE) matrix for $\widetilde{X}_{i}(t_{k+l})$ is written
as:

$\displaystyle\Sigma(l)$
$\displaystyle=\text{MSE}(X_{i}(t_{k+l})-\widetilde{X}_{i}(t_{k+l}))=\sum_{j=0}^{l-1}{\Pi^{\star}}^{j}\Sigma{\Pi^{\star}}^{\top
j}.$ (9)

The parameters of the VAR(1) process are estimated using multivariate least
squares [22]. Estimated parameters is noted $\widehat{\Pi}$ and
$\widehat{\beta_{0}^{i}}$. The best linear predictor of $X_{i}(t_{k+l})$
becomes:

$\widehat{\widetilde{X}}_{i}(t_{k+l})=\widehat{\beta_{0}^{i}}+\widehat{\Pi}\widehat{\widetilde{X}}_{i}(t_{k+l-1}).$
(10)

Finally, the estimated MSE matrix of the 1-step forecast is computing using
equation (9) as:

$\widehat{\Sigma}(l)=\sum_{j=0}^{l-1}\widehat{\Pi}^{j}\widehat{\Sigma}\widehat{\Pi}^{\top^{j}}.$

###  4.2 Algorithm

Our algorithm combines two step.

####  4.2.1 First step: counterfactual estimates

  1. 1.

After filtering, only observations where patients are not treated are kept.

  2. 2.

For each time-varying covariates $j$, the matrix in which $t_{\max}^{i}$ is
the last time without treatment for the patient $i$ is designed, i.0.e the
discret observe time before $S_{i}$.

  3. 3.

The linear regression with the design matrix of Equation (11) and the vector
of predicted outputs
$(X_{1}^{j}(1),\cdots,X_{1}^{j}(t_{\max}^{1}),\cdots,X_{n}^{j}(1),\cdots,X_{n}^{j}(t_{\max}^{n}))$
is computed.

  4. 4.

Thus we get the estimates of $(\Sigma_{jj})$, $(\Pi_{j,q})_{1\leq q\leq
d_{X}+d_{Z}}$ and $\beta_{0}^{i}$.

  5. 5.

Counterfactuals are modelled following Equation (10) using estimates of $\Pi$
and $\beta_{0}^{i}$, and are denoted $\widehat{\widetilde{X_{i}}}(t_{k})$ for
$i=1,\ldots,n$ and $t_{k}\geq t_{i}$.

$\begin{pmatrix}1&Z_{1}^{1}&\cdots&Z_{1}^{d_{Z}}&{\scriptstyle
X_{1}^{1}(0)}&\cdots&{\scriptstyle X_{1}^{d_{X}}(0)}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&Z_{1}^{1}&\cdots&Z_{1}^{d_{Z}}&{\scriptstyle
X_{1}^{1}(t_{\max}^{1}-1)}&\cdots&{\scriptstyle
X_{1}^{d_{X}}(t_{\max}^{1}-1)}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&Z_{n}^{1}&\cdots&Z_{n}^{d_{Z}}&{\scriptstyle
X_{n}^{1}(0)}&\cdots&{\scriptstyle X_{n}^{d_{X}}(0)}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&Z_{n}^{1}&\cdots&Z_{n}^{d_{Z}}&{\scriptstyle
X_{n}^{1}(t_{\max}^{n}-1)}&\cdots&{\scriptstyle
X_{n}^{d_{X}}(t_{\max}^{n}-1)}\\\ \end{pmatrix}$ (11)

####  4.2.2 Second step: ATT estimate with modelled counterfactuals

The ATT estimate is implemented in different analysis softwares such as the
aalen function of the timereg package [23, 24]. We implemented our own method
to apply our correction.

  * •

The ATT is estimating by solving Equation (2) and (3.3).

##  5 Simulation study

In this section, ATT estimations with and without correction are compared
using simulations in the VAR(1) model. This simulation study is based on the
simulation setup proposed in [8] . Data are simulated to mimic a cohort
analysis where individuals are under risk of having an event of interest which
can be prevented or delayed by treatment. Several covariates are simulated
over time. Their values depend on treatment status, and reciprocally,
treatment initiation depends on these values, see Figure 1. The probability of
having the event is then affected by both the treatment status and the
covariate values.

Figure 1: Directed Acyclic Graph

###  5.1 Data generation

Simulation setup is motivated in [8] by an application on the effect of
antiretroviral therapy (HAART) on the time to AIDS onset or death. We have
extended these simulations with several time-varying covariates and baseline
covariates.

The data are simulated on time intervals $[t_{k},t_{k+1})$ with pre-specified
time points $t_{k}=0,1,....11$. On each time interval, we define event status
$N(t_{k})$, treatment status $D(t_{k})$, $d_{X}$ time-varying covariates,
$X(t_{k})=(X^{1}(t_{k}),\dots,X^{d}(t_{k}))$ and $d_{Z}$ baseline covariates,
$Z$. Dataset with one, three and six time-varying covariates are simulated.

At initiation, all patient are untreated. Initial time-varying covariates are
simulated as uniform distribution. Three baseline covariates are simulated to
mimic age as numbers uniformly distributed, gender as a binomial distribution
and a comorbidity score as a Poisson distribution.

On the interval where the treatment has not yet been initiated, each
covariates $X^{j}(t_{k})$ will decrease steadily according to the following
equation
$X(t_{k+1})=\kappa_{D0}X(t_{k})+\mathcal{N}(0_{\mathbb{R}^{d}},\Sigma)$ with
$\kappa_{D0}$ a $d_{X}\times d_{X}$ matrix and
$\Sigma=\\{\sigma_{jj}\\}_{1\leq j\leq k}$. While from treatment initiations,
$X^{j}(t_{k})$ increase over time following the equation
$X(t_{k+1})=X_{i}(t+1)=X_{i}(t)+\kappa_{D1}(\sqrt{1000}-X_{i}(t))+\mathcal{N}(0_{\mathbb{R}^{d}},\Sigma)$
with $\kappa_{0}$ a $d_{x}$-dimensional vector and $\kappa_{D1}$ a
$d_{X}\times d_{X}$ matrix. Treatment initiation $D(t_{k})$ is simulated
according to a Bernoulli distribution with parameter
$m\sum_{j=1}^{d}\lambda_{j}e^{\sum_{j=1}^{d}\lambda_{j}X^{j}(t_{k})})$. We
assume that once treatment is initiated, the patient remains on treatment
until the end of the follow-up.

Finally, treatment status and time-varying covariates act on the number of
experienced events, $N(]t_{k},t_{k+1}])$ experienced on $]t_{k},t_{k+1}]$ .
The event times are simulated according to a non-homogeneous Poisson process
with intensity function $\delta
D(t_{k})+\delta_{0}+\sum_{j=1}^{d_{Z}}\delta_{Zj}Z^{j}+\sum_{j=1}^{d}\delta_{Xj}X^{j}(t_{k})$.
This process is generated using acceptance-rejection method called thinning
[25]. All patients are followed up until time $t_{k}=11$. The relation between
covariates, treatment and outcome is summarized by the direct acyclic graph of
Figure 1.

At the same time, the "true" counterfactual values are simulated by following
the evolution of the covariates when $D(t_{k})=0$. The simulation equation is
as follows:
$\widetilde{X}(t_{k+1})=\kappa_{D0}\widetilde{X}(t_{k})+\mathcal{N}(0_{\mathbb{R}^{d}},\Sigma)$
whatever the treatment status. The simulation algorithm is presented in
Supplementary Materials, see Section 2.

The number of time-varying covariates varied in the simulated framework. The
following results are based on simulations with 1 time-varying covariate
affecting $D$ and $N$, and simulations with three and six covariates with
$d/3$ covariates affecting $D$ and $N$, $d/3$ covariates affecting only $D$
and $d/3$ covariates affecting only $N$. Finally, the parameter $\Sigma_{jj}$
was set to the values $0.4$, $0.8$, $1.2$ and $1.6$. The values of parameters
used in the simulations are given in the appendix, see Section
LABEL:simulation_algorithme. We ran the simulations 100 times with the same
parameters on data sets of 1000 individuals. All simulations were done in R
version 4.1. Codes used for the simulations are available on gitlab
https://gitlab.com/camille.nevoret/att-estimation.

###  5.2 Analysis

As presented in the section Algorithm, the first step in the simulation study
is the modelling of counterfactuals. Within the VAR(1) model, we estimated
$\Sigma$ and $\kappa_{D0}$, used to calculate the correction of the ATT
estimator.

The data set of the 100 simulation repetitions with the true counterfactual
values is used to evaluate the "true" ATT value by Monte Carlo. At each time
$t$, we note $d^{MC}(t)$ the "true" ATT. All ATT estimates for each of the 100
simulations are compared to this value. We use the mean integrated square
error (MISE) to calculated the error between ATT estimate and the true ATT.
For a given estimator $\hat{d^{l}}(t)$ for the l-th simulation, the MISE is
given by:

$\text{MISE}^{l}=\mathbb{E}\left[\int^{\tau=11}_{0}(\widehat{d}^{l}(t)-d^{MC}(t))^{2}dt\right].$

For each of the 100 simulations, we estimate ATT using five different
estimators.

  * •

The first estimator is calculated using the aalen function of the R timereg
package. Counterfactual data use for this estimate are the one simulated as
"real counterfactuals". Thus, there are no errors related to the modelling of
counterfactuals introduced in the ATT estimation. In this sense, the other
estimators, based on modelled counterfactuals cannot make more accurate
estimations than this one. It is thus uses as a reference.

  * •

Then, two uncorrected estimators are presented. The first one corresponds to
the one presented in [8] using aalen function and modelled counterfactuals.
The second one is obtain using algorithm presented section Algorithm.

  * •

Finally, two corrected estimators are presented. They use both the algorithm
presented section Algorithm and modelled counterfactuals. The first one use
the true time invariant covariance matrix, $\Sigma$, and the true
$\kappa_{D0}$. In real life data, they are estimated when modelling the
counterfactuals. This is done in the second estimator.

Wilcoxon tests for the difference in distribution are performed between the
corrected estimator with estimated parameters and each uncorrected estimators.

###  5.3 Simulation results

Table 1 shows the mean integrated square error for three ATT estimates: the
two uncorrected estimators and the corrected estimator using estimate
parameters $\widehat{\Sigma}$ and $\widehat{\kappa}_{D0}$. We observe that,
whatever the number of covariates, the true mean integrated square error
increases with $\Sigma$. Data in bold are the lowest MISEs. The symbol ⋆
represents a significative difference given by the Wilcoxon test.

For all simulation parameters, ATT estimates using the aalen function [23, 24]
and using our recoded algorithm are very close. Our corrected estimator is
always better than uncorrected estimated. We notice significant differences
between these estimators, for all simulation scenarios with one covariate.

Table 1: Mean MISE $\pm$ standard deviation for corrected and uncorrected estimators (⋆: Significant Wilcoxon test between corrected estimator with estimated parameters and other estimators, boldface character: the smaller MISE of the corrected and uncorrected estimators) | Corrected | Uncorrected estimator
---|---|---
| estimator |  With aalen | Recoded
|  with $\widehat{\Sigma}$ | function | estimate
1 Covariate
$\Sigma_{k,k}=0.4$ | $0.034\pm 0.038$ | $0.045^{\star}\pm 0.047$ | $0.044^{\star}\pm 0.047$
$\Sigma_{k,k}=0.8$ | $0.136\pm 0.131$ | $0.178^{\star}\pm 0.155$ | $0.177^{\star}\pm 0.155$
$\Sigma_{k,k}=1.2$ | $0.412\pm 0.384$ | $0.533^{\star}\pm 0.445$ | $0.527^{\star}\pm 0.443$
$\Sigma_{k,k}=1.6$ | $0.787\pm 0.730$ | $1.026^{\star}\pm 0.853$ | $1.014^{\star}\pm 0.848$
3 Covariates
$\Sigma_{k,k}=0.4$ | $0.040\pm 0.046$ | $0.045\pm 0.051$ | $0.045\pm 0.051$
$\Sigma_{k,k}=0.8$ | $0.224\pm 0.237$ | $0.257\pm 0.265$ | $0.256\pm 0.264$
$\Sigma_{k,k}=1.2$ | $0.518\pm 0.548$ | $0.608\pm 0.604$ | $0.604\pm 0.603$
$\Sigma_{k,k}=1.6$ | $1.381\pm 1.299$ | $1.598\pm 1.434$ | $1.587\pm 1.429$
6 Covariates
$\Sigma_{k,k}=0.4$ | $0.075\pm 0.080$ | $0.081\pm 0.095$ | $0.081\pm 0.096$
$\Sigma_{k,k}=0.8$ | $0.320\pm 0.360$ | $0.338\pm 0.376$ | $0.337\pm 0.377$
$\Sigma_{k,k}=1.2$ | $0.901\pm 1.012$ | $0.949\pm 1.043$ | $0.947\pm 1.045$
$\Sigma_{k,k}=1.6$ | $2.448\pm 2.316$ | $2.567\pm 2.404$ | $2.560\pm 2.406$

Table 2 allows to compare our corrected ATT estimate with estimated
parameters, which is the one that can be used on real data, to corrected the
ATT estimate with known parameters and the estimate using the "real
counterfactuals". It appears that MISE of corrected estimates are very close
which shows that counterfactual modelling lead to good estimates of $\Sigma$
and $\kappa_{D0}$. We can also note large standard deviations of the order of
the means of the MISE. These large standard deviations do not seem to be
inherent to our estimation method. Indeed, we observe the same type of
deviation for the estimate using "real counterfactuals".

Finally, for all simulation parameters tested, the corrected estimators appear
to behave like the estimator based on the "real counterfactuals" and perform
better than the uncorrected estimators.

Table 2: Mean MISE $\pm$ standard deviation for corrected estimators and estimator obtained using "real counterfactuals" | Corrected estimator | Estimate
---|---|---
| estimator | using "real
|  with $\widehat{\Sigma}$ |  with $\Sigma$ | counterfactuals"
1 Covariate
$\Sigma_{k,k}=0.4$ | $0.034\pm 0.038$ | $0.034\pm 0.038$ | $0.034\pm 0.034$
$\Sigma_{k,k}=0.8$ | $0.136\pm 0.131$ | $0.136\pm 0.131$ | $0.060\pm 0.066$
$\Sigma_{k,k}=1.2$ | $0.412\pm 0.384$ | $0.412\pm 0.384$ | $0.186\pm 0.168$
$\Sigma_{k,k}=1.6$ | $0.787\pm 0.730$ | $0.783\pm 0.726$ | $0.315\pm 0.436$
3 Covariates
$\Sigma_{k,k}=0.4$ | $0.040\pm 0.046$ | $0.040\pm 0.046$ | $0.041\pm 0.044$
$\Sigma_{k,k}=0.8$ | $0.224\pm 0.237$ | $0.224\pm 0.237$ | $0.170\pm 0.237$
$\Sigma_{k,k}=1.2$ | $0.518\pm 0.548$ | $0.516\pm 0.546$ | $0.283\pm 0.379$
$\Sigma_{k,k}=1.6$ | $1.381\pm 1.299$ | $1.371\pm 1.291$ | $0.652\pm 0.956$
6 Covariates
$\Sigma_{k,k}=0.4$ | $0.075\pm 0.08$ | $0.075\pm 0.079$ | $0.071\pm 0.060$
$\Sigma_{k,k}=0.8$ | $0.320\pm 0.360$ | $0.320\pm 0.359$ | $0.330\pm 0.471$
$\Sigma_{k,k}=1.2$ | $0.901\pm 1.012$ | $0.898\pm 1.011$ | $0.775\pm 1.033$
$\Sigma_{k,k}=1.6$ | $2.448\pm 2.316$ | $2.437\pm 2.306$ | $1.563\pm 2.146$

##  6 Application

Now, let us apply method to real life data from the open access MIMIC-III
database [17]. This database contains deidentified data of patients admitted
to the critical care unit of the Beth Israel Deaconess Medical Center in
Boston, Massachusetts between 2001 and 2012. Over this period, 53,423 distinct
admissions of adult patients were identified. Data includes demographics data,
laboratory and microbiology test results, diagnosis and procedure codes,
bedside monitoring data (vital signs, waveforms, trends…), fluids,
medications.

The focus here is on patients with sepsis. These patients are identified
according to the algorithm defined on [26], based on specific ICD-9 codes and
procedures codes. We try to estimate the effect of vasopressor on death. The
following vasopressors, norepinephrine, epinephrine, phenylephrine,
vasopressin, dopamine, isuprel, are sought in the fluids administered to
patients during their stay. Several prognostic factors are consider, systolic
and diastolic blood pressures as time-dependent factors and gender, age,
lactate, creatinine and SOFA score [27] at inclusion in intensive care as
baseline factors. Blood pressure data are filled in every hour. Patients are
followed from entry into intensive care until discharge or death.

Time-dependent covariates are smoothed using moving average of 5 values and
standardized to get more stable trajectories. Observed blood pressures at
inclusion for each patients are considered as baseline covariates.
Counterfactual covariate trajectories under the scenario of no treatment are
estimated according to a VAR model as presented in section Section_VAR.
Estimate starts at treatment time and stops at death or discharge. Baseline
variables mentioned above are used as adjusted covariates.

Figures 3 and 2 shows observed and counterfactual covariates from treatment
initiation. The graph is truncated at 200 hours to keep enough patients. The
red lines represent counterfactual covariates under the assumption that the
patients would not have been treated. These are unobserved trajectories. The
grey lines represent observed covariates when patients are treated.
Graphically, we see, for treated individuals, an increase in systolic and
diastolic blood pressure which is the clinical effect of vasopressors.

As a reminder, the objective here is to evaluate the effect of vasopressors on
death. ATT is estimated using Gran’s estimator and ours. We include in these
models the above mentionned covariates: smoothed standardised systolic and
diastolic blood pressure grouped into hourly intervals, gender, age, lactate,
creatinine, SOFA score, systolic and diastolic blood pressure on entry into
intensive care. These estimations are compared to naive estimation and ATE
estimation. Naive estimation is obtained using Aalen model with baseline
covariates and observed time-dependent covariates. ATE is estimated using
marginal structural additive model with the same covariates. Time-dependent
covariates where replaced by stabilised inverse probability of treatment
weights.

Cumulative treatment effect is showed in Figure 4 for vasopressor versus no
treatment. Treatment has a positive contribution to the intensity (increasing
curve) all the more important in the first days. In other words, treatment
have not a protective effect neither on the eligible population (ATE) nor on
the patients actually treated (ATT). The direction of the effect can mainly be
explained because prescription of vasopressor is directly correlated to the
severity of septic shock and thus to death. Moreover these methods do not
allow to take into account the dose which may have an effect on death.

Finally, we observe that the cumulative effect estimate from the naive model
and from the MSM seems to be linear, while this effect seems to stabilize for
the Gran’s estimators and ours. We observe that the correction has an
increasing effect with time, which is explained by the increase in the number
of patients treated over time.

Figure 2: Standardized diastolic blood pressure (black: observed trajectoried,
red: modelled trajectorie) Figure 3: Standardized systolic blood pressure
(black: observed trajectoried, red: modelled trajectorie) Figure 4: Cumulative
treatment effect estimated using Naive model, Marginal structural model Gran’s
estimator with and without correction.

##  7 Discussion

In this article, we proposed a corrected estimator of the ATT under time-
varying confounders. We were inspired by the method developed by Gran et al.
[8]. This method consists of two steps. The first one is the modelling of
time-dependent covariates if the patients had not been treated and the second
is the estimation of the treatment effect. The correction we propose is based
on the error that is made when modelling the counterfactuals. Indeed, as in
any modelling, there is a difference between the "real" data and those
modelled. Our model allows the study of:

  1. 1.

terminal outcomes such as death or repeated events/counting processes such as
rehospitalizations or asthma attacks ;

  2. 2.

multiple time varying covariates and baseline covariates.

All simulation scenarios showed that our corrected estimator provides better
estimates of the ATT than those obtained with uncorrected estimators.

Estimates of ATT and ATE answer the two different questions. The ATT estimator
addresses the causal question "Is the treatment efficient in patients in whom
it is actually prescribed?" rather than "Is the treatment efficient in the
whole population" which is addressed by ATE. It may be interesting to estimate
these two measures in the same study to get additional information about the
effect of the treatment. Their comparison can be used to assess the effect of
current treatment policy such as the choice of patients receiving treatment
and the treatment initiation time. The methods used to estimate these measures
can also be used to select optimal treatment regimes. This is particularly the
case for g-estimation method [28], Inverse Probability of Censoring Weighting
model [29] or a non-parametric structural equation model and g-computation
[30].

Our correction is introduced int the cases where the data are continuous and
the counterfactuals are modelled by a VAR(1) model. This model is simple to
implement and allows an explicit writing of the correction. Other methods
could be used such as models more specific to diseases and time-dependent
covariates such as differential equations. However, it should be noted that
the more complex the model, the more difficult it will be to make explicit the
modelling error and thus the correction of the ATT estimator.

Note that in this paper, only continuous covariates are studied. Depending on
the databases analyzed, the time-dependent covariates may be binary or count
data such as the number of nursing or general practitioner visits over each
time interval studied. Specific methods to model the counterfactuals for this
type of variable should be developed. It would also be interesting to be able
to consider different types of data in the same model.

## References

  * [1] Daniel R, Cousens B, De Stavola MG et al. Methods for dealing with time-dependent confounding. Stat Med 2013; 32: 1584-1618.
  * [2] Karim ME, Petkau J, Gustafson P et al. Comparison of statistical approaches dealing with time-dependent confounding in drug effectiveness studies. Stat Methods Med Res 2018; 27: 1709-1722.
  * [3] Thomas LE, Yang S, Wojdyla D et al. Matching with time-dependent treatments: A review and look forward. Stat Med 2018; 39: 2350-2370.
  * [4] Fewell Z, Hernán MA, Wolfe F et al. Controlling for Time-dependent Confounding using Marginal Structural Models. The Stata Journal 2004; 4(4): 402-420.
  * [5] Hogan JW and Lancaster T. Instrumental variables and inverse probability weighting for causal inference from longitudinal observational studies. Stat Methods Med Res. 2004; 13(1):17-48.
  * [6] Wang A, Nianogo RA, Arah OA. G-computation of average treatment effects on the treated and the untreated. BMC Med Res Methodol. 2017 Dec;17(1):1-5.
  * [7] Vansteelandt S and Joffe M. Structural nested models and G-estimation: the partially realized promise. Statistical Science. 2014 Nov;29(4):707-31.
  * [8] Gran JM, Hoff R, Røysland K et al. Estimating the treatment effect on the treated under time-dependent confounding in an application to the Swiss HIV Cohort Study.  J R Stat Soc Ser C Stat Soc. 2018 Jan;67(1):103-25.
  * [9] Aalen OO. A linear regression model for the analysis of life times. Stat Med. 1989 Aug;8(8):907-25.
  * [10] Aalen OO, Røysland K, Gran JM et al. Can we believe the DAGs? A comment on the relationship between causal DAGs and mechanisms. Stat Methods Med Res. 2016 Oct;25(5):2294-314.
  * [11] Tchetgen EJ, Walter S, Vansteelandt S et al. Instrumental variable estimation in a survival context. Epidemiology. 2015 May;26(3):402.
  * [12] Vansteelandt S, Martinussen T, Tchetgen ET. On adjustment for auxiliary covariates in additive hazard models for the analysis of randomized experiments. Biometrika. 2014 Mar 1;101(1):237-44.
  * [13] Martinussen T, Vansteelandt S, Gerster M et al. Estimation of direct effects for survival data by using the Aalen additive hazards model. J R Stat Soc Ser B Stat Soc. 2011 Nov;73(5):773-88.
  * [14] Strohmaier S, Røysland K, Hoff R et al. Dynamic path analysis–a useful tool to investigate mediation processes in clinical survival trials. Stat Med. 2015 Dec 20;34(29):3866-87.
  * [15] Aalen OO, Cook RJ, Røysland K. Does Cox analysis of a randomized survival study yield a causal treatment effect?. Lifetime Data Anal. 2015 Oct;21(4):579-93.
  * [16] Laan MJvd and Robins JM. Unified methods for censored longitudinal data and causality. Springer series in statistics, Springer, 2003.
  * [17] Johnson AE, Pollard TJ, Shen L et al. MIMIC-III, a freely accessible critical care database. Scientific data. 2016 May 24;3(1):1-9.
  * [18] Robins JM, Rotnitzky A. Recovery of information and adjustment for dependent censoring using surrogate markers. InAIDS epidemiology. 1992 (pp. 297-331). Birkhäuser, Boston, MA.
  * [19] Hernán MA, McAdams M, McGrath N et al. Observation plans in longitudinal studies with time-varying treatments. Stat Methods Med Res. 2009 Feb;18(1):27-52.
  * [20] Gaïffas S and Guilloux A. High-dimensional additive hazards models and the Lasso. Electronic Journal of Statistics. 2012 Jan;6:522-46.
  * [21] Lütkepohl H. Introduction to multiple time series analysis. 2nd ed. Springer-Verlag, 1993.
  * [22] Lütkepohl H. New introduction to multiple time series analysis. Springer-Verlag GmbH, 2005.
  * [23] Martinussen T, Scheike TH. Dynamic regression models for survival data. Statistics for biology and health, Springer, 2006\.
  * [24] Scheike TH and Zhang MJ. Analyzing competing risk data using the R timereg package. J Stat Softw. 2011 Jan 4;38(2).
  * [25] Lewis PW and Shedler GS. Simulation of nonhomogeneous Poisson processes by thinning. Nav Res Logist Q. 1979 Sep;26(3):403-13.
  * [26] Martin GS, Mannino DM, Eaton S et al. The epidemiology of sepsis in the United States from 1979 through 2000. N Engl J Med. 2003 Apr 17;348(16):1546-54.
  * [27] Vincent JL, Moreno R, Takala J et al. The SOFA (Sepsis-related Organ Failure Assessment) score to describe organ dysfunction/failure. 1996; 707-710.
  * [28] Barrett JK, Henderson R, Rosthøj S. Doubly robust estimation of optimal dynamic treatment regimes. Stat Biosci. 2014 Nov 1;6(2):244-60.
  * [29] Shen J, Wang L, Taylor JM. Estimation of the optimal regime in treatment of prostate cancer recurrence from observational data using flexible weighting models. Biometrics. 2017 Jun;73(2):635-45.
  * [30] Díaz I, Williams N, Hoffman KL et al. Nonparametric causal effects based on longitudinal modified treatment policies. J Am Stat Assoc. 2021 Aug 29:1-6.