Modified treatment policies, natural value of treatment, longitudinal g-computation formula.

1 Introduction↩︎

Modified treatment policies (MTPs), also known as dynamic interventions that depend on the natural value of treatment [1]–[6], provide a flexible way to define, identify, and estimate causal effects in real-world situations such as continuous or multivariate exposures [7], limited overlap or practical positivity violations [8], grace periods or delays in treatment initiation [9], longitudinal studies with time-varying treatments, censoring, and competing risks [10], resource-constrained treatment policies [11], and policy interventions with spillovers or network interference [12]. [13] provide a thorough review of the literature on methods involving the natural value of treatment. Longitudinal modified treatment policies were originally formulated broadly enough to allow the intervention at time \(t\) to depend on the history of natural treatment values up to time \(t\) [3]. The natural treatment value at time \(t\) is the treatment that would have been observed had the intervention been stopped immediately before treatment at time \(t\) was assigned. Thus, in the original formulation, the policy at time \(t\) may depend both on the contemporaneous natural treatment value and on earlier natural treatment values. In contrast, subsequent work on inverse probability weighted [4] and targeted learning estimators [5] focused on the narrower class of policies in which the intervention at time \(t\) depends only on the contemporaneous natural treatment value. This restriction rules out policies whose treatment rule depends on the natural values of treatment at previous timepoints, such as when the intervention modifies the timing of a naturally occurring treatment process rather than only the contemporaneous treatment value. For example, in longitudinal policy evaluation, one may wish to estimate the effect of delaying the natural date of policy adoption, such as shifting the enactment time of a naloxone access law [14]. In this case, treatment at time \(t\) depends on the natural policy status at an earlier time, or equivalently on the natural adoption time. For longitudinal modified treatment policies that depend only on the contemporaneous natural value of treatment, [5] developed targeted machine learning estimators based on efficient influence function that achieve \(\sqrt{n}\)-consistency under slow convergence rates for the nuisance regressions, including the outcome regressions and treatment mechanisms. These estimators rely on the fact that, for this class of interventions, the generalized g-computation formula of [3] admits a standard iterated conditional expectation representation [15], [16]. We show that this representation fails for longitudinal modified treatment policies that depend on the natural history of treatment. We therefore introduce a modified iterated conditional expectation formula on an augmented data structure and prove that it equals the corresponding g-computation functional. We then derive the efficient influence function of the resulting nonparametric functional and characterize conditions for \(\sqrt{n}\)-consistent estimation. These results yield two extensions of the estimators of [5]: a targeted minimum loss-based estimator and a sequentially doubly robust estimator. Like their original counterparts, our estimators are asymptotically efficient under a product-rate condition requiring certain regression errors to converge to zero at \(\sqrt{n}\)-rate in \(L_2\) norm. We assess finite-sample performance in simulations and illustrate the methods in an application estimating the effect of delaying a risky pain treatment by one month on 12-month incidence of opioid use disorder or overdose diagnosis among 12,745 Medicaid-insured patients with lumbar spinal stenosis.

2 Notation and definition of causal effects↩︎

Let \(Z_1,\ldots, Z_n\) denote i.i.d. observations with \(Z=(L_1, A_1, L_2, A_2, \ldots, L_\tau, A_\tau, Y)\sim \mathsf{P}\), where \(L_t\) denotes time-varying covariates, \(A_t\) denotes a discrete exposure variable, and \(Y\) denotes an outcome of interest measured at the end of study follow-up. We let \(\mathsf{P}_n\) to denote the empirical distribution of \(Z_1,\ldots\,Z_n\), and assume \(\mathsf{P}\) is an element of the non-parametric statistical model defined as all continuous densities on \(Z\) with respect to a dominating measure \(\nu\). We let \(\mathsf{P}f = \int f(z)\mathrm{d}\mathsf{P}(z)\) for a given function \(f(z)\), and let \(\mathsf{E}\) denote the expectation with respect to \(\mathsf{P}\), i.e., \(\mathsf{E}\{f(Z)\} = \int f(z)\mathrm{d}\mathsf{P}(z)\). We use \(||f||^2\) to denote the \(L_2(\mathsf{P})\) norm \(\int f^2(z)\mathrm{d}\mathsf{P}(z)\). \(\bar X_t = (X_1,\ldots, X_t)\) will be used to denote the history of a variable, \(\underline X_t = (X_t,\ldots, X_\tau)\) to denote the future of a variable, and \(H_t = (\bar A_{t-1}, \bar L_t)\) to denote the history of all variables up until just before \(A_t\). When useful, we will also use notation \(\bar X_{t:k} = (X_t,\ldots X_k)\). We let \(\mathsf{g}_t(a_t \mid h_t)\) denote the probability mass function of \(A_t\) conditional on \(H_t=h_t\), and use calligraphic font to denote the support of a random variable, e.g., \(\mathcal{A}_t\) denotes the support of \(A_t\). Causal effects are defined relying on counterfactual logic. Counterfactual random variables will be formally defined using a non-parametric structural equation model [17]. Specifically, for each timepoint \(t\), we assume the existence of deterministic functions \(f_{L_t}\), \(f_{A_t}\), and \(f_Y\) such that \(L_t=f_{L_t}(A_{t-1}, H_{t-1}, U_{L,t})\), \(A_t=f_{A_t}(H_t, U_{A,t})\), and \(Y=f_Y(A_\tau, H_\tau, U_Y)\), where \(U=(U_{L,t}, U_{A,t}, U_Y:t\in \{1,\ldots,\tau\})\) is a vector of exogenous variables with a distribution satisfying assumptions that will be given in §3. Unless defined otherwise, sums of the type \(\sum_{k=s}^ta_k\) and products of the type \(\prod_{k=s}^ta_k\) for \(s>t\) are defined as zero and one.

2.1 Longitudinal modified treatment policies that only depend on the contemporaneous natural value of treatment↩︎

We first revisit the definition of LMTPs used in [5] to contrast it with the more general definition we use in this manuscript (§2.2). [5] defined LMTP effects in terms of hypothetical interventions where the equation \(A_t=f_{A_t}(H_t, U_{A,t})\) is removed from the structural model, and the exposure is assigned as follows.

For a user-defined function \(d_1\), at the first timepoint, \(t=1\), the post-intervention treatment is defined as a transformation \(A_1^d=d_1(A_1, L_1)\) of the natural value of treatment at that time, \(A_1\), and the history of data collected up until \(A_1\), \(L_1\). At time \(t=1\), the natural value of treatment is equal to the observed value of treatment, but this will not be the case for subsequent timepoints. The hypothetical intervention, \(d_1\), generates an intervened-on treatment at \(t=1\), which we have denoted \(A_1^d\), and counterfactual \(L_2(A_1^d)\). These are then part of the counterfactual history and the natural value of treatment at time \(t=2\), denoted with \(H_2(A_1^d) = (L_1, A_1^d, L_2(A_1^d))\) and \(A_2(A_1^d)\).

At time \(t=2\), the post-intervention treatment is defined as \(A_2^d=d_2(A_2(A_1^d), H_2(A_1^d))\). In particular, note that this intervention depends only on the contemporaneous natural value of treatment at time \(t=2\), \(A_2(A_1^d)\), and not on the natural value of treatment at \(t=1\), \(A_1\). Instead, the intervention depends on the intervened-on value of treatment at \(t=1\), namely \(A_1^d=d_1(A_1, L_1)\). This intervention at time \(t=2\) generates an intervened-on \(A_2^d\) and updated counterfactual history, which are then used to assign the post-intervention exposure at time \(t=3\) as \(A_3^d=d_3(A_3(\bar A_2^d), H_3(\bar A_2^d))\). This process is iterated forward in time until interventions at all times \(t\in[\tau]\) have been performed, at which point we obtain the counterfactual outcome \(Y(\bar A^d)=f_Y(A_\tau^d, H_\tau(\bar A_{\tau-1}^d), U_Y)\). The general definition of LMTPs that depend on the current natural value of treatment is as follows:

The data generating process that leads to the above counterfactual outcome is illustrated in Figure [fig:swig-current] with \(\tau=2\) and omitting \(L_1\). [6] provide a large list of interventions of practical interest that can be described using this framework. Importantly, delay interventions do not fit this definition. We elaborate on this point in §2.3.

2.2 Longitudinal modified treatment policies that depend on the natural history of treatment↩︎

In this paper, we generalize the above definition to allow dependence of the intervention at time \(t\) on the history of the natural value of treatment—instead of just on the contemporaneous natural value of treatment. Specifically, counterfactual outcomes are generated as follows.

The intervention at \(t=1\) is defined identical to §2.1. At time \(t=1\), the post-intervention variable is defined as a transformation \(A_1^d=d_1(A_1, L_1)\) of the natural value (which, again, is also the observed value) of treatment \(A_1\) and \(L_1\). This intervention generates a counterfactual history and treatment at time \(t=2\), denoted with \(H_2(A_1^d) = (L_1, A_1^d, L_2(A_1^d))\) and \(A_2(A_1^d)\), respectively.

At time \(t=2\), the two definitions diverge. At this timepoint the post-intervention exposure is generated as \(A_2^d=d_2(\bar A_2(A_1^d), H_2(A_1^d))\). The difference between this and the definition of \(A_2^d\) in §2.1 is the overbar for the history of the natural value \(\bar A_2(A_1^d)\), which is defined as \(\bar A_2(A_1^d)=(A_1, A_2(A_1^d))\). In particular, note that this intervention can depend on two values of treatment at \(t=1\), the natural value of treatment, \(A_1\), and the intervened value of treatment, \(A_1^d=d_1(A_1,L_1)\). This intervention at time \(t=2\) generates a counterfactual history and natural value of treatment at time \(t=3\), \(H_3(\bar A_2^d) = (L_1, A_1^d, L_2(A_1^d), A_2^d, L_3(\bar A_2^d))\) and \(A_3(\bar A_2^d)\), similar to §2.1. These variables are then used to assign the post-intervention exposure at time \(t=3\) as \(A_3^d=d_3(\bar A_3(\bar A_2^d), H_3(\bar A_2^d))\), where \(\bar A_3(\bar A_2^d)=(A_1, A_2(A_1^d), A_3(\bar A_2^d))\). This process is iterated forward in time until interventions at all times \(t\in[\tau]\) have been performed, at which point we obtain the counterfactual outcome \(Y(\bar A^d)=f_Y(A_\tau^d, H_\tau(\bar A_{\tau-1}^d), U_Y)\). Below is the formal definition, which coincides with the general definition given in §5.2 of [3]. Figure 2 gives a graphical illustration of the differences between this definition and the definition of §2.1.

2.3 Interventions that depend on the natural history of treatment↩︎

In this section we discuss two scientifically and clinically relevant interventions that require the G-LMTP framework, highlighting issues that arise when more narrowly defined LMTPs do not have access to the history of natural treatment values once a modification has occurred. Specificially, the G-LMTP framework is required for policies that depend on aspects of the natural treatment trajectory that may be overwritten by the intervention.

2.3.1 Delay interventions↩︎

In vaccine studies, one may be interested in assessing the effect of a delay in receiving a vaccine booster. Likewise, one may be interested in estimating the effect of delaying a potentially harmful intervention such as invasive mechanical ventilation in acutely ill patients. Here we briefly show that such delay policies cannot be expressed within the LMTP framework described by [5], but can be expressed by the G-LMTP framework described in this paper. Consider a simple setting with only two time points and a binary intervention at each timepoint. A delay in exposure of one timepoint would correspond to an intervention that assigns treatment at time \(t=1\) as \(A_1^d= d_1(A_1, L_1) = 0\). At time \(t=2\), the desired delay intervention should assign the treatment as \(A_2^d= A_1\). This however, is not possible in the setting depicted by Figure [fig:swig-current], as there is no dependence of the MTP on the prior natural value of treatment under that setup. In other words, \(A_2^d= d_2(A_2(A_1^d), L_2(A_1^d), A_1^d, L_1)\) must be defined solely in terms of the natural value of treatment at time 2 and the intervened history, i.e., \(d_1(A_1, L_1)\). If \(d_1\) modifies \(A_1\) (e.g., from 1 to 0), the function \(d_2\) no longer has access to the original \(A_1\), and thus cannot shift it forward in time. In contrast, under an G-LMTP, \(A_2^d\) may be defined as a function of \(A_1\) even if \(A_1^d\ne A_1\), and the policy can be explicitly constructed to reassign the original treatment at time 1 to occur at time 2 instead.

To further elaborate, consider a long-standing question in critical care regarding the effects of delaying invasive mechanical ventilation. Let \(A_t \in \{0,1,2\}\) denote the level of oxygen support at time \(t\), where \(A_t = 0\) indicates no support, \(A_t = 1\) indicates non-invasive support, and \(A_t = 2\) indicates invasive mechanical ventilation. For an arbitrary vector \(\bar a_t\), define the function\[d_t(a_t, \bar a_{t-1}) = \begin{cases} 1 & \text{if } a_t = 2 \text{ and } a_s \leq 1 \text{ for all } s < t, \\ a_t & \text{otherwise,} \end{cases}\] which represents a policy that delays mechanical ventilation by substituting it with non-invasive support unless the patient has already received invasive mechanical ventilation (IMV), where the G-LMTP is defined by \(A_t^d= d_t(A_t(\bar A_{t-1}^d), \bar A_{t-1}(\bar A_{t-2}^d))\). Consider \(\tau=3\) and a patient with observed treatment history \((A_1, A_2, A_3) = (1, 2, 2)\). At time \(t = 1\), the patient would receive \(A_1^d= 1\). At time \(t = 2\), since the natural treatment is \(A_2(1) = 2\) and the patient has not yet received invasive ventilation, the policy prescribes \(A_2^d= 1\). At time \(t = 3\), the value of \(A_3^d\) depends on the counterfactual oxygen support requirement under the delayed history: if \(A_3(1) = 2\), then the policy assigns \(A_3^d= 1\); otherwise, \(A_3^d= A_3(1)\). That is, the possible paths under this intervention are \((1, 1, 2)\), \((1, 1, 1)\), and \((1, 1, 0)\), depending on the natural value of \(A_3(1)\).

In contrast, if the intervened-on treatment is assigned as an LMTP that depends on the current natural value of treatment, i.e., \(A_t^d= d_t(A_t(\bar A_{t-1}^d), \bar A_{t-1}^d)\), corresponding to the estimators developed in [5], the intervention would effectively prevent intubation for all patients rather than intend to delay it. Continuing with the example above, a patient with observed treatment history \((A_1, A_2, A_3) = (1, 2, 2)\) would have \(A_1^d= 1\) as before. At time \(t = 2\), the LMTP rule applies \(d_2\) to \(A_2(1)=2\) and \(A_1^d= 1\), yielding \(A_2^d= 1\). However, at time \(t = 3\), the intervention uses \(A_3(\bar A_2^d) = A_3(1,1)\) and \((A_2^d, A_1^d) = (1, 1)\). The only consistent regimens under the intervention are therefore \(\bar A^d= (1, 1, 1)\) and \(\bar A^d= (1, 1, 0)\), depending on the value of \(A_3(1,1)\). In other words, the LMTP cannot detect that invasive ventilation had been indicated and subsequently overridden, and so the rule continues to suppress it indefinitely.

2.3.2 Dose escalation interventions↩︎

Another example of a class of interventions that requires access to the history of natural treatment values arises when the treatment is a dose and the scientific question concerns modification of the clinician’s natural dose-escalation behavior. For example, suppose \(A_t\) denotes the dose of a medication at visit \(t\), and consider a policy that leaves the clinician’s natural dose unchanged unless the clinician would naturally increase the dose by more than \(\delta\) since the previous visit. If the natural increase exceeds \(\delta\), the policy caps the increase at \(\delta\).

Assume that \(A_t\) is ordered and let \(\delta>0\) denote the maximum allowed increase. For arbitrary treatment values \(a_t\) and \(a_{t-1}\), define, for \(t\geq 2\), \(d_t(a_t, a_{t-1}) = a_{t-1}+\delta\) if \(a_t-a_{t-1}>\delta\), and \(d_t(a_t, a_{t-1}) = a_t\) otherwise. The rule at \(t=1\) may be defined as \(d_1(a_1,a_0)=a_1\). The corresponding G-LMTP is \(A_t^d= d_t(A_t(\bar A_{t-1}^d), A_{t-1}(\bar A_{t-2}^d))\). This policy compares the natural dose at time \(t\) with the natural dose at time \(t-1\), both evaluated under the counterfactual history induced by the previous interventions. It therefore modifies only large natural escalations, while leaving smaller natural increases, dose decreases, and dose stability unchanged.

This intervention cannot generally be represented as an LMTP that depends only on the contemporaneous natural value of treatment. Under such an LMTP, the treatment at time \(t\) may be written as \(A_t^d= d_t(A_t(\bar A_{t-1}^d), A_{t-1}^d).\) The intervened previous dose \(A_{t-1}^d\), does not generally correspond to the previous natural dose \(A_{t-1}(\bar A_{t-2}^d)\). Once a large escalation has been capped at time \(t-1\), the contemporaneous LMTP no longer has access to the dose the clinician would naturally have assigned at that time. As a result, it cannot determine whether the clinician’s next natural dose represents a large increase relative to the previous natural dose.

To illustrate, take \(\delta=10\) and suppose \(A_1=10\). Since the policy does not modify the dose at time \(1\), both the G-LMTP and the corresponding current-time LMTP assign \(A_1^d=10\). Suppose that, under this intervened history, the clinician would naturally assign \(A_2(A_1^d)=40\). Both policies would then cap the increase at \(10\), yielding \(A_2^d=20\).

Now consider the next timepoint. Suppose that, under the common intervened history through time \(2\), the clinician would naturally assign \(A_3(A_1^d,A_2^d)=45.\) At this point the two policies differ. The G-LMTP compares this current natural dose with the previous natural dose, \(A_2(A_1^d)=40\), and therefore leaves the dose unchanged because \(45-40\leq \delta\). Thus the G-LMTP assigns \(A_3^d=45\). In contrast, a current-time LMTP that only has access to the intervened history would compare the current natural dose with the previously administered intervened dose \(A_2^d=20\), and would assign \(A_3^d=30\). This is a different intervention: it caps increases relative to the previously administered intervened dose, rather than capping the clinician’s natural dose-escalation decision.

3 Identification of causal effects↩︎

For fixed values \(\bar a_t\) and \(\bar l_t\), we recursively define \(a_t^d=d_t(\bar a_t, h^d_t)\), where \(h^d_t=(\bar a_{t-1}^d, \bar l_t)\). G-LMTPs will be identified under the following assumptions. Alternative assumptions under a different causal model were given by [3] and [4].

We have the following identification theorem, which allows us to compute the parameters \(\theta_{\scriptsize lmtp}\) as a function of the observed data distribution.

The proof of this theorem is identical to the first part of the proof of Theorem 1 in [5] and is given in the Supplementary Materials. Theorem 1 in [5] additionally shows that current-time LMTPs admit a sequential regression representation of the identification formula, which greatly facilitates their estimation because it allows for the straightforward generalization of estimators developed for static interventions [16], [18]. In contrast, constructing a sequential regression representation for G-LMTPs requires more care. Consider a simple example with only two timepoints, no baseline covariates, and an intervention that only depends on the history of the natural value of treatment, i.e., \(A_1^d=d_1(A_1)\) and \(A_2^d = d_2(A_1, A_2(A_1^d))\). The outcome expectation under this intervention is identified as \[\int E[Y\mid A_2=d_2(a_2, a_1), L_2=l_2, A_1=d_1(a_1)]\mathrm{d}\mathsf{P}(a_2, l_2\mid A_1=d_1(a_1))\mathrm{d}\mathsf{P}(a_1).\] The integral with respect to \((a_2,l_2)\) cannot be written as a conditional expectation because the outcome regression involves setting \(A_1\) to two different values \(a_1\) and \(d_1(a_1)\). However, progress can be made by decoupling the value of \(a_1\) used in the conditioning set \(\mathsf{P}(a_2, l_2\mid A_1=d_1(a_1))\) from that used in \(d_2(a_2, a_1)\) (see §5 and Figure 2 for more discussion). Specifically, define \(m_2(a_2, h_2)=E[Y\mid A_2=a_2, L_2=l_2, A_1=a_1]\). Then, for each \(s_1\) in the support of \(A_1\), define \(m_1(s_1,a_1)=E[m_2(d_2(s_1, A_2), H_2)\mid A_1=a_1]\), and note that \(\theta=E[m_1(A_1, d_1(A_1))]\). This can be generalized as in the following proposition.

To illustrate this proposition, consider an example with \(\tau=3\) timepoints, no baseline covariates, and an intervention that only depends on the history of the natural value of treatment but not on covariates, i.e., \(A_1^d=d_1(A_1)\), \(A_2^d = d_2(A_1, A_2(A_1^d))\), and \(A_3^d = d_3(A_1, A_2(A_1^d), A_3(\bar A_2^d))\). According to Theorem 1, the outcome expectation under this intervention would be identified as \[\int E[Y\mid d_3(\bar a_3), l_3,d_2(\bar a_2), l_2, d_1(a_1)]\mathrm{d}\mathsf{P}(a_3, l_3\mid d_2(\bar a_2), l_2, d_1(a_1))\mathrm{d}\mathsf{P}(a_2, l_2\mid d_1(a_1))\mathrm{d}\mathsf{P}(a_1).\] The sequential regression algorithm in Proposition 1 would proceed as follows. The recursion is initialized as \(m_3(a_3, h_3) = E[Y\mid A_3=a_3, H_3=h_3]\) and \(q_3(\bar s_2, a_3, h_3) = m_3(d(\bar s_2, a_3), h_3)\), which is equal to \(E[Y\mid d_3(\bar s_2, a_3), l_3, a_2, l_2, a_1]\). Then, we start the iteration by computing \(m_2(\bar s_2, a_2, h_2) = E[q_3(\bar s_2, A_3, H_3)\mid A_2=a_2, H_2=h_2]\), which according to our definitions is equal to \[\int E[Y\mid d_3(\bar s_2, a_3), l_3,a_2, l_2, a_1]\mathrm{d}\mathsf{P}(a_3, l_3\mid a_2, l_2, a_1).\] The iteration then proceeds by computing \[\begin{align} q_2(s_1, a_2, h_2) &= m_2((s_1, a_2), d_2(s_1, a_2), h_2)\\ &=\int E[Y\mid d_3(s_1, a_2, a_3), l_3,d_2(s_1, a_2), l_2, a_1]\mathrm{d}\mathsf{P}(a_3, l_3\mid d_2(s_1, a_2), l_2, a_1) \end{align}\] and the regression \(m_1(s_1, a_1)=E[q_2(s_1, A_2, H_2)\mid A_1=a_1]\) which is equal to \[\int E[Y\mid d_3(s_1, a_2, a_3), l_3,d_2(s_1, a_2), l_2, a_1]\mathrm{d}\mathsf{P}(a_3, l_3\mid d_2(s_1, a_2), l_2, a_1)\mathrm{d}\mathsf{P}(a_2, l_2\mid a_1).\] The prediction \(q_1(A_1) = m_1(A_1, d(A_1))\) is averaged across values of \(A_1\) to obtain the identification formula. In other words, we augment the data using \(s_t\) to track the possible values that \(A_t\) can take in \(d_{t+1}\), while reserving the original notation \(a_t\) to compute integrals with respect to \(P(a_{t+1}, l_{t+1}\mid \ldots)\) (see Figure 2).

4 Optimality theory↩︎

Efficiency theory in this paper focuses on functions \(d\) that do not depend on \(\mathsf{P}\) (recall that the function is deterministic but allowed to take a randomizer as argument). For each \(t< k\), and each \(s_t\) define the random variable \(D_{t,k}(\bar s_k)=\mathbb{1}\{A_{t+1}=d_{t+1}(\bar s_{t+1}, H_{t+1}), \ldots A_k=d_k(\bar s_k, H_k)\}\). For \(t=\tau-1,\ldots,0\), define \(\bar\varphi_{t+1}(\bar s_t, Z) = \varphi_{t+1}(\bar s_t, Z) - m_t(\bar s_t, A_t, H_t)\), where \(m_0=\theta\) and \[\label{eq:eifgamma} \begin{align} \varphi_{t+1}(\bar s_t, Z)&= \sum_{k=t+1}^\tau\sum_{\bar s_{(t+1):k}}D_{t,k}(\bar s_k)\left(\prod_{u=t+1}^k \frac{g_u(s_u\mid H_u)}{g_u(A_u\mid H_u)}\right)\{q_{k+1}(\bar s_k, A_{k+1}, H_{k+1}) - m_k(\bar s_k, A_k, H_k)\}\\ &+ q_{t+1}(\bar s_t, A_{t+1}, H_{t+1}). \end{align}\tag{1}\] When necessary, we use the notation \(\varphi_t(\cdot;\eta)\) or \(\varphi_t(\cdot;\underline\eta_t)\) to highlight the dependence of \(\varphi_t\) (or \(\bar\varphi_t\)) on \(\underline\eta_t=(g_t,m_t,\ldots,g_\tau, m_\tau)\). We also use \(\eta\) to denote \((g_1,m_1,\ldots,g_\tau, m_\tau)\), and define \(\varphi_{\tau+1}(s_\tau, Z;\eta) = Y\). In the following, we let \(\eta'=(g_1',m_1',\ldots,g_\tau', m_\tau')\) denote some value of \(\eta\). This value will typically represent the probability limit of a given estimator \(\hat{\eta}\). The function \(\varphi_t\) satisfies the following property, which will be crucial to establish consistency of some estimators of \(\theta\) under multiple robustness assumptions, and to construct estimators of \(\theta\) under slow convergence rates for estimation of \(g_t\) and \(m_t\).

This lemma shares important connections to the von Mises-type expansions used in some of the semi-parametric inference literature [19]–[21]. Inspection of this lemma teaches us that an estimator of \(m_t\) constructed by regressing \(\varphi_{t+1}(\bar s_t, Z_i;\hat{\eta})\) for each \(s_t\), has a doubly robust point-wise error equal to \(\mathsf{Rem}_t(\bar s_t, a_t, h_t;\hat{\eta})\). Furthermore, this implies the following corollary:

If the modified treatment policy depends only on the contemporaneous natural value of treatment, i.e., if \(d_t(\bar a_t, h_t^d) = d_t(a_t, h_t^d)\), then all the above results reduce to the results presented by [5]. If the intervention only depends on some of the history of natural values of treatment, then the argument \(\bar s_{t-1}\) in the function \(q_t\) can be modified to only include values on which the intervention depends. We elaborate on this point in the Supplementary Materials §12.

5 Estimation and statistical inference↩︎

In this section we assume that preliminary estimators \(\hat{g}_t\), and \(\hat{m}_t\) are available. These estimators may be obtained from flexible regression methods from the machine learning literature or ensembles thereof [22]. As implied by Theorem 2, the consistency of these estimators will determine the consistency of our estimators of the parameter \(\theta\). We will use Super Learner [22], an ensemble regression algorithm known that builds a convex combination of regression algorithms in a user-given library, with weights chosen to minimize the cross-validated prediction error.

A naïve plug-in estimator for \(m_t:t=1,\ldots,\tau\) may be obtained using pooled regressions in an augmented dataset. Let \(K_t = |\bar{\mathcal{A}}_t|\). Define the augmented dataset \(\mathcal{D}_\tau = \left\{(i, Z_i, \bar s_\tau) : i \in \{1, \dots, n\},\; \bar s_\tau \in \bar{\mathcal{A}}_\tau \right\}\). That is, \(\mathcal{D}_\tau\) is a dataset with \(n K_\tau\) rows and \(\tau\) additional columns corresponding to possible exposure sequences \(\bar s_\tau\). Start by running a regression of \(Y_i\) on \((A_{\tau,i},H_{\tau,i})\) to obtain an estimator of \(m_\tau(a_\tau, h_\tau) = E[Y \mid A_\tau = a_\tau, H_\tau = h_\tau]\) among the original data. Generate predictions \(q_\tau(\bar s_{\tau-1}, A_{\tau}, H_{\tau}) = m_\tau(d_\tau((\bar s_{\tau-1}, A_\tau), H_\tau), H_\tau) = E[Y \mid A_\tau = d_\tau((\bar s_{\tau-1}, A_\tau), H_\tau), H_\tau]\) using the augmented data \(\mathcal{D}_\tau\). Because predictions are generated using the augmented data, \(q_\tau(\bar s_{\tau-1}, A_{\tau}, H_{\tau})\) is a \(n K_\tau\) length vector. Recursively for \(t=\tau-1,\ldots, 1\):

1. Add \(q_{t+1}(\bar s_{t}, A_{t+1}, H_{t+1})\) as a new column to \(\mathcal{D}_{t+1}\).

2. Construct \(\mathcal{D}_{t}\) by subsetting \(\mathcal{D}_{t+1}\) to unique combinations of \((i, \bar s_{t})\), such that \(\mathcal{D}_{t}\) has \(n K_t\) rows.

3. Obtain an estimator \(m_t(\bar{s}_t,a_t, h_t) = \mathsf{E}[q_{t+1}(\bar s_{t}, A_{t+1}, H_{t+1}) \mid \bar{s}_{t}, A_t = a_t, H_t = h_t]\) by regressing \(q_{t+1}(\bar s_{t}, A_{t+1}, H_{t+1})\) on \((\bar{s}_{t}, A_{t}, H_{t})\) among \(\mathcal{D}_t\).

4. Generate predictions \(q_t(\bar{s}_{t-1}, A_t, H_t) = m_t((\bar{s}_{t-1},A_t),d_t((\bar{s}_{t-1}, A_t), H_t), H_t).\)

An example of this approach is shown in Figure 2. Note that this approach exchanges \(\tau K_\tau\) regressions in a size-\(n\) dataset for \(\tau\) regressions in a sequence of \(n, n K_{\tau -1}, \dots, n K_1\) sized datasets which may be more computationally efficient. A similar pooling strategy has also previously been used for estimating so-called time-smooth counterfactual curves [23]. Here, we give the algorithm in full generality for interventions that may depend on the entire history \(\bar s_t\). When the intervention depends only on selected components of \(\bar s_t\), the augmented data need only carry those components. For example, delay policies require only a rolling lagged window of natural treatment values with length is determined by the magnitude of the delay. We give concrete examples in the Supplementary Materials §12.

None

Although this naïve estimator is simple to implement, its error is first order in the error of each regression estimate at time \(t\). This means that it is not possible to use this algorithm to construct an estimator that is generally \(\sqrt{n}\)-consistent under data-adaptive regression. To achieve this, we will leverage two approaches that rely on the canonical gradient of Proposition 1: targeted minimum loss-based estimation (TMLE), and sequential doubly robust (SDR) regression. Both of these estimators are \(\sqrt{n}\)-consistent, efficient, and asymptotically normal under assumptions, but the SDR has the advantage over TMLE that it has a simpler error structure therefore leading to potential gains in finite sample bias and consistency, whereas the advantage of the TMLE over SDR is that it is a substitution estimator therefore offering protection in cases where the parameter space is bounded and the true value of the parameter is close to the boundary. The discussion about the difference between TMLE and SDR estimators is given elsewhere [5], [16], [24] and is therefore omitted here.

5.1 Targeted minimum loss-based estimator↩︎

We use sample splitting to avoid imposing Donsker conditions on the regression estimators [25]–[27]. Let \({\cal V}_1, \ldots, {\cal V}_J\) denote a random partition of the index set \(\{1, \ldots, n\}\) into \(J\) prediction sets of approximately the same size. For each \(j\), the associated training sample is given by \({\cal T}_j = \{1, \ldots, n\} \setminus {\cal V}_j\). We let \(\hat{\eta}_{j}\) denote the estimator of \(\eta\) obtained by training the corresponding prediction algorithm using only data in the sample \({\cal T}_j\). Further, we let \(j(i)\) denote the index of the validation set which contains observation \(i\).

The targeted minimum loss-based estimator \(\hat{\theta}_{\scriptsize tmle}\) is computed as a substitution estimator that uses an estimate \(\tilde{m}_{1,j(i)}\) constructed to solve the cross-validated efficient influence function estimating equation \(\mathsf{P}_n\{\varphi_1(\cdot, \tilde{\eta}_{j(\cdot)}) - \hat{\theta}_{\scriptsize tmle}\}= 0\). The construction of \(\tilde{m}_{1,j(i)}\) is motivated by the observation that, for an appropriately defined variable \(\omega(\bar s_t,a_t, h_t)\) the efficient influence function of \(\theta\) can be expressed as a sum of terms of the form \(\omega(\bar s_t,a_t, h_t)\{q_{t+1}(\bar s_t, a_{t+1}, h_{t+1}) - m_t(\bar s_t, a_t, h_t)\}\), which, for fixed \(\bar s_t\), takes the form of score functions \(\omega(W)\{M - \mathsf{E}(M\mid W)\}\) for appropriately defined variables \(M\) and \(W\) and some weight function \(\omega\). If \(\mathsf{E}(M\mid W)\) is estimated within a weighted generalized linear model with canonical link that includes an intercept, then the weighted MLE estimate solves the score equation \(\sum_i\omega(W_i)\{M_i - \hat{\mathsf{E}}(M_i\mid W_i)\}=0\). TMLE uses this observation to iteratively tilt preliminary estimates of \(m_t\) towards a solution of the efficient influence function estimating equation. The algorithm is defined as follows:

1. Initialize \(\tilde{\eta} =\hat{\eta}\) and \(\tilde{m}_{\tau+1,j(i)}=\tilde{q}_{\tau+1,j(i)} = Y_i\).

2. For \(k=1,\ldots,\tau\), and for each \(\bar s_k\), compute the weights \(\omega_{k,i}(\bar s_k) = D_{0,k}(\bar s_k)\left(\prod_{u=1}^k \frac{g_u(s_u\mid H_u)}{g_u(A_u\mid H_u)}\right)\).

3. For \(t=\tau,\ldots,1\), for each \(\bar s_t\):

   • Fit the generalized linear tilting model \(\mathop{\mathrm{link}}\tilde{m}_t^\epsilon(\bar s_t, A_{t,i},H_{t,i}) = \epsilon + \mathop{\mathrm{link}} \tilde{m}_{t,j(i)}(\bar s_t, A_{t,1}, H_{t,i})\), where \(\mathop{\mathrm{link}}(\cdot)\) is the canonical link. The parameter \(\epsilon\) may be estimated by running a generalized linear model of the pseudo-outcome \(\tilde{q}_{t+1,j(i)}(\bar s_t, A_{t+1,i}, H_{t+1,i})\) with only an intercept term, an offset term equal to \(\mathop{\mathrm{link}}\tilde{m}_{t,j(i)}(\bar s_t, A_{t,i},H_{t,i})\), and weights \(\omega_{t,i}(\bar s_t)\), using all the data points in the sample. An outcome bounded in an interval \([a,b]\) may be analyzed with logistic regression (i.e., \(\mathop{\mathrm{link}}=\mathop{\mathrm{logit}}\)) by mapping it to an outcome \((0,1)\) through the transformation \((Y - a) / (b - a)(1-2\epsilon) + \epsilon\) for some small value \(\epsilon>0\). This approach has robustness advantages compared to fitting a linear model as it guarantees that the predictions in the next step remain within the outcome space see [28].

   • Let \(\hat{\epsilon}\) denote the maximum likelihood estimate, and update the estimates as \[\begin{align} \mathop{\mathrm{link}}\tilde{m}_{t, j(i)}^{\hat{\epsilon}}(\bar s_t, A_{t,i},H_{t,i})&= \hat{\epsilon} + \mathop{\mathrm{link}}\tilde{m}_{t, j(i)}(\bar s_t, A_{t,i},H_{t,i})\\ \mathop{\mathrm{link}}\tilde{q}_{t, j(i)}^{\hat{\epsilon}}(\bar s_{t-1}, A_{t,i}, H_{t,i})&= \hat{\epsilon} + \mathop{\mathrm{link}}\tilde{m}_{t, j(i)}((\bar s_{t-1}, A_{t,i}),d((\bar s_{t-1}, A_{t,i}), H_{t,i}), H_{t,i}). \end{align}\]

   • Update \(\tilde{m}_{t, j(i)}= \tilde{m}_{t, j(i)}^{\hat{\epsilon}}\), \(t = t-1\), and iterate.

4. The TMLE is defined as \(\hat{\theta}_{\scriptsize tmle}=\frac{1}{n}\sum_{i=1}^n\tilde{q}_{1,j(i)}(A_{1,i},H_{1,i})\).

The iterative procedure and the score equation argument above guarantee that the above estimates solve the efficient influence function estimating equation, which is crucial to prove that \(\sqrt{n}(\hat{\theta}_{\scriptsize tmle}-\theta)\) converges to a random variable \(N(0;\mathop{\mathrm{\mathsf{Var}}}[\bar\varphi_1(Z;\eta)])\), provided that \(\mathsf{Rem}_0(\hat{\eta})=o_P(n^{-1/2})\) and that the true and estimated weights \(\omega_{t}(s_i)\) are bounded. The proof of this result follows standard arguments in the analysis of TMLE presented elsewhere [5], [26] and is therefore omitted.

5.2 Sequential regression estimator using SDR unbiased transformations↩︎

In this section we use the multiply robust unbiased transformation in expression (3 ) to obtain an estimate of \(m_t\). We say that \(\varphi_{t+1}\) is a multiply robust unbiased transformation [16], [29]–[32] for \(m_t\) due to the following proposition, which is a straightforward consequence of Theorem 2.

The above motivates the construction of the sequential regression estimator by regressing an estimate of the data transformation \(\varphi_{t+1}(\bar s_t, Z;\eta)\) on \((A_t,H_t)\) for each \(\bar s_t\), starting at \(\varphi_{\tau+1}(Z;\eta)=Y\). For preliminary cross-fitted estimates \(\hat{g}_{1,j(i)},\ldots,\hat{g}_{\tau, j(i)}\), the estimator is defined as follows:

1. Initialize \(\varphi_{\tau+1}(\bar s_\tau, Z_i;\underline{\check\eta}_{\tau,j(i)})= Y_i\) for \(i=1,\ldots,n\).

2. For \(t=\tau,\ldots,1\) and each \(s_t\in \mathcal{A}_t\):

   • Compute the pseudo-outcome \(\check Y_{t+1,i} = \varphi_{t+1}(\bar s_t, Z_i;\underline{\check\eta}_{t,j(i)})\) for all \(i=1,\ldots,n\).

   • For \(j=1,\ldots,J\):

     • Regress \(\check Y_{t+1,i}\) on \((A_{t,i,}H_{t,i})\) using any regression technique and using only data points \(i\in \mathcal{T}_{j}\).

     • Let \(\check m_{t,j}\) denote the output, update \(\underline{\check\eta}_{t, j} = (\hat{g}_{t,j}, \check m_{t,j},\ldots,\hat{g}_{\tau,j}, \check m_{\tau,j})\), and iterate.

3. Define the sequential estimator regression as \(\hat{\theta}_{\scriptsize sr}=\frac{1}{n}\sum_{i=1}^n\varphi_1(Z_i;\check\eta_{j(i)})\).

Like its TMLE counterpart, this estimator satisfies \(\sqrt{n}(\hat{\theta}_{\scriptsize sr}-\theta)\) converges to a random variable \(N(0;\mathop{\mathrm{\mathsf{Var}}}[\bar\varphi_1(Z;\eta)])\), provided that \(\mathsf{Rem}_0(\hat{\eta})=o_P(n^{-1/2})\). This can be used to construct asymptotically valid confidence intervals and hypothesis tests, for example a \((1-\alpha)100\%\) confidence interval may be constructed as \(\hat{\theta}_{\scriptsize sr}\pm z_{1-\alpha/2}\hat{\sigma}/\sqrt{n}\), where \(\hat{\sigma}^2\) is the empirical variance of \(\bar\varphi_1(Z_i;\check\eta_{j(i)})\).

6 Numerical illustration↩︎

We use Monte Carlo simulation to illustrate the performance of the TMLE and SDR estimators. For \(n \in \left\{250, 500,1000,5000,10000 \right\}\), data were simulated from the data generating mechanism given in Figure 3.

None

This data generating mechanism corresponds to a survival analysis scenario with time-varying confounding and conditionally degenerate exposure assignment. At each sample size, we generated \(1000\) data-sets, and for each data-set, using both the SDR and TMLE, we estimated the expected value of \(Y_5\) under a hypothetical intervention where exposure (\(A_t\)) was delayed by one time-point. Outcome regressions were estimated using a generalized linear model with all possible two-way interactions; the propensity score at each time-point was estimated using a main-term generalized linear model with the logit link. Because the robustness results of the SDR and TMLE have been thoroughly discussed in the existing literature, (see e.g., [5], [8], [10]) we limit the simulation study to the ideal scenario where all nuisance parameters are correctly specified and estimated at fast enough rates. Performance was summarized using absolute bias, mean squared error (MSE), \(n \times \text{MSE}\), and nominal 95% confidence interval coverage. The true value under the intervention was calculated using Proposition 1 with correctly specified models and a “super-population” of one-million observations. The R code for the simulation can be found at https://github.com/CI-NYC/htlmtp-simulations.

The results of the simulation study are presented in Table 1. As expected, both the TMLE and SDR demonstrate improving metrics with increasing sample size. Interestingly, both estimators have greater than nominal confidence interval indicating conservative variance estimation at sample sizes greater than \(n = 1000\). Further, at \(n = 250\), the bias-variance tradeoff of the estimators appears slightly different with the SDR trading a slight increase in bias for smaller variance (as indicated by \(n \times \text{MSE}\)), while the TMLE was opposite. However, we caution against drawing comparative conclusions regarding finite sample performance from such a simulation study as estimator performance can vary greatly with the data-generating mechanism [33].

7 Motivating application↩︎

We applied our proposed approach to estimate the effect of delaying a risky pain treatment by one month on 12-month incidence of opioid use disorder (OUD) or overdose diagnosis among non-dual-eligble, Medicaid-insured patients, aged 19-63 years, with lumbar spinal stenosis, N=12745. Day 0 was the first day of a pain treatment claim, which occurred within 30 days of a new-onset diagnosis of lumbar spinal stenosis (ICD-10 code M4806). Patients must have been continuously enrolled in Medicaid for the 6 months prior to day 0 (a “washout period") to determine cohort eligibility, which, in addition to the above age and eligibility criteria, also included being opioid-naive, with no history of OUD, no cancer diagnosis, not pregnant, and community-dwelling. Additionally, they could not have been hospitalized in the 1 month prior to day 0. This cohort was drawn from Medicaid T-MSIS Analytic Files (TAF), Demographics, Other Services, Inpatient, and Pharmacy claims, 2016-2019 from the 25 states that implemented the ACA prior to 2014: AR, AZ, CA, CO, CT, DE, HI, IA, IL, KY, MA, MD, MI, MN, ND, NH, NJ, NM, NV, NY, OH, OR, RI, VT, WA, WV (excluded Maryland due to poor diagnostic code quality).

We were interested in examining this treatment delay question among patients with a lumbar spinal stenosis diagnosis, because these patients typically receive multiple months of pain treatment and frequently receive pain treatments where concerns have been raised about their OUD risks. Starting from day 0, we discretized follow-up into 1-month intervals, ending at month 12. During each month, we evaluated 1) which pain treatments the patient had a claim for, 2) whether or not the patient was censored, and 3) presence of an OUD or overdose diagnosis code or receipt of medication for OUD. We collected information on the following baseline covariates during the washout period: age; sex; race/ethnicity; receipt of: TANF benefits, SSI benefits; veteran status; household size; marital status; english as primary language; diagnosis codes for: anxiety, bipolar, depression, ADHD, other psychiatric illness, alcohol use disorder, or other substance use disorder; mental health counseling claims; healthcare utilization variables: number of inpatient hospitalizations in the washout period prior to one month before day 0, number of outpatient visits, number of emergency department visits. We considered the following set of pain treatment variables: i) opioid analgesic prescription of \(\le7\) days and dose \(\le50\) morphine milligram equivalents (MME), ii) opioid analgesic prescription of \(>7\) days and dose \(\le50\), iii) opioid analgesic prescription of dose \(>50\)MME, iv) gabapentoid prescription, v) benzodiazepine prescription, vi) muscle relaxant prescription, vii) duloxetine prescription, viii) anti-inflammatory prescription, ix) interventional procedure, x) chiropractic visit, xi) physical therapy visit. Patients could have, and frequently did have, more than one treatment during a given month. For the purpose of this illustrative analysis we dichotomized treatment into the presence vs absence of any of the following higher-risk treatments, where this binary treatment variable, \(A_t=1\) if any of treatments ii-v were present in month \(t\) and \(A_t=0\) otherwise. Patients were censored in month \(t\) if: they had disenrolled from Medicaid by the end of month \(t\), they became dual-eligible for Medicare, or 2019 ended.

The study was approved by the Columbia University Institutional Review Board. All code to replicate cohort and variable creation and statistical analyses is given in the GitHub repository: https://github.com/CI-NYC/m4806-diagnosis.

We were interested in the effect of delaying a higher-risk pain treatment by one month as compared to the observed average risk, assuming all censoring were prevented, \(\mathsf{E}[Y(\bar{A}^d, \bar{C}=\bar{1})- Y(\bar{C}=\bar{1})]\), where \(d_t(a_t, \bar a_{t-1}) =0\) if \(a_t = 1 \text{ and } a_s =0 \text{ for all } s < t\) and \(d_t(a_t, \bar a_{t-1})=a_t\) otherwise. We estimated the effect of this delay using both the TMLE and SDR estimators with 2-fold cross-fitting. The nuisance parameters \(g_t\) and \(m_t\) were estimated using the Super Learner [22] with a GLM and multivariate adaptive regression splines (MARS) as candidate algorithms. Results are shown in Table 2. Using both the TMLE and SDR to adjust for right-censoring, we estimated the 12-month probability of remaining free from an OUD or overdose diagnosis (survivial-type probability) as 0.963 (95% CI: 0.96, 0.966). Under a hypothetical one-month delay of risky pain treatments, the TMLE estimated the 12-month probability of remaining free from an OUD or overdose diagnosis as 0.985 (95% CI: 0.92, 1), where the upper bound was truncated at one. The SDR produced a point estimate of 1.27, outside the parameter space. This is a well-known limitation of the SDR when the true value lies near the boundary due to it not being a substitution estimator (§5), and we omit the confidence interval accordingly.

8 Conclusion↩︎

This paper proposes targeted learning estimators for longitudinal modified treatment policies whose hypothetical interventions may depend on the history of natural treatment values, extending previous work on modified treatment policies that depend only on the contemporaneous natural value of treatment. This extension allows the double/de-biased machine learning estimation of the effects of policies that modify the timing or trajectory of treatment decisions. The proposed augmented sequential regression representation, together with the corresponding efficient influence function, yields practical TMLE and SDR estimators for this broader class of interventions. Several directions of research remain open. We have focused on discrete exposures, for which the relevant natural treatment histories can be represented by summing over possible exposure sequences. Many applications motivating modified treatment policies involve continuous, multivariate, or mixed discrete-continuous treatments, including medication doses and combinations of therapies. Future work should develop efficient estimators for G-LMTPs in these settings, avoiding explicit enumeration of treatment histories and replacing conditional mass-function ratios with appropriate density-ratio or Riesz representers. Other useful directions include scalable algorithms for long follow-up, diagnostics for support violations induced by history-dependent policies, and extensions to settings with censoring, competing risks, interference, or continuous time.

Acknowledgments↩︎

The authors are grateful to Alex Luedtke for insightful conversations that helped clarify the motivation for the generalization developed in this paper, as well as Anton Hung for the creation of the dataset used in our motivating application.

Iván Díaz was supported by a Patient-Centered Outcomes Research Institute (PCORI) Project Program Award (ME-2024C3-43044). Paweł Morzywołek was partially supported through a National Institutes of Health (NIH) Award (5UM1AI068635-20), the Pacific Institute for the Mathematical Sciences (PIMS), UW eScience Institute and the Pioneer Center for SMARTbiomed. Kara Rudolph was supported by a National Institute on Drug Abuse (NIDA) Award (R01DA053243).

Supplementary Materials for Modified treatment policies that depend on the natural history of treatment.↩︎

9 Identification (Theorem 1)↩︎

Recall that for fixed values \(\bar a_t\), \(\bar h_t\), we denote \(a_t^d=d(\bar a_t, h^d_t)\), where \(h^d_t=(\bar a_{t-1}^d, \bar l_t)\). Define the random variable \[Z_t(\bar a_t)= Y\{\bar a_t^d, \underline d_{t+1}(\bar a_t, \underline A_{t+1}(\bar a_t, \underline A_{t+1}^d), H_{t+1}(\bar a_t, \underline A_{t+1}^d))\}.\] Note that \(Z_0(a_0)=Y(\bar A^d)\), and that \(Z_t(\bar a_t)=Z_{t+1}(\bar a_{t+1})\) in the event \(A_{t+1}=a_{t+1}\) and \(H_{t+1} = h_{t+1}^d\), and that \(Z_\tau(\bar a_\tau) = Y\) in the event \(A_\tau=a_\tau\) and \(H_\tau = h_\tau^d\). We have \[\begin{align} E[Y(\bar A^d)] & =\int_{\mathcal{A}_1, \mathcal{L}_1} \mathsf{E}[Z_0(a_0)\mid A_1=a_1, L_1=l_1]\mathrm{d} \mathsf{P}(a_1,l_1)\notag\\ & =\int_{\mathcal{A}_1, \mathcal{L}_1} \mathsf{E}[Z_1(a_1)\mid A_1=a_1, L_1=l_1]\mathrm{d} \mathsf{P}(a_1,l_1)\notag\\ & =\int_{\mathcal{A}_1, \mathcal{L}_1} \mathsf{E}[Z_1(a_1)\mid A_1=a_1^d, L_1=l_1]\mathrm{d} \mathsf{P}(a_1,l_1)\notag\\ & =\int_{\mathcal{\bar A}_2, \mathcal{\bar L}_2} \mathsf{E}[Z_1(a_1)\mid A_2=a_2, L_2=l_2,A_1=a_1^d,L_1=l_1]\mathrm{d}\mathsf{P}(a_2,l_2\mid a_1^d, l_1)\mathrm{d}\mathsf{P}(a_1,l_1)\notag\\ & =\int_{\mathcal{\bar A}_2, \mathcal{\bar L}_2} \mathsf{E}[Z_2(\bar a_2)\mid A_2=a_2^d, L_2=l_2,A_1=a_1^d,L_1=l_1]\mathrm{d}\mathsf{P}(a_2,l_2\mid a_1^d, l_1)\mathrm{d}\mathsf{P}(a_1,l_1)\notag\\ & =\int_{\mathcal{\bar A}_2, \mathcal{\bar L}_2} \mathsf{E}[Z_2(\bar a_2)\mid A_2=a_2^d, H_2=h_2^d]\mathrm{d}\mathsf{P}(a_2,l_2\mid a_1^d, l_1)\mathrm{d}\mathsf{P}(a_1,l_1)\notag\\ &\,\,\, \vdots\notag\\ & =\int_{\mathcal{\bar A}_\tau, \mathcal{\bar L}_\tau} \mathsf{E}[Z_\tau(\bar a_{\tau})\mid A_\tau=a_\tau^d, H_\tau=h_\tau^d] \prod_{k=1}^\tau\mathrm{d}\mathsf{P}(a_k,l_k\mid a_{k-1}^d, h_{k-1}^d)\notag\\ & =\int_{\mathcal{\bar A}_\tau, \mathcal{\bar L}_\tau} \mathsf{E}[Y\mid A_\tau=a_\tau^d, H_\tau=h_\tau^d] \prod_{k=1}^\tau\mathrm{d}\mathsf{P}(a_k,l_k\mid a_{k-1}^d, h_{k-1}^d)\label{eq:eq1} \end{align}\tag{4}\] Where the first and fourth equalities follow by law of iterated expectation. 0◻

This lemma follows trivially after noticing that under the assumed NPSEM, \(Z_t(\bar a_{t-1})\) is a deterministic function of \((\underline U_{L,t+1}, \underline U_{A,t+1})\). 0◻

10 Efficient influence functions (Corollary 1)↩︎

In this proof we will use \(\Theta(\mathsf{P})\) to denote a parameter as a functional that maps the distribution \(\mathsf{P}\) in the model to a real number. A function \(\varphi_1(Z;\eta)\) is the EIF of a parameter functional \(\Theta(P)\) if it satisfies \[\frac{\mathrm{d}}{\mathrm{d}\epsilon}\Theta(P_\epsilon)\bigg|_{\epsilon=0} = \mathsf{E}[\varphi_1(Z;\eta) s(Z)],\label{eq:defeif}\tag{5}\] where \(\mathsf{P}_\epsilon\) is a smooth parametric submodel that locally covers the non-parametric model, with score \[s(Z)=\left(\frac{\mathrm{d}\log \mathsf{P}_\epsilon}{\mathrm{d} \epsilon}\right)\bigg|_{\epsilon = 0}\] such that \(\mathsf{P}_{\epsilon=0}=\mathsf{P}\). We start by conjecturing an EIF, and then prove it is in fact the EIF using the above definition.

To conjecture an EIF we find the EIF in a model where we assume that the measure \(\nu\) is discrete so that integrals can be written as sums. Recall that for fixed values \(\bar a_t\), \(\bar h_t\), we denote \(a_t^d=d(\bar a_t, h^d_t)\), where \(h^d_t=(\bar a_{t-1}^d, \bar l_t)\). We can express the non-parametric MLE of \(\theta\) as \[\Theta(\mathsf{P}_n)=\sum_{\bar a_\tau, \bar l_{\tau+1}}l_{\tau+1}\mathsf{P}_n f_{l_{\tau+1},a_\tau^d,h_\tau^d}\prod_{k=1}^\tau\frac{\mathsf{P}_nf_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}_nf_{a_k^d, h_k^d}}\label{nonpest},\tag{6}\] where we remind the reader of the notation \(\mathsf{P}f =\int f \mathrm{d}\mathsf{P}\). Here \(f_{x}(X)=\mathbb{1}(X=y)\) for any variable \(X\) and value \(x\). In the above display we used the assumption that \(d\) does not depend on \(\mathsf{P}\), and therefore does not need to be estimated.

We will use the fact that the efficient influence function in a non-parametric model corresponds with the influence curve of the NPMLE. This is true because the influence curve of any regular estimator is also a gradient, and a non-parametric model has only one gradient. The Delta method shows that if \(\hat{\Theta}(\mathsf{P}_n)\) is a substitution estimator such that \(\theta=\hat{\Theta}(\mathsf{P})\), and \(\hat{\Theta}(\mathsf{P}_n)\) can be written as \(\hat{\Theta}^*(\mathsf{P}_nf:f\in\mathcal{F})\) for some class of functions \(\mathcal{F}\) and some mapping \(\Theta^*\), then the influence function of \(\hat{\Theta}(\mathsf{P}_n)\) is equal to \[\mathop{\mathrm{IF}}_\mathsf{P}(Z)=\sum_{f\in\mathcal{F}}\frac{\mathrm{d}\hat{\Theta}^*}{\mathrm{d}\mathsf{P}f}\{f(O)-\mathsf{P}f\}.\]

Applying this result to (6 ) with \(\mathcal{F}=\{f_{l_{\tau+1},a_\tau^d,h_\tau^d}, f_{a_s, l_s, a_{s-1}^d, h_{s-1}^d}, f_{a_s^d, h_s^d}: l_{\tau+1}, a_\tau, h_\tau; s=1,\ldots,\tau\}\) and rearranging terms gives \(\mathop{\mathrm{IF}}_P(Z) = \varphi_1(Z;\eta)\). In particular we have \[\begin{align} \frac{\mathrm{d}\Theta(\mathsf{P})}{\mathrm{d}\mathsf{P}f_{l_{\tau+1},a_\tau^d,h_\tau^d}} &= l_{\tau+1}\prod_{k=1}^\tau\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}, \end{align}\] and \[\begin{align} &\frac{\mathrm{d}\Theta(\mathsf{P})}{\mathrm{d}\mathsf{P}f_{a_s, l_s, a_{s-1}^d, h_{s-1}^d}} = \sum_{l_{\tau+1}, \underline a_{s+1}, \underline h_{s+1}}l_{\tau+1}\frac{\mathsf{P} f_{l_{\tau+1},a_\tau^d,h_\tau^d}}{\mathsf{P}f_{a_s, l_s, a_{s-1}^d, h_{s-1}^d}}\prod_{k=1}^\tau\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\\ &= \left(\sum_{\underline a_{s+1}, \underline h_{s+1}}m_\tau(a_\tau^d, h_\tau^d)\frac{\mathsf{P} f_{a_\tau^d,h_\tau^d}}{\mathsf{P} f_{a_s^d,h_s^d}}\prod_{k=s+1}^\tau\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\right)\prod_{k=1}^{s-1}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\\ &= \left(\sum_{\substack{a_{s+1},\ldots,a_{\tau-1}\\ l_{s+1}, \ldots, l_{\tau-1}}}\left[\sum_{a_\tau, l_\tau}m_\tau(a_\tau^d, h_\tau^d)\frac{\mathsf{P}f_{a_\tau, l_\tau, a_{\tau-1}^d, h_{\tau-1}^d}}{\mathsf{P}f_{a_{\tau-1}^d, h_{\tau-1}^d}}\right]\frac{\mathsf{P} f_{a_{\tau-1}^d,h_{\tau-1}^d}}{\mathsf{P} f_{a_s^d,h_s^d}}\prod_{k=s+1}^{\tau-1}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\right)\prod_{k=1}^{s-1}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\\ &= \left(\sum_{\substack{a_{s+1},\ldots,a_{\tau-1}\\ l_{s+1}, \ldots, l_{\tau-1}}}\left[m_{\tau-1}(\bar a_{\tau-1}, a_{\tau-1}^d, h_{\tau-1}^d)\right]\frac{\mathsf{P} f_{a_{\tau-1}^d,h_{\tau-1}^d}}{\mathsf{P} f_{a_s^d,h_s^d}}\prod_{k=s+1}^{\tau-1}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\right)\prod_{k=1}^{s-1}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\\ &= m_s(\bar a_s, a_s^d, h_s^d)\prod_{k=1}^{s-1}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}, \end{align}\] as well as \[\begin{align} \frac{\mathrm{d}\Theta(\mathsf{P})}{\mathrm{d}\mathsf{P}f_{a_s^d,h_s^d}} &= -\sum_{l_{\tau+1}, \underline a_{s+1}, \underline h_{s+1}}l_{\tau+1}\frac{\mathsf{P} f_{l_{\tau+1}, a_\tau^d, h_\tau^d}}{\mathsf{P}f_{a_s^d, h_s^d}}\prod_{k=1}^\tau\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\\ &=-\sum_{l_{\tau+1}, \underline a_{s+1}, \underline h_{s+1}}l_{\tau+1}\frac{\mathsf{P} f_{l_{\tau+1}, a_\tau^d, h_\tau^d}}{\mathsf{P}f_{a_s^d, h_s^d}}\prod_{k=s+1}^\tau\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\prod_{k=1}^s\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}\\ &= -m_s(\bar a_s, a_s^d, h_s^d)\prod_{k=1}^{s}\frac{\mathsf{P}f_{a_k, l_k, a_{k-1}^d, h_{k-1}^d}}{\mathsf{P}f_{a_k^d, h_k^d}}. \end{align}\] Notice now that \[\begin{align} \sum_{l_{\tau+1},a_\tau,h_\tau}&\frac{\mathrm{d}\Theta(\mathsf{P})}{\mathrm{d} f_{l_{\tau+1},a_\tau^d,h_\tau^d}} f_{l_{\tau+1},a_\tau^d,h_\tau^d}(Z) \\ &=\sum_{a_\tau,h_\tau}\mathbb{1}\{A_\tau = a_\tau^d, H_\tau = h_\tau^d\}L_{\tau+1}\prod_{k=1}^\tau\frac{\mathsf{P}(A=a_k,L_k= l_k, A_{k-1}=a_{k-1}^d, H_{k-1}=h_{k-1}^d)}{\mathsf{P}(A_k=a_k^d, H_k=h_k^d)}\\ &=L_{\tau+1}\sum_{\bar a_\tau}\mathbb{1}\{A_\tau = d(\bar a_\tau, H_\tau),\ldots, A_1=d(a_1, L_1)\}\prod_{k=1}^\tau\frac{g_k(a_k \mid H_k)}{g_k(A_k\mid H_k)}\\ &=L_{\tau+1}\frac{\bar g_\tau^d(A_\tau, H_\tau)}{\bar g_\tau(A_\tau, H_\tau)}, \end{align}\] as well as \[\begin{align} &\sum_{a_s,h_s}\frac{\mathrm{d}\Theta(\mathsf{P})}{\mathrm{d}\mathsf{P}f_{a_s,l_s, a_{s-1}^d, h_{s-1}^d}} f_{a_s, l_s, a_{s-1}^d, h_{s-1}^d}(Z) \\ &=\sum_{\bar a_{s-1}}\mathbb{1}\{A_{s-1} = d(\bar a_{s-1}, H_{s-1}),\ldots, A_1=d(a_1, L_1)\}m_s((\bar a_{s-1}, A_s), d((\bar a_{s-1}, A_s), H_s), H_s)\prod_{k=1}^{s-1}\frac{g_k(a_k \mid H_k)}{g_k(A_k\mid H_k)}, \end{align}\]and \[\begin{align} \sum_{a_s,h_s}&\frac{\mathrm{d}\Theta(\mathsf{P})}{\mathrm{d}\mathsf{P}f_{a_s^d, h_s^d}} f_{a_s^d, h_s^d}(Z) \\ &=-\sum_{\bar a_s}\mathbb{1}\{A_s = d(\bar a_s, H_s),\ldots, A_1=d(a_1, L_1)\}m_s(\bar a_s, A_s, H_s)\prod_{k=1}^s\frac{g_k(a_k \mid H_k)}{g_k(A_k\mid H_k)}, \end{align}\]Putting these results together and summing over \(s\) yields the conjectured \(\mathop{\mathrm{IF}}_P(Z)=\varphi_1(Z;\eta)-\Theta(P)\). Then, we can confirm that \(\varphi_1(Z;\eta)\) satisfies (5 ). To see this, note that Lemma 2 implies that the functional satisfies the following expansion: \[\Theta(\mathsf{P}_\epsilon) = \Theta(\mathsf{P}) + \int \{\varphi_1(Z;\eta) - \Theta(\mathsf{P})\} \mathrm{d}\mathsf{P}_\epsilon - \mathsf{Rem}_0(\eta_\epsilon,\eta),\] Differentiating with respect to \(\epsilon\) and evaluating at \(\epsilon=0\) yields \[\begin{align} \frac{\mathrm{d}}{\mathrm{d}\epsilon}\Theta(\mathsf{P}_\epsilon)\bigg|_{\epsilon = 0} &= \int\{\varphi_1(Z;\eta) - \Theta(\mathsf{P})\} \left(\frac{\mathrm{d}\mathsf{P}_\epsilon}{\mathrm{d}\epsilon}\right)\bigg|_{\epsilon = 0} - \frac{\mathrm{d}}{\mathrm{d} \epsilon}\mathsf{Rem}_0(\eta_\epsilon,\eta)\bigg|_{\epsilon = 0} \\ &= \int\{\varphi_1(Z;\eta) - \Theta(\mathsf{P})\} \left(\frac{\mathrm{d}\log \mathsf{P}_\epsilon}{\mathrm{d} \epsilon}\right)\bigg|_{\epsilon = 0}\mathrm{d}P - \frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathsf{Rem}_0(\eta_\epsilon,\eta)\bigg|_{\epsilon = 0}, \end{align}\] and the result follows after noticing that \[\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathsf{Rem}_0(\eta_\epsilon,\eta)\bigg|_{\epsilon = 0}=0.\] 0◻

11 Sequential double robustness (Lemma 2)↩︎

We first prove an auxiliary result.

The right hand side can be written as \[\begin{align} \sum_{a_k, l_k}&\sum_{s_k}\mathbb{1}\{a_k=d(\bar s_k, h_k)\}\frac{g_k(s_k\mid h_k)}{g_k(a_k\mid h_k)} m_{k,1}(\bar s_{k}, a_{k}, h_{k})g_k(a_k\mid h_k)p(l_k\mid a_{k-1}, h_{k-1})\\ &=\sum_{a_{k}, l_{k}}\sum_{s_{k}}\mathbb{1}\{a_{k}=d(\bar s_{k}, h_{k})\}g_{k}(s_{k}\mid h_{k}) m_{k,1}(\bar s_{k}, a_{k}, h_{k})p(l_{k}\mid a_{k-1}, h_{k-1})\\ &=\sum_{l_{k}, s_{k}}g_{k}(s_{k}\mid h_{k}) m_{k,1}(\bar s_k, d(\bar s_{k}, h_{k}), h_{k})p(l_{k}\mid a_{k-1}, h_{k-1})\\ &=\sum_{l_{k}, a_{k}}g_{k}(a_{k}\mid h_{k}) q_{k,1}(\bar s_{k-1}, a_{k}, h_{k})p(l_{k}\mid a_{k-1}, h_{k-1})\\ &= \mathsf{E}[q_{k,1}(\bar s_{k-1}, A_k, H_k)\mid A_{k-1}=a_{k-1}, H_{k-1}=h_{k-1}]. \end{align}\] 0◻
We now proceed with the proof of the main result (Lemma 2).

This lemma follows from recursive application for \(k=t+1,\ldots,\tau\) of the following relation: \[\begin{align} &\textcolor{blue}{\mathsf{E}[q_{k,1}}\textcolor{blue}{(\bar s_{k-1}, A_k, H_k) - q_k(\bar s_{k-1}, A_k, H_k)\mid A_{k-1}, H_{k-1}]} =\\ &\mathsf{E}\left\{\sum_{s_k}D_k(\bar s_k)\left[\frac{\mathsf{g}_k(s_k\mid H_k)}{\mathsf{g}_k(A_k\mid H_k)}-\frac{\mathsf{g}_{k,1}(s_k\mid H_k)}{\mathsf{g}_{k,1}(A_k\mid H_k)}\right][m_{k,1}(\bar s_{k}, A_k, H_k) - m_k(\bar s_{k}, A_k, H_k)]\,\,\bigg|\,\,A_{k-1},H_{k-1}\right\}\\ -&\mathsf{E}\left\{\sum_{s_k}D_k(\bar s_k)\frac{\mathsf{g}_{k,1}(s_k\mid H_k)}{\mathsf{g}_{k,1}(A_k\mid H_k)}[q_{k+1, 1}(\bar s_k, A_{k+1}, H_{k+1}) - m_{k, 1}(\bar s_k, A_k, H_k)]\,\,\bigg|\,\,A_{k-1},H_{k-1}\right\}\\ +&\mathsf{E}\left\{\sum_{s_k}D_k(\bar s_k)\frac{\mathsf{g}_{k,1}(s_k\mid H_k)}{\mathsf{g}_{k,1}(A_k\mid H_k)}\textcolor{blue}{\mathsf{E}[q_{k+1, 1}(\bar s_k, A_{k+1}, H_{k+1}) - q_{k+1}(\bar s_k, A_{k+1}, H_{k+1})\mid A_k, H_k]}\,\,\bigg|\,\,A_{k-1},H_{k-1}\right\}, \end{align}\] which follows from rewriting it as \[\begin{align} a-b &= \sum (c-d)(a_1 - b_1)\\ &- \sum d(a_2-b_2)\\ &+ \sum d(a_3-b_3), \end{align}\] after noticing that \(\sum d a_2=\sum da_3\), \(\sum d b_2=\sum da_1\), \(\sum d b_1=\sum db_3\), \(a=\sum c a_1\), and \(b=\sum c b_1\), where the latter two follow from the above proposition. 0◻

12 Delay intervention with three time points↩︎

Proposition 1 shows that the augmented history \(\bar s_t\) need not be carried in full when the intervention depends only on selected components of the natural treatment history. Delay interventions provide a simple illustration. For an absorbing binary treatment status, a delay by \(\kappa\) time points may be written as \[d_t^{(\kappa)}(\bar a_t,h_t^d)= \begin{cases} 0, & t\leq \kappa,\\ a_{t-\kappa}, & t>\kappa. \end{cases}\] Thus, the rule at time \(t\) does not generally require the entire history \(\bar a_t\), only the corresponding value of treatment \(a_{t-\kappa}\). Equivalently, the sequential regression implementation only needs to augment the data by the relevant history. We discuss two examples.

None

12.1 Time horizon considerations.↩︎

When natural treatment initiation occurs close to the end of follow-up, the delayed treatment initiation time may fall outside the study horizon. For example, in a vaccine study with \(\tau=10\) monthly decision times and a five-month delay, individuals naturally vaccinated in months \(6,\ldots,10\) would remain unvaccinated throughout the observed ten-month horizon under the delayed intervention. This does not create a problem with our technical framework, but it may affect interpretation because, for late vaccines, the delayed intervention is equivalent within the study horizon to no vaccination.

References↩︎