Doubly Robust and Efficient Estimation of Marginal Structural Models for the Hazard Function

1.1 Organization of paper

This paper is organized as follows. Section 2, describes the data structure and defines the parameter using a nonparametric causal framework. Section 3 first reviews the non-targeted substitution estimator (a.k.a. G-computation estimator) and the IPTW estimator for the parameters of interest, and then presents the efficient influence curve, and describes the proposed TMLE estimator. Section 4 evaluates the performance of these three estimators in a simulation study mimicking an observational cohort study. This paper concludes with a summary.

2.1 Causal model and causal parameters

The causal assumptions implied in the notation of O can be made explicit by a nonparametric structural equations model (NPSEM, Pearl [4]):

This causal model assumes that each variable X in the observed data structure is an unknown deterministic function fX of certain observed variables, which we refer to as the parents of X and denote by Pa(X), and some unmeasured exogenous random factors UX. This causal model defines a random variable with distribution PO,U on a unit.

In practice, additional exclusion restriction assumptions informed by subject matter knowledge can be represented by limiting the variables in Pa(X). In the HIV example, at each time t, the investigators may specify that CD4 counts, death, and regimen modification all depend on the patient’s entire observed past and unmeasured factors. But if they know that the decision to switch regimen is based only on the most recent CD4 measurement, then Pa(At1) may be restricted to exclude all earlier CD4 measurements.

As we alluded to in Section 1, research questions can be formalized in terms of intervention rules on the treatment variables. An intervention rule d is a function that deterministically assigns treatment at time t based on covariate history according to At=dLt. This rule may be static – assigning the same treatment option at regardless of covariate history, i. e. At=dLt≡at; or it may be dynamic – assigning different treatment options to different covariate histories. We will use the boldface notation At=dLt to denote the vector Ak=dLkk=0,…t; this means that each of the variables Ak, from k=0 to t, was assigned value according to rule d and its covariate history. In the HIV example, a static rule would be to always switch regimen at m months after virologic failure, and a dynamic rule would be to switch ART at the first time t when the patient’s CD4 counts drop below a pre-specified threshold θ.

Given a set d of intervention rules of interest, investigators are often concerned about their comparative causal effect on the outcome process Yt. To be more precise, consider an ideal experiment where all subjects are assigned treatment under an intervention rule d and right censoring is also prevented; the covariates, on the other hand, take the value that they may in response to rule d. We call the variables At the intervention variables, as they are the ones that are subject to manipulation. This ideal experiment can be formalized in the NPSEM by setting the equations for At to At=dLt, and replacing the input At−1 in fLt with dLt−1, which are treatment assignments from time 0 to t−1 under rule d and history Lt−1. As a result, the only random endogenous variables of the system are the covariates; we use Lt(d) to denote the time varying covariates that result under the intervention rule d, in particular Yt(d) is the indicator of death under such a regimen. The comparative causal effect of rule d1 vs rule d2 can be assessed by comparing the distributions of the outcome processes Yt(d1) vs Yt(d2).

Suppose we wish to understand how the outcome process Yt(d) changes as a function of d, t and some baseline covariate V⊂L0. To this end, we study the intervention-specific hazard function

on the space d×τ×v, where d is the set of rules we wish to compare, τ={1,…K+1} contains all the follow-up times of interest, and v is the outcome space of V.

Studying the entire function (d,t,V)↦λ(d,t,V) is often difficult due to feasibility of computation, interpretation of results, and other challenges associated with high-dimensionality. Instead, we can study a more tractable, simplified model/summary of this function which captures how λ changes as a function of (d,t,V) in only a few summarizing parameters. More specifically, following the marginal structural working model framework in Neugebauer and van der Laan [25], let mψ(d,t,V) be a working MSM for λ(d,t,V), parameterized by ψ∈RJ. To be concrete, from here on, we will use a generalized linear hazard MSM with logit link,

where ϕ(d,t,V) is the vector of linear predictors. However, it is important to note that the methods presented here can be generalized to other working MSM.

In addition to specifying the working MSM, we also specify a user-supplied kernel weight function h(d,t,V), and the standard log-likelihood loss. Formally, the causal parameter of interest is defined as

In words, our causal parameter of interest is the vector of coefficients of the MSM weighted by h. Under the working MSM framework, this parameter can be understood as the best (assessed by the loss function) weighted (by h) approximation (given by the working model mψ) of the hazard function λ. Embedded in the definition of the parameter are the stabilizing weights h(d,t,V). To mitigate positivity violations, h(d,t,V) can be chosen to down weight a rule d which has little data support at time t within strata of marginal covariates V. But the effectiveness of this stabilization will be limited when d uses time-varying tailoring variables.

Continuing our HIV example, suppose we wish to assess how delay in switching regimen affect mortality. Let dm denote the rule that dictates switching regimen at m months after virologic failure. That is, At=dm(Lt)=0 for t<m and At=dm(Lt)=1 for t≥m, given Yt=0. Let d be the set of all possible switching times we wish to compare, i. e. d={dm:m=0,…,K,∞}. If we are only interested in the marginal hazard, then we set V=∅. Or, if we wish to assess the hazard function stratified by CD4 at the time of virologic failure, then we can set V=1CD40<500cells/μℓ. For simplicity, let us continue this example with the former option, and let the weights be h(d,t)=1. For this example, we choose the logistic working model to be mψ(dm,t,V)=expitψ0+ψ1t+ψ2max(t−m,0), where max(t−m,0)=0 if the patient has not switched by time t−1, and max(t−m,0)=t−m encodes how long the patient has been on his second line regimen if he has already switched.

2.2 Statistical parameter

Thus far, we have used the NPSEM to formulate the parameter of interest Ψ(PO,U) in eq. (2). This parameter is a function of the distributions of Lt(d) and Yt(d), which are generated within an ideal experiment. Unfortunately, such ideal experiments are not always possible in real life. Then, what are the sufficient assumptions on the data-generating process under which the parameter Ψ(PO,U) can be estimated using the observed data?

To answer this question, we review the expression in eq. (2) and note that P(Yt−1(d)=0|L0)=1−E(Yt−1(d)|L0) and

Consequently, to identify Ψ(PO,U) as a function of the data generating distribution P0, it suffices to establish the identification of the time-dependent causal dose-response curve {E(Yt(d)|L0):d∈d,t∈τ}.

To this end, we make the sequential randomization assumption (Robins [3]):

and the positivity assumption

for every d∈d, k∈τ, and V∈v with h(d,t,V)>0. Assumption (4) specifies that at each k, the variable Ak is randomized conditional on its parent variables. In our HIV example, this can be satisfied if the parent variables of Ak contain all common determinants of mortality and the decision to switch regimen. Assumption (5) requires that under P0, for each rule d and its compatible covariate history for which the weight h(d,k,V) is non-zero, there is non-zero probability that the subject’s treatment will continue to follow this rule. In our HIV example with h(d,t,V)=1, given one has not yet switched regimen, there should be non-zero probability that the patient will switch at any given time m in d. Assumption (5) would be violated if, say, d contains the rule to switch at 6 months but no patient in the study population actually switches at 6 months post-failure. This assumption would also be violated if, say, all patients with switch regimen at month 1.

Under assumptions (4) and (5), we obtain identification of the intervention-specific mean:

We are reminded that ℓ1,j−1≡ℓ1,...,ℓj−1. This is generally known as the G-computation formula (Robins [3]) for the intervention-specific mean. This in turn identifies the intervention-specific hazard function:

Consequently, the causal parameter of interest Ψ(PO,U) in eq. (2) is identified as a function of the observed data distribution given by

Note that the sequential randomization assumption (4) is an untestable assumption on the data-generating process under which we can claim the causal parameter (2) equals the statistical parameter (8). It provides a guiding principle to judge the causal interpretation of the statistical parameter. By contrast, the positivity assumption (5) is a testable assumption on the data-generating distribution P0 under which the statistical parameter (8) is well-defined.

The remainder of this paper focuses on statistical inference for ψ0.

3.1 Notations

We use Pn to denote the empirical distribution of n i.i.d. copies of O∼P0. Given a function O↦f(O), Pnf denotes the empirical mean Pnf≡1n∑i=1nf(Oi). More generally, for any P∈m, Pf≡EPf(O).

For a generic P∈m. We use QL0 to denote the marginal distribution of L0. Generalizing the functional in eq. (6), for a given t∈τ and d∈d, denote at P∈m:

and Qt+1d,t(⋅)≡Yt. In the above notation, the superscript t signals the expectant outcome variable Yt and the subscript k signals the length of the conditioning covariate history. Under assumptions (4) and (5), Qkd,t(P0)(Lk−1)=E(Yt(d)|Lk−1,Ak−1=d(Lk−1)), which is the conditional mean of Yt(d) given an observed past that has followed the rule d up to time k−1. For simplicity, we may sometimes write Qkd,t instead of Qkd,t(P) when referring to the functional evaluated at a generic P∈m, and Qk,0d,t(Lk−1) instead of Qkd,t(P0)(Lk−1) when the functional is evaluated at P0∈m. Bang and Robins [18] made a key observation that these functionals satisfy the relation

For a counting outcome process, these functionals are also monotonic in t, i. e. for a given k, Qkd,t≤Qkd,t+1 for all t≥k. We denote the components corresponding to the non-intervention variables as Q≡QL0,Qkd,t:,t∈τ,k≤t,d∈d. For the intervention variables, we denote the treatment allocation probabilities P(Ak|Lk,Ak−1) as g(Ak|Lk,Ak−1), and the product ∏k=1tg(Ak|Lk,Ak−1) as g(At|Lt). For our purposes, the couple (Q,g) readily specifies a distribution P∈m, so sometimes we may abuse notation and write P=(Q,g). At the data generating distribution P0, we adopt the subscripts Q0 and g0.

3.2 G-computation Estimator

The G-computation formula in eq. (6) readily delivers a non-targeted substitution estimator (as opposed to the targeted substitution estimator that is TMLE), which is generally known as the parametric G-computation (Gcomp) estimator. More precisely, using the notations in Section (3.1), the statistical parameter Ψ(P0) in eq. (8) can be expressed as Ψ(Q0). Therefore, fitting the hazard MSM using a non-targeted estimator Qn of Q0 yields a non-targeted substitution estimator Ψ(Qn) of Ψ(Q0).

Recall that Q0 consists of the marginal distribution of L0, QL0(P0), and the conditional means of Yt(d), Qkd,t(P0) defined in eq. (9). To estimate the marginal distribution QL0(P0), we can use the empirical distribution of L0, denoted QnL0. To estimate the conditional means Qkd,t(P0), one approach is to estimate the conditional densities for each Lt given its parents and use the definition in eq. (9). While this density-based approach ensures that the monotonicity in t of Qkd,t is preserved, the dimension of the nuisance parameter and the computational cost grows with the dimension of Lt and the number of time points. One way to implement this approach is to use simplifying parametric modeling assumptions on these nuisance parameters. We refer to Taubman et al. [26] and Young et al. [27] for expositions of this technique.

Another approach to estimate the conditional means Qkd,t(P0) is to minimize estimation of nuisance parameters by exploiting the recursive relation eq. (10) noted by Bang and Robins [18], which allows one to use the many available regression techniques (parametric or data-adaptive) in the literature. We can implement this regression-based approach by running the following two-level algorithm. Recall that in our notation Qkd,t(P0), d is the intervention of interest, t is the time of the expectant outcome, and k is the length of the conditioning history.

1. Starting at t=K+1, estimate the conditional means Qkd,t(P0): k≤t,d∈d as follows:

   1. Initiate at k=t: Recall that Qtd,t(P0)≡EP0Yt|Lt−1,At−1=d(Lt−1) and Qt+1d,t≡Yt. We obtain an estimator Qt,nd,t of Qtd,t(P0) by regressing Yt on the observed values Lt−1 and At−11 among observations that remain uncensored by time t−1, and evaluate the fitted function at the observed Lt−1 and the intervened values At−11=d(Lt−1), for these uncensored observations. For a subject that has died before time t, the value is deterministically set to 1. We can organize the data by having one row for each patient i and each regimen d.

   2. At each subsequent k=t−1,…,1: At the previous step, we have thus far obtained an estimator Qk+1,nd,t of Qk+1d,t(P0). To obtain an estimator Qk,nd,t of Qkd,t(P0), we regress Qk+1,nd,t(Lk) on the observed values Lk−1 and Ak−11 among those observations that remained uncensored by time k−1, and evaluate the fitted function at the observed Lk−1 and the intervened values Ak−11=d(Lk−1) on these uncensored observations. For observations that had died before time k, the value is deterministically set to 1.

   3. After iterating step (b) in order of decreasing k, we have obtained a sequence of means Qk,nd,t:k≤t,d∈d.

2. Repeat step 1 in order of decreasing t, from t=K to t=1. At the end, we will have obtained estimators Qk,nd,t:t∈τ,k≤t,d∈d

At the end of this algorithm, we have the sequence of conditional means. {Q1,nd,t(L0): t∈τ,d∈d} for each of the n observations. Pooling together these estimates Q1,nd,t(L0) over all d and t, we have one row per patient i, rule d∈d and t∈τ. The G-computation estimator ψnGcomp≡Ψ(Qn) is obtained by fitting a weighted logistic regression of Q1,nd,t(L0)−Q1,nd,t−1(L0)1−Q1,nd,t−1(L0) according to the MSM, with weights h(d,t,V)(1−Q1,nd,t−1(L0)).

Consistency of ψnGcomp relies on consistency of Qn. In the density-based approach, this means consistent estimation of all the conditional densities of Lt given its parents; in the regression approach, this means consistent estimation of the conditional means Qkd,t(P0). In either case, correct specification of Q0 under a finite dimensional parametric model is rarely possible in practice. Alternatively, we may use machine learning algorithms, such as Super Learner. This option is more enticing, especially when used with the regression-based approach, since there are more data-adaptive techniques available to estimate a conditional mean via regression. However, unlike the Inverse Probability of Treatment Weighted Estimator (Section 3.3) and the Targeted Maximum Likelihood Estimator (Section 3.4) which satisfy a central limit theorem under appropriate empirical process conditions, analogous results for ψnGcomp are not available, regardless of model size. In addition, a non-targeted estimator Qn of Q0 is obtained by minimizing a global loss function for Q0, not for Ψ(Q0). This means, in particular, that the bias-variance tradeoff in Qn is optimized for the high-dimensional nuisance parameter Q0, instead of a much lower-dimensional parameter of interest Ψ(Q0). The proposed targeted estimator in Section 3.4 aims to address these two issues by providing a substitution estimator that is asymptotically linear (under appropriate regularity conditions), and optimizes the bias-variance tradeoff of Qn towards Ψ(Q0) via an updating step.

3.3 Inverse probability of treatment weighted estimator

Inverse probability of treatment weighted estimation is the standard methodology for estimating the parameters of a marginal structural model in the presence of time-varying confounding [1]), due to its ease of implementation and its asymptotic linearity, which allows for construction of Wald-type confidence intervals.

To begin, we first note that the parameter in eq. (8) can be rewritten in an IPTW form:

The IPTW estimator ψnIPTW is obtained by fitting a weighted logistic regression of Yt according to the MSM, with weights h(d,t,V)I(At−1=d(Lt−1))gnAt−1=d(Lt−1)|Lt−11−Yt−1, where gn is an estimator for g0.

The asymptotic theory of the IPTW estimator is well understood in the literature. We refer the reader to Robins [12], van der Laan and Robins [28] and van der Laan and Petersen [7], where the last reference specifically addresses dynamic intervention rules. In summary, ψnIPTW described above satisfies the estimating equation PnDIPTW(ψ,gn)=0, where

with h˜(d,t,V)≡h(d,t,V)ϕ(d,t,V). The estimator ψnIPTW is well-defined if it is a unique minimizer, and it is a consistent estimator of ψ0 provided gn is a consistent estimator of g0. Moreover, ψnIPTW is asymptotically linear with influence curve

where

The sample covariance matrix of the estimated influence curve is given by ΣnIPTW≡MIPTW(ψnIPTW,gn,Pn)−1DIPTW(ψnIPTW,gn). The variance of n(ψnIPTW−ψ0) can be estimated using ΣnIPTW. Consequently, we can construct Wald confidence intervals of level (1−α) as ψnIPTW±ξ1−α/2(ΣnIPTW/n)1/2, where ξ1−α/2 is the (1−α/2)-quantile of the standard normal distribution. These confidence intervals ought to be conservative when g0 is estimated using a correctly specified parametric model. But no such theoretical guarantee for data-adaptive estimation of g. Note that this influence curve based variance estimate assumes that the weights h(d,t,V) are known functions; when these weights are estimated, this variance estimate should be interpreted as estimating the variance of an estimator of MSM parameters defined by the estimated weights.

Though both G-computation estimator and IPTW estimator properly account for time-varying confounding, the popularity of IPTW over G-computation estimator is apparent from its straightforward implementation and its theoretical validity for confidence intervals. Moreover, the treatment probabilities g0 may arguably be easier to specify correctly than the sequential conditional means Qkd,t(P0). However, the IPTW estimator may be generally more susceptible to near positivity violations due to the inverse probability weighting in eq. (11). Often, the kernel functions h(d,t,V) are chosen to stabilize these inverse probability weights. This remedy, however, is less effective when the rules d are dynamic and use time varying covariates as tailoring variables, since the marginal function h(d,t,V) cannot depend on the time varying covariates Lt.

3.4 Targeted maximum likelihood estimator

As discussed earlier, consistency of the IPTW estimator relies on consistency of gn, while consistency of the G-computation estimator relies on consistency of Qn. In this section, we propose a semiparametric efficient estimator that is robust against misspecification of either Q0 or g0. These theoretical promises hinge on the use of one important ingredient – the efficient influence curve for ψ0.