Considerations for Estimating Causal Effects of Informatively Timed Treatments

TL;DR

This paper highlights the importance of accounting for informative timing in treatment decisions within epidemiological studies and demonstrates how g-methods can be used to correct biases caused by timing issues.

Contribution

It formalizes the problem of informative timing, illustrates biases from ignoring it, and provides practical guidance on using g-methods for adjustment in such scenarios.

Findings

Ignoring timing can bias causal effect estimates.

Adjusting for waiting times as time-varying confounders reduces bias.

Implementation guidance with software examples is provided.

Abstract

Epidemiological studies are often concerned with estimating causal effects of a sequence of treatment decisions on survival outcomes. In many settings, treatment decisions do not occur at fixed, pre-specified followup times. Rather, timing varies across subjects in ways that may be informative of subsequent treatment decisions and potential outcomes. Awareness of the issue and its potential solutions is lacking in the literature, which motivate this work. Here, we formalize the issue of informative timing, problems associated with ignoring it, and show how g-methods can be used to analyze sequential treatments that are informatively timed. As we describe, in such settings, the waiting times between successive treatment decisions may be properly viewed as a time-varying confounders. Using synthetic examples, we illustrate how g-methods that do not adjust for these waiting times may be biased and how adjustment can be done in scenarios where patients may die or be censored in between treatments. We draw connections between adjustment and identification with discrete-time versus continuous-time models. Finally, we provide implementation guidance and examples using publicly available software. Our concluding message is that 1) considering timing is important for valid inference and 2) correcting for informative timing can be done with g-methods that adjust for waiting times between treatments as time-varying confounders.

Tables

Table 1. Table 1 : Bias and MSE of various point estimator methods for P ​ ( T 1 , 1 > 15 ) P(T^{1,1}>15) . Bias is calculated as average difference between estimate and truth value. %Bias is calculated as Bias divided by the truth value. We also report average width and coverage probability of the 95% confidence interval. Intervals were computed using bootstrap via 1000 replicates using the percentile method. Results are across 1,000 simulation data sets, each of sample size n = 2000 n=2000 . We use 500 bootstrap replicates. More details on the simulation are provided in the Appendix.

Setting	Method	%Bias	rel. MSE	Width	Coverage
No censoring	Adjusted IPTW	0.294	1	0.05	.942
	Unadjusted IPTW	8.76	2.51	0.05	.803
	Naive	45.50	38.5	0.06	.064
Censoring	Adjusted IPTW	0.03	1	0.07	.947
	Unadjusted IPTW	214.00	484.00	0.40	.119
	CC-IPTW	22.22	4.26	0.04	.012

Table 2. Table 2 : Summary of synthetic data consisting of 600 600 patients moving through two treatment courses. At each course, the table lists count and percentage (of at-risk set) of patients receiving ACT ( A k = 1 A_{k}=1 ) vs. nACT ( A k = 0 A_{k}=0 ) and patients who have normal ( L k = 0 L_{k}=0 ) versus low ( L k = 1 L_{k}=1 ) ejection fraction (EF). 416 patients reach the second treatment course. Of the remaining 184 subjects, 129 either die and 55 are censored before initiating the second course. Median waiting time until the second course is 4.2 months and median overall survival time is 7 months.

Data Summary Statistics
	Course $k=1$	Course $k=2$
	$n=600$	$n=416$
Treatment, $A_k$
ACT	$235\left(60.8\%\right)$	$256\left(61.5\%\right)$
nACT	$365\left(39.2\%\right)$	$160\left(38.5\%\right)$
EF, $L_k$
Normal	$317\left(52.8\%\right)$	$181\left(43.5\%\right)$
Low	$283\left(47.2\%\right)$	$235\left(56.4\%\right)$

Equations

$$\begin{array}{rl}{A}_1& \perp {\underline{W}}_1^{a_1},{\delta }_1^{a_1}\mid L_1=l_1,\\ A_2& \perp W_2^{a_1,a_2}\mid W_1=w_1,A_1=a_1,L_1=l_1,L_2=l_2,{\delta }_1=1.\end{array}$$

$$\begin{split}P(T^{a_{1},a_{2}}>&\tau)=\int_{\tau}^{\infty}\Big{(}\int f_{10}(t\mid l_{1},a_{1})h_{1}(l_{1})dl_{1}\\ &+\int\Big{[}\int_{0}^{t}\Big{\{}\int f_{W_{2}}(t-w_{1}\mid w_{1},\bar{a}_{2},\bar{l}_{2},\delta_{1}=1)h_{2}(l_{2}\mid a_{1},l_{1},w_{1})dl_{2}\Big{\}}f_{11}(w_{1}\mid l_{1},a_{1})dw_{1}\Big{]}h_{1}(l_{1})dl_{1}\Big{)}dt\\ \end{split}$$

$$f_{1q}(w_1\mid l_1,a_1)=\underset{dw\to 0}{\lim }\frac{P(w_1\le W_1<w_1+dw,{\delta }_1=q\mid L_1=l_1,A_1=a_1)}{dw}.$$

$$\hat{P}(T^{1,1}>\tau )=\frac{{\sum }_{i|{\delta }_{1i}=0}{A}_{1i}I(W_{1i}>\tau ){\hat{\omega }}_i+{\sum }_{i|{\delta }_{1i}=1}{A}_{1i}{A}_{2i}I(W_{1i}+W_{2i}>\tau ){\hat{\omega }}_i}{{\sum }_{i|{\delta }_{1i}=0}{A}_{1i}{\hat{\omega }}_i+{\sum }_{i|{\delta }_{1i}=1}{A}_{1i}{A}_{2i}{\hat{\omega }}_i}.$$

$${\hat{\pi }}_i=\hat{P}(A_1=1\mid L_1=l_{1i}).$$

$${\hat{\pi }}_i=\hat{P}(A_2=1\mid W_1=w_{1i},A_1=1,L_2=l_{2i},L_1=l_{1i},{\delta }_1=1)\hat{P}(A_1=1\mid L_1=l_{1i}).$$

$$\begin{array}{rl}{C}_1& \perp W_1^{a_1},W_2^{a_1,a_2}\mid A_1,L_1,\\ C_2& \perp W_2^{a_1,a_2}\mid A_1,A_2,L_1,L_2,W_1,{\delta }_1=1.\end{array}$$

$$S_{C1}(w_1\mid A_1,L_1)=P(C_1>w_1\mid A_1,L_1)>0.$$

$$S_{C2}(w_2\mid A_1,A_2,L_1,L_2,W_1,{\delta }_1=1)=P(C_2>w_2\mid A_1,A_2,L_1,L_2,W_1,{\delta }_1=1)>0.$$

$$\hat{P}(T^{1,1}>\tau )=\frac{{\sum }_{i|{\delta }_{1i}\ne 1}{A}_{1i}{\delta }_{1i}^{E}I(W_{1i}>\tau ){\hat{\omega }}_i+{\sum }_{i|{\delta }_{1i}=1}{A}_{1i}{A}_{2i}{\delta }_{2i}I(W_{1i}+W_{2i}>\tau ){\hat{\omega }}_i}{{\sum }_{i|{\delta }_{1i}\ne 1}{A}_{1i}{\delta }_{1i}^E{\hat{\omega }}_i+{\sum }_{i|{\delta }_{1i}=1}{A}_{1i}{A}_{2i}{\delta }_{2i}{\hat{\omega }}_i}.$$

$${\hat{\eta }}_i={\hat{S}}_{C1}(w_{1i}\mid A_1=1,L_1=l_{1i}).$$

$${\hat{\eta }}_i={\hat{S}}_{C2}(w_{2i}\mid A_1=1,A_2=1,L_1=l_{1i},L_2=l_{2i},W_1=w_{1i},{\delta }_1=1){\hat{S}}_{C1}(w_{1i}\mid A_1=1,L_1=l_{1i}).$$

$$r_j(V_j,Y_j;{\bar{a}}_2)=\{\begin{array}{ll}{a}_2,& \text{if}{V}_j=1,Y_j=0\\ \varnothing ,& \text{otherwise}\end{array}$$

$$\begin{split}P(T^{a_{1},a_{2}}>\tau)&=P(Y_{1}^{d_{1}}=0,Y_{2}^{\bar{d}_{2}}=0,\dots,Y_{J}^{\bar{d}_{J}}=0)\\ &=\prod_{j=1}^{J}\big{(}1-P(Y_{j}^{\bar{d}_{j}}=1\mid Y_{j-1}^{\bar{d}_{j-1}}=0)\big{)}\\ \end{split}$$

$$D_j\perp {\underline{Y}}_{j+1}^{{\bar{d}}_j}\mid {\bar{X}}_j={\bar{x}}_j,{\bar{D}}_{j-1}={\bar{d}}_{j-1},{\bar{V}}_j={\bar{v}}_j,{\bar{Y}}_j=0$$

$$A_1\perp {\underline{Y}}_2^{{\bar{d}}_1}\mid L_1=l_1,V_1=1,{\bar{Y}}_1=0$$

$$A_2\perp {\underline{Y}}_{j+1}^{{\bar{d}}_j}\mid L_1=l_1,L_2=l_2,A_1=a_1,S=j,Y_j=0$$

$$\hat{P}(Y_j^{{\bar{d}}_j}=1\mid Y_{j-1}^{{\bar{d}}_{j-1}}=0)=\frac{{\sum }_{i=1}^{n}I(D_{i,j}=d_{i,j})Y_{i,j}(1-Y_{i,j-1}){\hat{\omega }}_{i,j}}{{\sum }_{i=1}^{n}I(D_{i,j}=d_{i,j})(1-Y_{i,j-1}){\hat{\omega }}_{i,j}}$$

$${\hat{\pi }}_{i,j}=\prod _{u=1}^j\hat{P}(D_u=d_{i,u}\mid {\bar{D}}_{u-1}={\bar{d}}_{i,u-1},{\bar{X}}_u={\bar{X}}_{iu},{\bar{V}}_u={\bar{V}}_{iu},{\bar{Y}}_u=0)$$

$${\hat{\pi }}_{i,j}=\hat{P}(A_1=1\mid L_{i,1},V_1=1,Y_1=0)$$

$${\hat{\pi }}_{i,j}=\hat{P}(A_1=a_1\mid L_{i,1},V_1=1,Y_1=0)\hat{P}(A_2=a_2\mid L_{i,2},L_{i,1},A_1=a_1,S=j,V_1=1,Y_j=0)$$

$$P(T^{a_1,a_2}>\tau )={\int }_{\tau }^{\infty }{f}_{T^{a_1,a_2}}(t)dt$$

$$\begin{array}{rl}{f}_{T^{a_1,a_2}}(t)& =\sum _{u=0}^1{f}_{T^{a_1,a_2}}(t,{\delta }_1^{a_1}=u)\\ & =f_{T^{a_1,a_2}}(t,{\delta }_1^{a_1}=0)+f_{T^{a_1,a_2}}(t,{\delta }_1^{a_1}=1)\end{array}$$

$$\begin{array}{rl}{f}_{T^{a_1,a_2}}(t,{\delta }_1^{a_1}=0)& =\int f_{W_1^{a_1},{\delta }_1^{a_1}\mid L_1}(t,{\delta }_1^{a_1}=0\mid L_1=l_1)h_1(l_1)d{l}_1\\ & =\int f_{W_1^{a_1},{\delta }_1^{a_1}\mid L_1}(t,{\delta }_1^{a_1}=0\mid A_1=a_1,L_1=l_1)h_1(l_1)d{l}_1\\ & =\int f_{W_1,{\delta }_1\mid L_1}(w_1,{\delta }_1=0\mid A_1=a_1,L_1=l_1)h_1(l_1)d{l}_1\end{array}$$

$$\begin{array}{rl}{f}_{T^{a_1,a_2}}(t,{\delta }_1^{a_1}=0)& =\int f_{10}(t\mid A_1=a_1,L_1=l_1)h_1(l_1)d{l}_1\end{array}$$

$$\begin{split}f_{T^{a_{1},a_{2}}}(t,\delta_{1}^{a_{1}}=1)&=\int f_{T^{a_{1},a_{2}}}(t,\delta_{1}^{a_{1}}=1\delta_{1}^{a_{1}}=1\mid L_{1}=l_{1})h_{1}(l_{1})dl_{1}\\ &=\int f_{T^{a_{1},a_{2}}}(t,\delta_{1}^{a_{1}}=1\mid A_{1}=a_{1},L_{1}=l_{1})h_{1}(l_{1})dl_{1}\\ &=\int\Big{[}\int_{0}^{t}f_{T^{a_{1},a_{2}}}(t\mid W_{1}=w_{1},A_{1}=a_{1},L_{1}=l_{1},\delta_{1}=1)\\ &\ \ \ \times f_{W_{1},\delta_{1}\mid A_{1},L_{1}}(w_{1},\delta_{1}=1\mid A_{1}=a_{1},L_{1}=l_{1})dw_{1}\Big{]}h_{1}(l_{1})dl_{1}\\ \end{split}$$

$$\begin{split}f_{T^{a_{1},a_{2}}}(t,\delta_{1}^{a_{1}}=1)&=\int\Big{[}\int_{0}^{t}f_{T^{a_{1},a_{2}}}(t\mid W_{1}=w_{1},A_{1}=a_{1},L_{1}=l_{1},\delta_{1}=1)\\ &\ \ \ \times f_{11}(w_{1}\mid A_{1}=a_{1},L_{1}=l_{1})dw_{1}\Big{]}h_{1}(l_{1})dl_{1}\\ \end{split}$$

$$\begin{array}{rl}{f}_{T^{a_1,a_2}}(t\mid w_1,a_1,l_1,{\delta }_1=1)& =\int f_{T^{a_1,a_2}}(t\mid w_1,a_1,{\bar{l}}_2,{\delta }_1=1)h_2(l_2\mid a_1,l_1,w_1)d{l}_2\\ & =\int f_{T^{a_1,a_2}}(t\mid w_1,{\bar{a}}_2,{\bar{l}}_2,{\delta }_1=1)h_2(l_2\mid a_1,l_1,w_1)d{l}_2\\ & =\int f_W(t-w_1\mid w_1,{\bar{a}}_2,{\bar{l}}_2,{\delta }_1=1)h_2(l_2\mid a_1,l_1,w_1)d{l}_2\end{array}$$

$$\begin{split}f_{T^{a_{1},a_{2}}}(t,\delta_{1}^{a_{1}}=1)&=\int\Big{[}\int_{0}^{t}\Big{\{}\int f_{W_{2}}(t-w_{1}\mid w_{1},\bar{a}_{2},\bar{l}_{2},\delta_{1}=1)h_{2}(l_{2}\mid a_{1},l_{1},w_{1})dl_{2}\Big{\}}\\ &\ \ \ \times f_{11}(w_{1}\mid a_{1},l_{1})dw_{1}\Big{]}h_{1}(l_{1})dl_{1}\\ \end{split}$$

$$\begin{split}f_{T^{a_{1},a_{2}}}(t)&=\int f_{10}(t\mid A_{1}=a_{1},L_{1}=l_{1})h_{1}(l_{1})dl_{1}\\ &\ \ \ +\int\Big{[}\int_{0}^{t}\Big{\{}\int f_{W_{2}}(t-w_{1}\mid w_{1},\bar{a}_{2},\bar{l}_{2},\delta_{1}=1)h_{2}(l_{2}\mid a_{1},l_{1},w_{1})dl_{2}\Big{\}}f_{11}(w_{1}\mid a_{1},l_{1})dw_{1}\Big{]}h_{1}(l_{1})dl_{1}\\ \end{split}$$

Full text

1 Introduction

Epidemiological studies often target the causal effect of a sequence of treatment decisions on an event-time outcome. This is especially common in studies of chronic diseases, which are managed with a series of treatment decisions over time rather than administering a single point-treatment. These decisions are not made at fixed, pre-specified followup times. Rather, the timing varies across subjects in ways that are informative of subsequent treatment decisions and outcomes. G-methods provide valid causal inference for time-varying treatments in the presence of time-varying confounding [4, 19, 23, 24]. However, it is not obvious how these methods can be used in settings with informative, irregularly timed treatments. The purpose of this manuscript is to 1) illustrate biases due to ignoring timing and 2) describe how g-methods can be implemented to adjust for it.

As our central motivating example, will consider recent studies assessing the effects of chemotherapy treatments on survival among patients with pediatric acute myeloid leukemia (AML) [21, 6, 7]. These are secondary analyses of the AAML1031 and AAML1053 trials in which enrolled patients moved through multiple treatment courses. At the start of each course, covariates were measured and used to inform the decision to administer anthracyclene chemotherapy (ACT) or some other non-anthracycline chemotherapy (nACT) at that course. The scientific quantity of interest was the proportion of patients that would have survived at least 3 years under various ACT administration strategies. Note ACT inclusion was not randomized in the trial. Due to cardiotoxicity of ACT, ejection fraction (EF) was monitored at each course. Patients with high EF were more likely to be given ACT. An additional complication is that patients did not and, indeed, could not initiate each course at pre-set, protocol-defined times. This is because chemotherapy is physically taxing and a subsequent course could only be initiated when a subject was deemed sufficiently stable. According to the AAML1031 trial protocol [1], “Patients should progress to the next course of therapy as soon as clinically acceptable…Patients should have time for full or partial recovery from any toxicities experienced in the prior course(s).” Subjects recovering faster between courses may be systemically better off and tolerate their initial treatment well. As a result, they may both have better future potential survival and be less likely to change their previous treatment. In this sense, the waiting times between treatments can be considered as a kind of time-varying confounder - which we call “informative timing”.

As another example, consider the use of bisphosphonates in patients diagnosed with osteoporosis [10, 9]. Bisphosphonates are a class of drugs used to treat osteoporosis. Since long-term use is thought to lead to increased adverse events, it is common to switch patients to a second-line therapy such as teriparatide or denosumab [18]. We can formalize this as a two-course treatment - a first-line treatment decision followed by a second-line treatment decision. The target of inference is cumulative risk of hip fracture or atypical femoral fracture. As in the AML example, different patients switch after different durations of bisphosphonate therapy. The length of time on first-line therapy may be informative of which drug they switch to as well as outcomes.

Using AAML1031 as a motivator, we begin by formalizing the problem of informative timing as a type of time-varying confounding. Then, we describe continuous-time approaches for adjusting for informative timing via inverse-probability of treatment weighting (IPTW). Continuous-time causal approaches have been used for structural nested modeling [17], structural nested failure time modeling [16], and marginal structural modeling [25, 26]. In epidemiological research, discrete-time causal models are also common and some work has investigated the impact of discretizing continuous data [30, 5, 27]. We build on this literature in Section 5 by providing an equivalent discrete-time framing and adjustment for the issue of informative timing.

2 Potential Outcomes and Informative Timing

Consider a case with $K$ treatment courses indexed by $k=1,2,\dots,K$. For simplicity of exposition, we consider a case with just $K=2$ courses each of which have two possible treatment options which we denote as $A_{k}\in\{0,1\}$, respectively. For example, patients undergo $K=2$ courses of chemotherapy and at each course there are two possible chemotherapy agents, ACT ($A_{k}=1$) or nACT ($A_{k}=0$). Following the template of the AAML1031 trial, start of followup (time-zero) is the time of the first ($k=1$) treatment course for each patient. Let $L_{k}$ denote the set of confounders measured at course $k$ just before the treatment decision $A_{k}$. For concreteness, this could include EF in the AAML1031 context. At course $k=1$, $L_{1}$ is measured and used to inform ACT assignment, $A_{1}$. After this course, patients take some time to recover and initiate their second course $W_{1}$ time units later. If they reach course $k=2$, confounders $L_{2}$ are once again measured and, along with available information $(A_{1},L_{1},W_{1})$, is used to inform ACT decision $A_{2}$. After the initiation of the second treatment course, there is some waiting time $W_{2}$ until death, measured in units of time since start of course $k=2$. Of interest is a survival time outcome, $T>0$. For a subject who survived through two treatment courses, the survival time is the sum of the waiting times $T=W_{1}+W_{2}$. This process is illustrated in the the directed acyclic graph (DAG) in Figure 1(a). The backdoor path $A_{2}\leftarrow W_{1}\rightarrow W_{2}$ is what we refer to as “informative timing.” That is, two patients who received $A_{2}=1$ vs. $A_{2}=0$ in their second treatment course at different times $W_{1}$ are, in general, not exchangeable with respect to subsequent survival time $W_{2}$ - even if the match exactly on $A_{1}$, $L_{1},$ and $L_{2}$. This may be because subjects who move through their sequence faster (shorter $W_{1}$) are seen to, say, tolerate their previous treatment especially well and so are more likely to stay on their initial treatment.

In practice, deaths may occur before completion of the sequence so that $W_{1}$, is the first of either 1) the waiting time until death, $W_{T1}$, or 2) waiting time until the second treatment, $W_{A1}$. We denote this as $W_{1}=\min(W_{T1},W_{A1})$. The corresponding event type indicator is $\delta_{1}=I(W_{A1}<W_{T1})$, where $\delta_{1}=0$ indicates death following the first treatment course. Until Section 4, we assume there is no censoring. The observed data is given by $D=\{l_{1i},a_{1i},w_{1i},\delta_{1i}\}_{i=1}^{n}\cup\{l_{2i},a_{2i},w_{2i}\}_{i|\delta_{i1}=1}$. As is standard in this literature, we use the overbar notation to denote history (e.g. $\bar{L}_{3}=(L_{1},L_{2},L_{3})$) and underbar notation to denote future (e.g. $\underline{L}_{3}=(L_{3},L_{4},\dots,L_{K})$).

To formalize the causal structure, we define the potential outcome. Let $W_{1}^{a_{1}}=\min(W_{A1}^{a_{1}},W_{T1}^{a_{1}})>0$ be the potential waiting time until the first event after the initiation of the first treatment course had a subject received first-course chemo $a_{1}$. Similarly, let $\delta_{1}^{a_{1}}=I(W_{A1}^{a_{1}}<W_{T1}^{a_{1}})$ be the potential event type - i.e. whether a patient would have died or not before their second course. Thus, the potential survival time for a subject who potentially would have underwent two courses is the sum of the waiting times $T^{a_{1},a_{2}}=W_{1}^{a_{1}}+W_{2}^{a_{1},a_{2}}$. For a subject who would have potentially died before reaching the second treatment $T^{a_{1},a_{2}}=W_{1}^{a_{1}}$. More compactly, $T^{a_{1},a_{2}}=W_{1}^{a_{1}}+\delta_{1}^{a_{1}}W_{2}^{a_{1},a_{2}}$.

It is helpful to note that the treatment sequence $(a_{1},a_{2})$ can be viewed as being generated by an underlying dynamic treatment rule (DTR), $r$, that maps the available history at course $k$, $H_{k}=(\bar{A}_{k-1},\bar{W}_{k},\bar{L}_{k},\bar{\delta}_{k})$, to either ACT or nACT. For example, $(a_{1}=1,a_{2}=1)$ may be generated by always treating with ACT in the first course $a_{1}=r(H_{1})=1$ for any $H_{1}$. Then, treating with $a_{2}=r(H_{2})=1$ only if $\delta_{2}=1$ - i.e. if they survive long enough to initiate the second treatment course - and do nothing otherwise. So, contrasts of treatment strategies $(a_{1},a_{2})$ can be seen as emulations of a target trial where we randomly assign patients at baseline to follow different rules, $r$, for as long as they remain alive. The timing at which this rule is executed to get $a_{2}=r(H_{2})$ occurs as it naturally would under previous treatment $a_{1}=r(H_{1})$. This mirrors the actual AAML1031 trial protocol which states that patients “progress to the next course of therapy as soon as clinically acceptable.” Throughout, the target estimand is $P(T^{1,1}>\tau)$ - the proportion of the target population that would have survived at least $\tau$ time-units had they followed treatment sequence $\bar{a}_{2}=(1,1)$.

3 Adjustment for Informative Timing

Since $P(T^{1,1}>\tau)$ involves a distribution of potential survival time, causal identification assumptions are required to map this to the distribution of the observed data. With informative timing, one such required assumption is sequentially exchangeable first and second course treatment decisions

[TABLE]

The first statement above asserts that once we condition on $L_{1}$, it should not be the case that, for example, treated patients are more likely to die before subsequent treatment than untreated patients. The second statement says that among those who have reached the second treatment ($\delta_{1}=1$), conditioning on timing ($W_{1}$), previous treatment ($A_{1}$), and confounder history $(L_{1},L_{2})$, is sufficient to render the second treatment assignment independent of potential waiting time until death. Positivity assumptions are also required. Specifically, we require that $P(A_{1}=1\mid l_{1})>0$ for all $l_{1}$ in the support support and that $P(A_{2}=1\mid w_{1},a_{1},\bar{l}_{2},\delta_{1}=1)>0$ for all histories $(w_{1},a_{1},\bar{l}_{2})$ in the support. Intuitively, this second positivity assumption asserts that for each possible value of the history vector $(w_{1},a_{1},\bar{l}_{2})$ among those at-risk for second treatment ($\delta_{1}=1$), there must be some patients treated with $A_{2}=1$ and some with $A_{2}=0$. Otherwise, we would not be able to learn causal effects within that level of the history ($w_{1},a_{1},\bar{l}_{2}$). Since $w_{1}$ continuous and $\bar{l}_{2}$ may include continuous covariates (as in most application), modeling/smoothing is unavoidable since it is rare a treated and untreated unit will have the same exact history ($w_{1},a_{1},\bar{l}_{2}$). Given these assumptions, a g-formula for the quantity of interest is given by

[TABLE]

This is a special case of Equation 8.99 in [28] - a derivation is provided in the appendix. Note it is a sum of two terms. The second captures survival from patients who survive through both treatment courses while the first term captures survival for patients who never reach their second course. $h_{1}$ and $h_{2}$ denote the density functions of the time-varying covariate while $f_{W_{2}}$ is the density of the waiting time from course 1 to death. Unlike in typical settings with regularly spaced, protocol-specified treatment times, the g-formula above contains the sub-density function of the waiting time until each event of type $q\in\{0,1\}$,

[TABLE]

The appearance of sub-density functions, commonly encountered in the competing risk literature, is a consequence of death and subsequent treatment “competing” to be the next event. They can be modeled directly (e.g. using accelerated failure time models) or indirectly by modeling the corresponding cause-specific hazards [14, 2]. An estimate $\hat{P}(T^{1,1}>\tau)$ can obtained by plugging in model estimates, $\hat{h}_{1},\hat{h}_{2},\hat{f}_{W_{2}},\hat{f}_{1q}$ into the g-formula and evaluating the integrals via Monte Carlo - an approach known as “g-computation.” This outcome modeling approach is especially complex as a competing risk setup is required to capture patients’ transition through a sequence of treatments and death in continuous time. In contrast, the Hájek inverse-weighted estimator in (2) is much simpler as it only requires models for the two treatment decisions:

[TABLE]

Notice there are two terms in both the numerator and denominator. The first sums over contributions from patients who only received one treatment before death and the second sums over contributions from patients who received two treatments before death. The weights are given by $\hat{\omega}_{i}=1/\hat{\pi}_{i}$ where for subjects with $\delta_{1i}=0$, the propensity score estimate is

[TABLE]

For subjects with $\delta_{1i}=1$ the propensity score estimate is

[TABLE]

Here, the conditional probability of the second treatment is estimated using only patients who survived long enough to receive it ($\delta_{1}=1$), and this is again justified by the second line of the exchangeability condition in (1). This estimator allows for the treatment probability at the second course to depend on time elapsed since the first treatment, $W_{1}$. The top panel of Table 1 displays results from a simulation in which patients move through $K=2$ courses and initiate their second course at some random time, $W_{1}$, that depends on prior treatment and covariates and also informs subsequent treatment, $A_{2}$, and survival, $W_{2}$. The “Unadjusted IPTW” estimator refers to the estimator in (2) where the propensity score weights do not adjust for $W_{1}$. Notice that it produces a biased points estimate (8.76%) and 95% intervals with under-coverage (80.3%). On the other hand, adjusting for $W_{1}$ in (2) - labeled “Adjusted IPTW” - alleviates this bias. The “Naive” estimator that subsets to subjects who survived through both courses and were treated in both is obviously severely biased. The bias in the naive estimator shows the importance of properly accommodating deaths between treatments as well as adjustment for informative timing.

4 Censoring Between Treatment Courses

In the case where patients drop out of a study, survival time is censored. Accordingly, the observed waiting time after the first treatment is $W_{1}=\min(W_{T1},W_{A1},C_{1})$, where $C_{1}$ is the waiting time until censoring. Now, there are three possible events of type $\delta_{1}\in\{1,0,-1\}$, indicating a subsequent treatment event, death event, or censoring event respectively. Similarly, for those who received a second treatment, the observed waiting time is the first of either the waiting time from second treatment to death, $W_{T2}=T-W_{1}$, or the waiting time from second treatment to censoring, $C_{2}$. We denote this as $W_{2}=\min(W_{T2},C_{2})$. The associated indicator is $\delta_{2}=I(W_{T2}<C_{2})$. The observed data, then, is then $D=\{l_{1i},a_{1i},w_{1i},\delta_{1i}\}_{i=1}^{n}\cup\{l_{2i},a_{2i},w_{2i},\delta_{2i}\}_{i|\delta_{i1}=1}$. We also define an indicator for either death or treatment (i.e. not censoring) between treatment 1 and 2, $\delta_{1i}^{E}=I(\delta_{1i}\neq-1)$. The methods of the previous section can be modified under an additional non-informative censoring assumption

[TABLE]

This states that, after each treatment, censored patients should not be any more or less likely to experience the next event any sooner/later than uncensored subjects conditional on the available information. In addition to exchangeability, we also require additional positivity assumptions about the censoring mechanism. Specifically, there must be some positivity probability of remaining uncensored long enough to realize a second event,

[TABLE]

for each $w_{1}$ and all levels of $(A_{1},L_{1})$. Here $S_{C1}(w_{1}\mid-)$ is the “survival” function of the waiting time to censoring. If this assumption did not hold, then there would be some subgroup of patients who would always be censored before an event - making it impossible to learn about counterfactual times within this subgroup. Similarly,

[TABLE]

among patients who reached the second treatment, there must be some chance of remaining uncensored long enough to reach a death event.

If these conditions hold, then we can estimate $P(T^{1,1}>\tau)$ as

[TABLE]

The estimator only includes contributions from patients who follow sequence $(A_{1}=1,A_{2}=1)$ for as long as they are at-risk. That is, $A_{1}\delta^{E}=1$ only for uncensored subjects who underwent just one treatment course before death and $A_{1}A_{2}\delta_{2}=1$ only for for uncensored subjects who underwent two treatment courses. The weights must now correct for the fact that this subset of patients are not representative of the target population both in terms of treatment and censoring probability and so they become $\hat{\omega}_{i}=1/(\hat{\pi}_{i}\cdot\hat{\eta}_{i})$. Here, $\hat{\pi}_{i}$ is the estimated propensity score as given in Equations (3) and (4) in the previous section. The term $\hat{\eta}_{i}$ is the estimated probability of subject $i$’s observed censoring history. For subjects who underwent one treatment before death, this is given by the conditional probability that the waiting time to censoring exceeds the observed waiting time, $w_{1i}$,

[TABLE]

Similarly, for subjects who underwent two treatment courses, this is given by the conditional probability that the waiting time from treatment 1 to censoring exceeds the observed waiting time until treatment 2, multiplied by the conditional probability of the waiting time from treatment 2 to censoring exceeding the observed waiting time from treatment 2 to death

[TABLE]

Note that computation of these weights require estimates of the censoring probabilities, $\hat{S}_{C1}$ and $\hat{S}_{C2}$. The bottom panel of Table 1 displays results from a simulation in which patients may die or be censored before undergoing a second treatment. Both “Adjusted IPTW” and “Unadjusted IPTW” methods implement (6). The former adjusts for $W_{1}$ when computing $\hat{\pi}_{i}$ and $\hat{\eta}_{i}$, whereas the latter does not. Again, failing to adjust for timing results in bias and under-coverage. The method “CC-IPTW” method implements the estimator (2) using only subjects who are uncensored. This essentially fails to weight for censoring and thus leads to biased estimates. The bias of CC-IPTW emphasizes the need to handle censoring by including contributions of both patients who die after completion of the sequence as well as those who die in between treatments. Ignoring any of these complexities may result in poor estimates.

5 Adjustment with Discrete-Time Models

We now describe the above concepts within a discrete-time framework. We show that causal estimands, identification results, and the need to adjust for informative timing are agnostic to the choice of discrete versus continuous time modeling and so the main takeaways are analogous. For simplicity, we consider a case with no censoring until the end.

Again, the target of inference is $P(T^{1,1}>\tau)$. We start by partitioning the time window $[0,\tau]$ into $J$ disjoint, equal-length intervals indexed by $j=1,2,\dots,J$ with interval $j=1$ being the baseline. Let $D_{j}\in\{\emptyset,0,1\}$, be the treatment decision made in interval $j$ where $\emptyset$, 0, and 1 denote no treatment, decision to treat with nACT, and decision to treat with ACT, respectively. We let $V_{j}\in\{0,1\}$ be an indicator of whether a treatment course occurred at interval $j$ and we let $X_{j}$ denote covariates measured in interval $j$. We will let $X_{j}=\emptyset$ denote that covariates were not assessed in interval $j$. Finally, we let $\{Y_{j}\}_{j=1}^{J}$ be a monotone indicator trajectory which discretizes the continuous survival time outcome, $T$. Specifically, $Y_{j}=0$ for all intervals up to the interval of death and $Y_{j}=1$ for all intervals on and after the interval of death, as shown in Figure 1(a)

The continuous-time analysis of the two-course treatment sequence of the previous sections correspond to an analysis of these discrete-time variables under the following constraints: First, covariates are assessed and treatments decisions are made only during a treatment course: i.e., $V_{j}=0\implies D_{j}=\emptyset$ and $X_{j}=\emptyset$ with probability 1 for all $j$. Second, all subjects initiate the first course at baseline so that $V_{1}=1$ for all subjects. Third, at a single, random future interval $S\in\{2,\dots,J+1\}$, $V_{S}=1$. That is, the treatment of interest consists of only two courses (one post-baseline course) and different patients may initiate this second course at different intervals, $S$. The value $S=J+1$ indicates that the patient did not initiate the second course within $\tau$. Thus there is a 1-1 correspondence between the random variable $S$ and the length-$J$ random vector $\bar{V}_{J}$. For example, $S=5$ implies $\bar{V}_{J}$ contains a 1 in the $1^{st}$ and $5^{th}$ entries, and is zero in all other entries - and vice versa. To link this back to the continuous-time notation explicitly, the first course is initiated at the first interval $j=1$ so $V_{j}=1$ and we measure covariates $X_{1}=L_{1}$ and make decision $D_{1}=A_{1}\in\{0,1\}$. The second course is initiated at $W_{1}$ which is in interval $S$ so that $V_{S}=1$. At this interval, covariates are collected $X_{S}=L_{2}$ and a second treatment decision is made, $D_{S}=A_{2}\in\{0,1\}$. This relationship between the discrete and continuous-time variables is illustrated in Figure 1(b).

Now we consider the discrete-time potential outcomes. Just as in Section 2, the vector of observed decisions, $\bar{d}_{J}$, can be framed as the output of a dynamic treatment rule (DTR), $r_{j}(V_{j},Y_{j};\bar{a}_{2})$, where treatment course initiation and survival status are tailoring variables. Specifically, at the first interval, $V_{1}=1$ and $Y_{1}=0$ and so $d_{1}=r_{1}(1,0;\bar{a}_{2})=a_{1}$. For later intervals $j>1$ the rule $r_{j}$ makes decision $a_{2}$ if a patient is alive and has their second course at interval $j$. If they are either dead or a treatment course does not occur at interval $j$, then no decision is made. Formally,

[TABLE]

In the DTR literature, $\bar{a}_{2}$ would be called the “parameter” of the rule. The sequence of rules $r=(r_{1},r_{2},\dots,r_{J})$ explicitly connects the treatment decisions in continuous-time $(a_{1},a_{2})$ with the discretized decision sequence $\bar{d}_{J}=(d_{1},d_{2},\dots,d_{J})$. That is, we should really write $\bar{d}_{j}^{r}$ as it is the sequence of decisions made by following rule $r$ up through interval $j$, but we avoid explicit indexing by $r$ for compactness. The DTR formulation clearly distinguishes between the aspects of the decision rule we can specify in a trial protocol from what we cannot specify in a protocol. In trials like AAML1031, while we can specify which decisions are made, $\bar{a}_{2}$, we cannot specify when they are made, $S$.

We let $Y_{j}^{\bar{d}_{j}}\in\{0,1\}$ denote the potential death indicator with zero indicating a patient would have remained alive through interval $j$ and one denoting dead in interval $j$ had they followed treatment sequence $\bar{d}_{j}$. The terminal nature of death implies that $Y_{j-1}^{\bar{d}_{j-1}}=1\implies Y_{j}^{\bar{d}_{j}}=1$ with probability 1 for all $j$. The target estimand of $P(T^{1,1}>\tau)$ in this discrete-time setting is given by

[TABLE]

where $P(Y_{j}^{\bar{d}_{j}}=1\mid Y_{j-1}^{\bar{d}_{j-1}}=0)$ is the marginal discrete-time hazard as discussed in [12]. Identification of the discrete-time hazard requires a condition analogous to (1)

[TABLE]

This is the usual sequential exchangeability assumption where $X_{j}$ and $V_{j}$ are treated as observed time-varying covariates. Graphically, this asserts that $(\bar{X}_{j},\bar{D}_{j-1},\bar{V}_{j},\bar{Y}_{j})$ d-separates $D_{j}$ and $\underline{Y}_{j+1}^{\bar{d}_{j}}$ on a single-world intervention graph (See Appendix A4). The assumption holds trivially except at $j=1$ and $j=S$. This is because $V_{j}=0$ for all other $j$, and so $D_{j}=\emptyset$ deterministically. At $j=1$, everyone initiates their first course $V_{1}=1$ and the covariates consist of $X_{1}=L_{1}$ and so (9) simplifies to

[TABLE]

That is, among those at-risk at the first interval ($Y_{1}=1$) who initiate the first treatment course at the first interval ($V_{1}=1$) the treatment decision, $D_{1}\in\{0,1\}$ is unrelated to future potential survival trajectory. Replacement of $D_{1}$ with $A_{1}$ is justified since they are linked by $r$, $D_{1}=r_{1}(1,0)=A_{1}$. Notice that this is equivalent to the first line of (1). Similarly, at the interval of the second treatment course, $j=S$, the covariate history $\bar{X}_{S}$ only consists of two non-null values, $(L_{1},L_{2})$. And $\bar{D}_{j-1}$ only consists of one non-null value, $A_{1}$. Thus, (9) simplifies to

[TABLE]

Replacing conditioning on $\bar{V}_{j}$ in (9) with $S$ is justified by the aforementioned 1-1 correspondence between them. This conditional independence statement asserts that among those at risk ($Y_{j}=0$) patients with different $A_{2}$ who have the same elapsed time since previous treatment, $S$, same covariates $L_{1}$, $L_{2}$, and same prior treatment $A_{1}$ are exchangeable. This is analogous to the second line of (1). Just as in (1) we require conditioning on the time elapsed since the first course $W_{1}$ to adjust for informative timing, here we require conditioning on the discrete-time analogue - the number of intervals elapsed since the last course $S$.

Under these identification assumptions (along with analogous sequential consistency and positivity assumptions), the Hájek inverse-probability of treatment weighted estimator for the discrete-time hazard at interval $j$ is

[TABLE]

See Equation 39 in [29] which presents this estimators in settings with competing events. Extensions to the DTR settings are given in [20] and [28] (Section 8.3). The indicator $I(D_{i,j}=d_{i,j})$ selects person-period contributions with decisions, $D_{i,j}$, consistent with the DTR $r$. The time-varying weights are $\hat{\omega}_{i,j}=1/\hat{\pi}_{ij}$, where the estimated time-varying propensity score is

[TABLE]

In our setting, this formula simplifies. Since $D_{j}=\emptyset$ whenever $V_{j}=0$ with probability 1 and $V_{j}=0$ at only $j=1$ and $j=S$, the product in $\hat{\pi}_{i,j}$ is comprised only of at most two factors, not $j$ factors. This is because all other factors in the product of $\hat{\pi}_{i,j}$ are $P(D_{j}=\emptyset\mid-,V_{j}=0)=1$. The weight for person-time intervals $1\leq j<S$ is

[TABLE]

For subject $i$ who has had both treatment courses by interval $k$, the weight person-time interval $j\geq S$ is

[TABLE]

Notice that these are the same (up to a discretization) as (3) and (4). Just as in (4) we control for the time elapsed since the first course $W_{1}$, here we control for the number of intervals elapsed since the last course, $S$. While $\hat{P}(Y_{j}^{\bar{d}_{j}}=1\mid Y_{j-1}^{\bar{d}_{j-1}}=0)$ can be computed directly from the formula above for each $j$, it is commonly computed via a saturated marginal structural cox model (MSM) [12]. That is, the discretized data are arranged in long person-interval format and a logistic regression model with interval-specific intercepts is fit for $Y_{j}$ using the $\hat{\omega}_{i,k}$ as weights. In the case of censoring, the time-varying weights $\hat{\omega}_{i,j}$ will be the inverse of the product of both the treatment probabilities and the probability of remaining uncensored. The latter is estimated via another pooled logistic regression on a monotone censoring indicator [12].

To compare both approaches, we simulated a dataset mimicking AAML1031, summarized in Table 2. For simplicity of illustration, we have a binary time-varying covariate corresponding to normal ($L_{k}=0$) versus low ($L_{k}=1$) ejection fraction (EF). In this dataset, 416 patients reach the second treatment course. The rest die or are censored beforehand. The distribution of the waiting time until this second course is given in Figure 2(a). Figure 2(b) depicts estimates of $P(T^{1,1}>\tau)$ at points $\tau\in 0,1,2,\dots,20$. It depicts estimates using the continuous-time formula in (6) where we include waiting time $W_{1}$ as a covariate in the second-stage propensity score and censoring models (in blue) and where we exclude it (in green). In black, we show estimates using the discrete-time adjustment in Equation (10). Both adjustment procedures lead to estimates close to truth and perform similarly. On the other hand, failing to adjust for $W_{1}$ yields an estimate that is off. More details on software implementation are provided in the Appendix and the code333https://github.com/stablemarkets/informative_timing is made available online.

6 Discussion

The concepts discussed here relate to the literature on irregular visit/observation processes. Continuous-time [15, 3, 22] and discrete-time approaches [13] have been developed both within and outside of a formal causal inference setting. In the causal context with discrete-time models, these consider the observation process (roughly – but not exactly – analogous to $\{V_{j}\}$ in Section 5) to be intervenable and target causal effects of a specified joint treatment-observation strategy $(\bar{d}_{J},\bar{v}_{J})$. Person-intervals are artificially censored when they deviate from the specified strategy. The remaining person-intervals are re-weighted using the joint probability of their treatment and observation process. However in motivating settings, such as the AAML1031 trial, considered here, the timing of the treatment courses are not intervenable so such contrasts are not clinically relevant. Such weighting schemes may lead to unstable estimates as there may be very few person-intervals consistent with a specified timing strategy $\bar{v}_{J}$ in such settings. Therefore, we instead work within an alternative framework [28, 8] which views waiting times between successive events as potential outcomes of previous treatment decisions rather than a process controllable by design. It is possible, however, that these procedures could also be equivalently motivated via a stochastic intervention that sequentially sets each component of $\bar{V}_{J}$ according to its natural, time-varying distribution. A formal connection is an interesting area of future inquiry.

Appendix

A1 - Identification Proofs

Let $W_{1}^{a_{1}}=\min(W_{A1}^{a_{1}},W_{T1}^{a_{1}})$, $\delta_{1}^{a_{1}}=I(W_{A1}^{a_{1}}<W_{T1}^{a_{1}})$, and $W_{2}^{a_{1},a_{2}}$. So there are $\delta_{1}(a_{1})+1$ potential number of treatment courses. The total potential survival time is $T^{a_{1},a_{2}}=W_{1}^{a_{1}}+\delta_{1}^{a_{1}}W_{2}^{a_{1},a_{2}}$. For some fixed $\tau>0$, the estimand of interest is

[TABLE]

where $f_{T^{a_{1},a_{2}}}(t)$ is the density of the potential survival time. By the total law of probability, this can be expressed as an average over the potential number of treatments

[TABLE]

The first term captures survival past $\tau$ among those who died before intiated the second course. The second term captures survival past $\tau$ among those who survived through their second treatment course. We will identify each term in (A.11), starting from the first term. For the first term note that since $\delta_{1}^{a_{1}}=0$, $T^{a_{1},a_{2}}=W_{1}^{a_{1}}$. Below, in the second equality we average over the distribution of confounders measured at the first course, $L_{1}$. The third equality follows from invoking conditional ignorability of $A_{1}$ given $L_{1}$ to condition on $a_{1}$

[TABLE]

The last equality follows from consistency. Now, in the last line is just $f_{W,\delta_{1}\mid L_{1}}(w_{1},\delta_{1}=0\mid A_{1}=a_{1},L_{1}=l_{1})$ is just the conditional sub-density of the waiting time from the first course to death, $W_{T1}$ which in the main text we had denoted by $f_{kq}(w_{k}\mid\bar{W}_{k-1}=\bar{w}_{k-1},\bar{A}_{k}=\bar{a}_{k},\bar{L}_{k}=\bar{l}_{k})$ with $k=1$ and $q=0$. So, the first term in (A.11) is

[TABLE]

Now we identify the second term in a similar fashion starting by averaging over $L_{1}$ and invoking exchangeability and consistency of $A_{1}$:

[TABLE]

The last line follows from averaging over the distribution of waiting times which, since $\delta_{1}=1$, ranges from [math] to $t$. Now note that $f_{W_{1},\delta_{1}\mid A_{1},L_{1}}$ is just the observed data sub-density function of the waiting time from the first course to the second course. So,

[TABLE]

Averaging over the conditional density of $L_{2}$, the first factor in brackets above is

[TABLE]

In the first line $h_{2}$ is the conditional density of $l_{2}$ among those with $\delta_{1}=1$. The second line follows from sequential exchangeability. The Third line follows from the fact that $T^{a_{1},a_{2}}=W_{1}^{a_{1}}+W_{2}^{a_{1},a_{2}}$.

So the second term in (A.11) is

[TABLE]

Now that we have both terms, substituting back into (A.11) yields the expression in the main text:

[TABLE]

The survival probability is attained by integrating the density from $\tau$ to infinity,

[TABLE]

A2 - Synthetic Data Generation Details

In this section we describe the data-generating mechanisms used in the simulation analyses and the worked example. For each simmulation scenario, we generate each data sets with $n=2000$ subjects in each data set. For the worked example we generate a small data set of $n=600$ subjects.

In the scenario without censoring, the data for a single subject are generated as follows:

$L_{1}\sim Ber(.5)$ 2. 2.

A_{1}\mid L_{1}\sim Ber\Big{(}\operatorname*{expit}(1-L_{1})\Big{)} 3. 3.

Generate time to death W_{T1}\mid A_{1},L_{1}\sim Exp\Big{(}rate=\exp(-3-A_{1}+L_{1})\Big{)} and time to next treatment W_{A1}\mid A_{1},L_{1}\sim Exp\Big{(}rate=\exp(-3+A_{1}+L_{1})\Big{)}. 4. 4.

If $W_{T1}<W_{A1}$, then the patient died so set survival time $T=W_{T1}$ and stop. Move on to next patient. Otherwise, the patient survived to treatment 2. Set $W_{1}=W_{A1}$ and generate data at second treatment course: 5. 5.

L_{2}\mid L_{1},W_{1}\sim Ber\Big{(}\operatorname*{expit}(.5+.15L_{1}-.15I(W_{1}>15))\Big{)} 6. 6.

A_{2}\mid L_{2},W_{1}\sim Ber\Big{(}\operatorname*{expit}(1-L_{2}+1.5I(W_{1}>15))\Big{)} 7. 7.

W_{2}\mid A_{2},L_{2},W_{1}\sim Exp\Big{(}rate=\exp(-3+A_{2}+L_{2}-2I(W_{1}>15))\Big{)} 8. 8.

The simulated survival time is thus $T=W_{1}+W_{2}$.

Notice that subjects may die before reaching the second treatment.

In the second scenario, we allow for censoring and so the data for a single subject are generated as follows:

$L_{1}\sim Ber(.5)$ 2. 2.

A_{1}\mid L_{1}\sim Ber\Big{(}\operatorname*{expit}(1-L_{1})\Big{)} 3. 3.

Generate time to death W_{T1}\mid A_{1},L_{1}\sim Exp\Big{(}rate=\exp(-3-A_{1}+L_{1})\Big{)} and time to next treatment W_{A1}\mid A_{1},L_{1}\sim Exp\Big{(}rate=\exp(-3+A_{1}+L_{1})\Big{)}. 4. 4.

Generate time to censoring C_{1}\sim Exp\Big{(}rate=-4-A_{1}+L_{1}\Big{)} 5. 5.

If $W_{A1}>\max(W_{T1},C_{1})$, then the subject did not reach the second treatment. Their observed time is $W_{1}=\min(C_{1},W_{T1})$. The associated indicator is $\delta_{1}=0$ if they died and $\delta_{1}=-1$ if they were censored. Stop and move on to next subject. If $W_{A1}<\min(W_{T1},C_{1})$, then the patient started the second treatment. Set the waiting time between treatments to $W_{1}=W_{A1}$ and the associated indicator to $\delta_{1}=1$. Move on to second stage: 6. 6.

L_{2}\mid L_{1},W_{1}\sim Ber\Big{(}\operatorname*{expit}(.5+.15L_{1}-.15I(W_{1}>15))\Big{)} 7. 7.

A_{2}\mid L_{2},W_{1}\sim Ber\Big{(}\operatorname*{expit}(1-L_{2}+1.5I(W_{1}>15))\Big{)} 8. 8.

Simulate time from treatment 2 to death:

[TABLE] 9. 9.

Simulate time from treatment 2 to censoring:

[TABLE]

The observed waiting time is $W_{2}=\min(W_{T2},C_{2})$. The associated indicator is $\delta_{2}=I(W_{T2}<C_{2})$.

Repeating the above steps 1,000 times yields a data set for $n=1000$ subjects - who may complete the two-course sequence, die/be censored before second treatment, or be censored before death after the second treatment. Notice that censoring time is allowed to depend on waiting times and confounders. For the worked example, we analyze a data set consisting of $n=600$ patients generated under the second scenario (with censoring).

A3 - Details of Worked Example

Here we discuss the software details behind the worked example in Section 5 of the main text.

To adjust for informative timing we implement the continuous-time approach in Equation (6) with weights $\hat{\omega}_{i}=1/\hat{\pi}_{i}\hat{\eta}_{i}$ that are formed using both inverse probability of treatment and censoring weights. We compute the treatment probabilities $\hat{\pi}_{i}$ by estimating a logistic regression using the glm function in R. The predict function is used to obtain the fitted probabilities, $\hat{\pi}_{i}$, from this estimated logistic regression. The survival probabilities required to compute $\hat{\eta}_{i}$ are computed by estimating an exponential proportional hazard (PH) model for the cause-specific hazard of censoring. In the software, censoring time is treated as as the event of interest and the first of either subsequent treatment or death as the “censoring” time. For example, with these data we specified an exponential PH model for the waiting time from the first treatment to censoring, $C_{1}$, as $\lambda(c_{1})=\lambda_{0}(c_{1})\exp(\beta_{1}A_{1}+\beta_{2}L_{1})$. The corresponding survival function is given by $S_{C_{1}}(c_{1}\mid A_{1},L_{1})=\exp(-\Lambda_{0}(c_{1})\exp(\beta_{1}A_{1}+\beta_{2}L_{1}))$. Here, $\lambda_{0}(c_{1})$ is the baseline hazard of censoring and $\Lambda_{0}(c_{1})=\int_{0}^{c_{1}}\lambda_{0}(w)dw$ is the cumulative baseline hazard. We use the flexsurv package in R to estimate this PH model with, again for simplicity, a parametric specification of the baseline hazard. flexsurv obtains estimates $\hat{\Lambda}_{0}(c_{2})$, $\hat{\beta}_{1}$, and $\hat{\beta}_{2}$ via maximum likelihood estimation. After estimation, we use the predict.flexsurvreg function to plug in these estimates into $S_{C_{1}}$ to obtain $\hat{S}_{C_{1}}(c_{1}\mid A_{1i},L_{1i})$. Alternatively, the function coxph in the survival package can be used to obtain estimates $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ via partial likelihood maximization, while a nonparametric Nelson-Aalen estimate of $\Lambda_{0}(c_{1})$ can be obtained via the basehaz function. The function predict.coxph can then be used to compute $\hat{S}_{C_{1}}(c_{1}\mid A_{1i},L_{1i})=\exp(-\hat{\Lambda}_{0}(c_{1})\exp(\hat{\beta}_{1}A_{1i}+\hat{\beta}_{2}L_{1i}))$. This provides a more robust alternative.

Using these weights, we estimate $P(T^{1,1}>\tau)$ at points $\tau\in 0,1,2,\dots,20$. Figure 2(b) depicts the estimates where we include waiting time $W_{1}$ as a covariate in the propensity score and censoring models (in blue) and where we exclude it (in green). The intervals represent 95% confidence intervals calculated via 300 bootstrap replicates using the boot package. If parametric approaches are used for estimation of the nuisance parameters, the variance could also be estimated using the sandwich-variance estimator.

As a comparator we also implement the discrete-time adjustment in Equation (10). This is done as follows. First, the window from time zero to the maximum observed time is partitioned into 1 month intervals up to month 20 - which is the last point on the survival curve we wish to evaluate. We arranged these data into long person-interval format and used a pooled logistic regression, implemented in glm, to estimate the discrete-time hazard of censoring. The cumulative product of these estimates hazards at each interval gives us the estimated probability of remaining uncensored in that interval. We then multiply these with the estimated treatment probabilities given in Section 5 which are obtained by separate logistic regressions. The inverse of this product gives us the time-varying weight. Since we are interested in evaluating survival under sequence $(a_{1}=1,a_{2}=1)$, we subset to uncensored person-intervals consistent with the dynamic treatment rule, $r$, that assigns $A_{k}=1$ whenever $V_{k}=1$. Then fit a logistic regression for outcome $Y_{k}$ with interval-specific intercepts weighted by the time-varying weights in the glm function. This yields an estimate of the discrete-time hazard $P(Y_{j}^{\bar{d}_{j}}=1\mid Y_{j-1}^{\bar{d}_{j-1}}=0)$, which we then multiply up to desired interval to get the estimated survival up to that interval, as shown in Equation (8), using the cumprod function. The results are plotted in black in Figure 2(b). As can be seen, these results are nearly identical to adjusted estimates using the continuous-time model.

A4 - Graphical Perspective on Sequential Exchangeability

In Equation (9) of the main text discussing the discretized analysis, we required the following sequential conditional exchangeability assumption to hold

[TABLE]

and described how it is analogous to the required exchangeability assumption in continuous-time. Recall that $V_{j}$ is an indicator of whether a subject initiated a treatment course at interval $j$. If a course is initiated, then an active treatment, $D_{j}\in\{0,1\}$, may be given but if it is not initiated then no active treatment is given $D_{j}=\emptyset$. In this sense, $D_{j}$ can be viewed as the output of a dynamic treatment rule where $V_{j}$ is the tailoring variable.

With static treatment rules, sequential exchangeability asserts that $(\bar{X}_{j},\bar{D}_{j-1},\bar{V}_{j},\bar{Y}_{j})$ d-separates $D_{j}$ and $\underline{Y}_{j+1}^{\bar{d}_{j}^{r}}$ in a single-world intervention graph (SWIG). In the case of dynamic treatment rules, violations can occur if there are 1) unblocked backdoor paths from $D_{j}$ to future values of the tailoring variable as well and 2) unblocked backdoor paths from $D_{j}$ to future $V_{j}$. These points are discussed at length in Chapter 19 of [11]. Interested readers should see especially Section 19.5. We provide a brief summary of these points for this particular application with informative timing.

To illustrate this, Figure 3 shows a possible DAG and the corresponding SWIG template for DTR $r=\{r_{1},r_{2},\dots,r_{J}\}$. That is, the SWIG represents the joint distribution of these random variables in a world in which we intervene to set treatment decisions based on a rule $r$ with tailoring variable process ${V_{j}}$. We have omitted observed covariate process, $\{X_{j}\}$ and only depicted the graphs for $j\in\{1,2\}$ for compactness and to focus attention on ${V_{j}}$. We assume everyone enters the study alive so that $Y_{1}:=1$ for all subjects. For $j=1,2$ exchangeability asserts that

[TABLE]

Since $V_{j}$ is a tailoring variable, new dashed paths, absent on the DAG Figure 3(a), now appear on the SWIG in Figure 3(b). Note in the second line above we omit the superscripts on $D_{2}$ and $V_{2}$ which are present on the SWIG by consistency, since $d_{1}^{r}$ is conditioned on. In the presence of unmeasured common causes of treatment decision and course initiation, such as $U_{1}$, $D_{1}$ and $Y_{2}^{\bar{d}_{2}^{r}}$ are not d-separated due to open backdoor path $D_{1}\leftarrow U_{1}\rightarrow V_{2}^{d_{1}^{r}}\dashedrightarrow r_{2}(V_{2}^{d_{1}^{r}})\rightarrow Y_{3}^{d_{1}^{r},d_{2}^{r}}$. Thus, sequential exchangeability is violated since $D_{1}\not\perp(Y_{3}^{d_{1}^{r},d_{2}^{r}},Y_{2}^{d_{1}^{r}})\mid V_{1}=v_{1},Y_{1}=0$. This is true even if there is no direct effect of treatment course initiation on outcomes (arrows from $V_{j}\rightarrow Y_{j+1}$) on the DAG.