1.1 Illustrative example

Throughout this paper, we will illustrate our presentation using an example from mental health research. The STAR*D (Sequenced Treatment Alternatives to Relieve Depression) trial is a multi-level, longitudinal pragmatic trial of treatment strategies for major depression. After an initial screening, patients are enrolled into level 1 of the study, where everyone is assigned citalopram. At the start of each subsequent level, if a patient still remains in the study, then he/she is randomized to an option within one of the chosen treatment strategies. At level 2, the patient and his/her provider can choose between augmenting the citalopram with multiple options or switching to a new regimen with multiple options. We are interested in the comparative effect of augmenting vs switching medication. Because these two strategies are not randomized, this analysis is analogous to an observational study.

Suppose we wish to assess the effect modification of level 2 treatment strategy (augmenting vs switching) by the depression symptoms measured prior to entering level 2. Covariates either measured at enrollment or level 1 exit, for our purpose, are both considered baseline modifiers. These symptom measures are obtained at clinical visits and the level 1 exit survey. It is reasonable to believe that the more depressed patients and those less satisfied with their level 1 treatment may be less inclined to follow up with clinical visits or to respond to surveys. In this case, a simple complete-case analysis of effect modification may only provide results applicable to patients with less severe symptoms or more satisfied with their level 1 treatment. If we assume that this missingness, while associated with the outcome, can be predicted using covariate history collected up to level 1 exit, then we can regain some of the generalizability. Indeed, we cast these symptom measures as counterfactual variables under an intervention on their missingness mechanism. This way, the target parameter is formulated in terms of an ideal experiment where the report of the symptom measures was always ensured (i.e., intervened to be non-missing).

1.2 Organization of the article

This paper is organized as follows. In Section 2, we formulate the causal inference problem and determine the identifiability of the desired causal parameters from the observed data distribution. In Section 3, we present the EIF for the parameter of interest under a saturated semiparametric model, as well as a projected influence function. The robustness conditions for the efficient and projected influence functions are also established. We then present the construction of two TMLE estimators, one using the EIF (we call it the full TMLE) and one using the projected influence function (we call it the projected TMLE). For comparison we presented an IPTW estimator and a non-targeted substitution estimator. In Section 4, a simulation study demonstrates robustness properties of the projected TMLE. In Section 5, we use the STAR*D trial to illustrate the application of the projected TMLE. A final summary concludes the article.

2.1 Data structure

Specifically, consider a longitudinal data structure

where W encodes baseline covariates, A_t encodes the exposures of interest and censoring indicators at time t,V encodes the time-varying history (between treatment A0 and A1) on which one wishes to condition, Lt encodes covariates (including time-varying confounders and the outcome process Yt) measured between At and At+1. Variable YK is the final outcome of interest. For the sake of discussion, assume that YK is either a binary or a bounded continuous variable (without losing generality, we may assume it is bounded in (0,1)). We note that in general, A_t encodes intervention variables of interest, i.e. variables whose distribution one would set in an ideal experiment. These include exposure variables (where one would randomize in an ideal experiment) as well as censoring variables (where one would prevent in an ideal experiment).

For later convenience, we introduce some useful notations. For a time-varying variable Xs, let Xs≡X1,…,Xs and Xs,t≡Xs,…,Xt for s<t. We will also use the shorthand X=XK. The equivalent notations in lower case will apply to the realizations of the corresponding variable. Finally, variables with degenerate indices, such as Lt at t=−1, are empty.

The time-ordering assumptions can be captured by a nonparametric structural equations model (NPSEM, [12]):

This framework assumes that each variable X in the observed data structure is an unknown deterministic function of all its preceding observed variables and some unmeasured exogenous random factors U. We will refer to the observed variables in the input of fX as the parents of X and denote their set as Pa(X).

Concretely, we describe the data example introduced in Section 1.1 with these notations. Our goal is to assess the effect modification of level 2 treatment strategy (augmenting vs switching) by the depression symptoms measured prior to entering level 2; yet many of these baseline effect modifiers were missing at random. By study protocol, all subjects will have either entered remission (treatment success), moved onto the next level (treatment failure), or dropped out (right censoring), by the end of 23 weeks of the level 2 period. For our study goal, the variables that we would enforce/randomize in an ideal experiment are the missingness status of the baseline effect modifier, the treatment assignment at the start of level 2, and the censoring status in each of the 23 weeks.

To this end, let V be the baseline effect modifier of interest. Let A0 be the indicator for non-missingness of V, and A1 be the treatment strategy received by a patient at level 2 (A1=0 for augmenting medication and A1=1 for switching medication). We use W to encode the baseline covariates that may affect the measurement status A0, and L0 to encode other baseline covariates (beside the modifier of interest) measured prior to treatment assignment A1. Hence, under this indexing, K=24 and the first week of level 2 starts at t=2. For t≥2,A_t is a counting process which drops to 0 if patient was censored by time t, and Lt are time-varying covariates such as visit statistics (duration in level thus far, visit frequency, etc), side-effect burden and symptom measures at time t. The variable Lt also contains two counting processes: the outcome process Yt, which is a binary indicator for entering remission by time t, and an exit process Et that jumps to 1 if a patient is moved to the next level, in which case the remission status will be zero for this level and the patient is considered non-censored (since the outcome was observed to be unsuccessful). Our final outcome of interest is YK — the remission status by the end of 23 weeks.

2.2 Parameter of interest

To assess the effect of the interventions, we study the intervention-specific counterfactual mean outcomes as a function of the interventions. These counterfactual mean outcomes can be obtained from an ideal experiment where one intervenes to assign A=a, and measures the resulting covariates and the outcome of interest. In our example, an ideal experiment would set A0=1 (always measure baseline modifier), A1=a1 (switching or augmenting), and A2,K=1 (always prevent drop outs). This data-generating process can be formalized using eq. (1) by setting At=at in the equations for A_t and in the parents of the non-intervention variables V,Lt and Y. We denote these resulting, possibly counterfactual, covariates as V(a0),Lt(at) and YK(a).

Now, we wish to assess effect modification of such treatment effects by changing values of V. That is, we ask “how does the differential effect of (a0,a∗1,K) vs (a0,a1,K), and of (a0,a1,K) vs (a0∗,a1,K), differ as a function of V, which is affected by the treatment assignment A₀?”. To answer this question, our parameters of interest are the so-called Counterfactual-History-Adjusted mean outcomes [1]

and our goal is to study how these mean outcomes change as a function of a0,a1,K and v. As discussed in Section 1, eq. (2) differs from the traditional History-Adjusted mean outcomes, EY(A0,a1,K)∣V=v,A0=a0, in that these traditional parameters condition on the observed treatment and modifier values. Hence within each such stratum, the conditional mean outcome may still depend on the observed treatment mechanism. As such, these parameters are affected by potential selection bias.

In our example, the conditional mean outcomes we wish to study are E(Y(1,a1,K)∣V(1)=v), with {a1,K=(a1,…,aK):a1∈{0,1};at=1,t≥2}. In words, it is the intervention-specific counterfactual mean outcome in an ideal experiment without missing measurements or censoring.

A challenge to study of eq. (2) is that the curse of dimensionality arises in applications with more than two time points and/or with categorical or continuous V. To address this issue, it is useful to summarize eq. (2) using a working MSM, mψ(a,v), which is parametrized by a finite dimensional ψ∈S⊂Rd. Since the outcome Y is binary (or bounded within the unit interval), concretely we consider a logistic MSM

where ϕ(a,v) is the vector of linear predictors in the generalized linear model. The methods presented here are easily modified to other forms of MSM.

Given a CHA-MSM eq. (3) for the conditional mean outcomes eq. (2), the true MSM parameter ψ can be interpreted as the best summary measure of ρ(a,v) as a function of a and v. Formally, let A denote the set of interventions of interest and V denote the set of modifier values of interest. For a given kernel weight function h(a,v), the true MSM parameter ψ in eq. (3) is defined as

In words, ψ yields the best weighted approximation of the counterfactual conditional dose-response curve ρ(a,v), according to the quasi-loglikelihood loss, kernel weights and working model mψ(a,v).

The rest of this paper is devoted to the identification and inference of eq. (4) from the observed data.

2.3 Statistical estimand

To identify eq. (4) from the data generating distribution P0, we make the Positivity Assumption (PA) and the Sequential Randomization Assumption (SRA, derived by [13]). Specifically, under the PA, there exists αt>0 such that αt≤P0(At=at∣Pa(At)), for all t and a∈A, almost everywhere. The SRA assumes that At⊥V(a0),{Lt(a):t}∣Pa(At). Under these conditions, the joint distribution V(a0),{Lt(a):t} is identifiable from the observed data distribution P0. In our STAR*D example, the plausibility of the SRA can be fortified by measuring sufficiently many confounders of the modifier’s missingness, the treatment selection, and the censoring mechanism.

By calculations shown in the Appendix, the SRA allows us to identify P(V(a0)=v) as

and the counterfactual mean outcome ρ(a,v) as

where, for t=0,…,K,

where lt,j−1=(lt,...,lj−1) for t<j-1. Equation (7) represents a conditional expectation of the final outcome YK given observed histories (Lt−1,V,W) and under an intervention a on A.

For simplicity, we suppressed the notation for the functionals γ,ρ,Qta, as each should be evaluated with respect to a given distribution. We will use γ(P),ρ(P),Qta(P) to express explicitly these functionals evaluated at a specific distribution P when needed, but will keep the shorthand notation otherwise.

In comparison to eq. (6), parameters in conventional settings where the modifier V is not counterfactual can be written as

The weight 1/P(A0=a0∣W)E1/P(A0=a0∣W)∣V=v,A0=a0 adjusts for potential selection bias introduced by treatment assignment A0=a0. Indeed, if A0 does not depend on W (in our STAR*D example, this means missingness is completely at random), then this weight equals 1, in which case ρ(a,v) is equivalent to the estimand in an analysis stratified by A0.

Combining eqs. (5) and (6), the causal MSM parameter ψ in eq. (4) identifies to

where the functionals γ,ρ are evaluated at P0. In the forthcoming sections, we study the statistical inference of Ψ(P0).

2.4 Notations

Before we proceed, let us introduce some more useful definitions and notations. Let M be a saturated semiparametric model containing our data generating distribution P0. The parameter of interest in eq. (8) is the map P↦Ψ(P), from M to Rd, evaluated at P₀.

Suppose we observe n i.i.d. copies of O∼P0. Let Pn denote the empirical distribution of this sample. For a function f of O, we will write Pnf≡1n∑i=1nf, and for a distribution P, we will write Pf=EPf(O).

Bang and Robins [14] noted the following recursive property when dealing with longitudinal intervention-specific mean outcomes:

for t=0,…,K. This will prove useful in our upcoming endeavor. We also adopt the notations QW(P) for the marginal distribution of W,QV(P) for the conditional distribution P(V∣W,A0), and Q(P)≡QW(P),QV(P),Qta(P):t=0,…K. We write g for the treatment allocation probabilities P(At∣At−1,Lt−1,V,W). When referring to a generic P∈M, we may sometimes write Q and g in place of Q(P) and g(P), similarly for their respective components. When referring to the functionals at the data-generating distribution P0, we may sometimes write Q0 and g0, in place of Q(P0) and g(P0). When there is no confusion, we will simply write Q and g under the implied distribution.

3.1 Inverse probability weighted estimator ψIPW

Define

for t=0,…,K. We note that these are functionals of g(P), but we have suppressed the notation here.

From Remark 1, a valid estimating function for Ψ is given by

Up to a normalizing matrix M(ψ)−1,DIPW(g,ψ) is a gradient for Ψ under a model Mg where g is known. Note that P0DIPW(g(P0),ψ)=U(ψ0) and the corresponding equation to zero is solved for ψ0=Ψ(P0). Therefore, this is an unbiased estimating function for ψ0 and ψIPW is an unbiased estimator if g=g(P0).

To implement the IPW estimator, we first obtain estimators of the probability of treatment propensity gˆ of g. For an observation in the dataset with A having a value in A, the variable CKA(gˆ) will have numerator equals 1 and probabilities in the denominator given by gˆAj∣Aj−1,Lj−1,V,W; otherwise CKA(gˆ)=0. The estimator ψIPW can be obtained by fitting a weighted logistic regression of YK on ϕ(A,V), with weights CKA(gˆ)h˜(A,V). This ψIPW satisfies PnDIPW(gˆ,ψIPW)=0, and it is consistent if gˆ consistently estimates g. Under standard regularity and empirical process conditions, ψIPW is asymptotically linear with influence function M(ψ0)−1DIPW(g,ψ0). The asymptotic covariance of n(ψIPW−ψ0) can be estimated by the sample covariance matrix ΣIPW of M(ψIPW)−1DIPWgˆ,ψIPW.

Due to its inverse weighting by products of treatment and censoring probabilities, this estimator may be sensitive to near positivity violations, wherein observations with small observed treatment propensity would be given excessive weight. This can be particularly pronounced as the number of time points increase. Substitution estimators like the ones we propose in sections 3.2 and 3.3 can slightly mitigate this problem by incorporating global information in the parameter map. However, the effect of near positivity violations can still take form of poor estimation in the nuisance parameters.

3.2 Non-targeted substitution estimators ψVsubs and ψA0subs

From eq. (8) and Remark 1, we can formulate Ψ as ΨQ0a,QV or ΨQ0a,g. These two formulations lead to two non-targeted substitution estimators, ψVsubs and ψA0subs. Note the latter option opens the door for G-computation estimators in applications with high-dimensional V.

We start by obtaining an estimator of Q0a(V,W) for a given a∈A. We will exploit the recursive property of Qta noted in eq. (9), so that we can avoid estimation of densities, and will instead employ regression-based procedures to estimate the conditional expectations.

1. Initiate at t=K: We first regress YK on LK−1,A,V,W among observations uncensored by time K, to obtain an estimator of EYK∣LK−1,A,V,W. This can be implemented using data-adaptive (e.g. SuperLearner by [15]) or parametric procedures. We then evaluate this fitted estimator at the observed (LK−1,V,W) with setting A=a, to obtain an estimator QˆKa(LK−1,V,W) for each observation uncensored at time K−1.

2. At each t=K−1,…,0, in decreasing order, we would have obtained an estimator Qˆt+1a(Lt,V,W) for observations uncensored by time t+1. We now regress this variable on Lt−1,At,V,W among observations uncensored by time t, and then evaluate the resulting fit at the observed (Lt−1,V,W) with setting At=at among observations uncensored at time t−1.

3. After running step 2 sequentially from t=K−1 down to t=0, we would have an estimator Qˆ0a(V,W) for each of the n observations. We can now use Qˆ0a(V,W) to estimate Ψ through the characterizing equations in Remark 1.

Consider the first representation ΨQ0a,QV. Let QˆV be an estimators of QV(v∣A0=a0,W). For each a∈A and each v∈V, we create a row of data consisting of Qˆ0a(v,W),ϕ(a,v),h(a,v),QˆV for each observation in the dataset. We then pool together these datasets characterized by (a,v). The estimator ψVsubs can be obtained by a weighted logistic regression of Qˆ0a(v,W) onto ϕ(a,v) on the pooled dataset, with weights h(a,v)QˆV. This ψVsubs is consistent if both Qˆ0a and QˆV are consistent. Here we use the empirical distribution to estimate the marginal distribution of W.

Consider now the alternative representation ΨQ0a,g0, from the last equality in Remark 1. Let gˆ be an estimator of g. For each a∈A, we create a row of data consisting of Qˆ0a(V,W),ϕ(a,V),h(a,V) and I(A0=a0)/gˆ(a0∣W) for each observation in the dataset. We then pool together these datasets characterized by a. The estimator ψA0subs can be obtained by a weighted logistic regression of Qˆ0a(V,W) onto ϕ(a,V) on the pooled dataset, with weights h(a,V)I(A0=a0)gˆ(a0∣W). This ψA0subs is consistent if both Qˆ0a and gˆ are consistent. Here we use the empirical distribution to estimate the marginal distribution of (W,A0,V).

While these substitution estimators use weighted regressions to obtain estimators for ψ0, they are less susceptible to extremely large weights than the IPW estimator described above. Indeed, QˆV is bounded within (0,1) while the inverse 1/gˆ utilizes product of only one conditional probability. If the Q components are correctly specified parametric models, then the substitution estimators are also asymptotically linear and converge to a normal random variable at n1/2. However, the parametric modeling assumptions are difficult to justify in scientific investigations. Data-adaptive alternatives may yield estimators which are not n1/2 consistent, which prohibits statistical inference.

3.3 Targeted maximum likelihood estimators

In a non-targeted substitution estimator, we can estimate the nuisance parameters of the data-generating distribution by stratification (nonparametric MLE), by fitting a parametric statistical model, or by using a machine learning algorithm. These estimates are then used to evaluate the parameter of interest. As the number of potential confounders increases, these methods may break down due to the curse of dimensionality, or yield a bias–variance trade off that is not optimal for the parameter of interest. The targeted learning approach aims at developing optimal (n1/2 consistent, asymptotically normal, and efficient) estimators of smooth low-dimensional parameters through the use of data-adaptive machine learning [16].

A targeted MLE adds a loss-based targeting step that aims to modify the initial estimators towards the parameter of interest, and provides potential robustness and semiparametric efficiency gains. The targeting is guided by an influence function for the target parameter eq. (8). As a result of this targeting step, the final estimates of the nuisance parameters and of the target parameter satisfy a user-chosen score equation. Under certain regularity conditions, this yields inference based on the Central Limit Theorem. For example, it has been shown that both TMLE and the estimating equation based estimators (using the EIF as estimating equation) are double robust asymptotically efficient estimators under a Donsker class condition and a second order reminder condition [16]. The Donsker class condition roughly states that the EIF at Qn and gn falls in a Donsker class such as multivariate real valued cadlag functions with uniformly bounded sectional variation norm. The second order condition states that the product, of (a) the L2 norm of the estimated outcome regressions minus their true counterpart with (b) the L2 norm of the estimated treatment/censoring mechanism minus their true counterpart, converges to zero at a a rate faster than n1/2. The advantage of the TMLE is that it does not rely on the existence of a unique solution to the estimating equation. The TMLE estimators in this section are asymptotically linear and efficient under these conditions. They also allow different combination of rates for Qn and gn. For example, if g is known, then only consistency is needed for asymptotic efficiency.

Another advantage of the targeted learning approach is that it does not rely on the existence of a unique solution of the estimating equation. Thus for marginal structural working models that are non-linear in parameters, the non-targeted estimating equations estimators may have multiple solutions or no solution. We refer to [2, 3] for the general methodology.