Structural Nested Mean Models to Estimate the Effects of Time-Varying Treatments on Clustered Outcomes

3.1 Structural nested mean models

For binary treatments, we may consider the joint treatment {A1,A2} across units within a cluster as one combined treatment at the cluster level with 4 levels, namely {0,0}, {1,0}, {0,1}, {1,1}. Then a SNMM with joint treatment can be formulated in the same way as one with a single treatment as a contrast between means of potential outcomes under joint treatment versus means of potential outcomes under no treatment as seen in eq. (1), where g(.) is a monotone link function. If there is an additional binary systemic treatment A†, the combined treatment {A†,A1,A2} will contain eight levels. In the rest of the paper, we focus on the SNMM with identity link such that g(x)=x. The method introduced, however, can be applied to other link functions, which will be discussed in the discussion section:

The vector-valued blip function γ∗(.) models the effect of removing joint treatment a from both unit-level outcomes within the cluster. It can be parameterized with the finite dimensional parameter ψ to form

As an example, let

In eq. (3), ψ∈R6; ψ1 and ψ1′ represent the effect of unit-specific treatment on the outcome in the same unit. ψ2 and ψ2′ represent the effect of treatment on the outcome in the opposite unit, ψ3 and ψ3′ represent the interaction between treatments in different units. We can also interpret ψ1 and ψ1′ as the direct effect of treatment on the same unit. ψ2 represents the spillover effect of unit 2’s treatment on unit 1 when a1=0, and ψ2+ψ3 represents the spillover effect of unit 2’s treatment on unit 1 when a1=1. Similarly, ψ2′ represents the spillover effect of unit 1’s treatment on unit 2 when a2=0, and ψ2′+ψ3′ represents the spillover effect of unit 1’s treatment on unit 2 when a2=1. This example imposes the assumption that effect modification by covariates is absent.

The dimension of causal parameters will increase sharply with increase in the number of units per cluster. Some additional model assumptions are considered to reduce the dimension of causal parameters. For the ophthalmology example, it is reasonable to assume symmetry between two eyes, namely the effect of the joint treatment is the same for each unit. In eq. (3), the symmetry assumption implies ψ1=ψ1′,ψ2=ψ2′,ψ3=ψ3′. In some cases, we may want to assume no interference, namely unit-specific treatment on one unit has no effect on outcomes in the other units of the same cluster. With two units per cluster, it refers to

In eq. (3), the no interference assumption implies ψ2=ψ3=ψ2′=ψ3′=0.

3.2 Identifiability

A1–A4 below are sufficient assumptions to identify the causal effect of joint treatment a relative to no treatment from observed data and therefore the causal parameter ψ.

1. Consistency. If A=a for a given cluster, then Ya=Y for that cluster.

2. Conditional exchangeability. Y0,0⊥ A | L=l.

3. Positivity. If fL l > 0, then fA|L(a|l)>0 for all treatments a.

4. Cluster-level SUTVA. The observation on one cluster is unaffected by the particular treatment assignment to the other clusters.

In the presence of interference between units, potential outcomes like Yjaj are not well-defined. However, under no interference, it is reasonable to relax A1 and define consistency under the standard assumption for non-clustered data: if Aj=aj for a given unit j, then Yjaj=Yj for that unit. A2 can be viewed as no unmeasured confounding. It states the independence between treatment-free potential outcomes and the joint treatment to a cluster conditional on the covariates within the cluster. A3 can be relaxed in SNMMs such that the assumption does not need to hold for all L with fL l>0. A4 is similar to the partial interference assumption described in Sobel [10]. Note that our method allows for interference between units within the same cluster, but the SUTVA assumption needs to hold for clusters.

3.3 Alternative modeling

An alternative way of modeling the joint effect of {a1,a2} is to use the concepts in SNMMs with time-varying treatments and sequentially remove the effect of each unit-specific treatment. In the following example, the effect of treatment on the opposite unit is removed first, followed by the effect of treatment on the same unit.

1. Remove the effect of treatment on the opposite unit aj′:

   (4a)EYjaj,aj′|L=l,A=a−EYjaj,0|L=l,A=a=γ1,j∗(aj′,l,aj)

2. Remove the effect of treatment on the same unit aj:

   (4b)EYjaj,0|L=l,Aj=aj−EYj0,0|L=l,Aj=aj=γ2,j∗(aj,l)

For example, the blip functions in eqs (4a) and (4b) can be parameterized as γ1,j∗(aj′,l,aj;ψ)=ψ2aj′+ψ3aj′aj and γ2,j∗(aj,l;ψ)=ψ1aj, respectively.

Assumptions A1 and A3 plus a sequential ignorability assumption are sufficient to identify the causal effect in eqs (4a) and (4b). Specifically, the sequential ignorability assumption states that Yjaj,0⊥Aj′|L,Aj and Yj0,0⊥Aj|L. Testing for no interference under this model is equivalent to testing whether the blip function in eq. (4a) equals 0. In this case, since aj and aj′ are simultaneous treatments, the order of blipping down is arbitrary. Parameters from different orders of blipping down have different interpretations. For this reason, the approach in eq. (1) is used for the rest of the paper. However, the alternative approach may be useful in other settings where the order of blipping down is meaningful. For instance, if clusters are defined by families, it may be reasonable to blip down the treatment on parents first and then that on children or vice versa.

3.4 Estimation

Robins [1] has derived a class of estimating equations for parameters of SNMMs based on semiparametric theory [1]. He defined Hk(ψ)=Yk−∑l=0k−1γl,k∗Lˉl,Aˉl;ψ as the fully blipped down outcome at time k, and H(m)=[Hm+1(ψ),Hm+2(ψ),....,HK(ψ)]T as the vector of blipped down outcomes after time m, where m=0,...,K−1. Under the sequential ignorability assumption, EH(m)|Lˉm,Aˉm=EH(m)|Lˉm,Aˉm−1. H˙(m) is defined as H˙(m)=H(m)−EH(m)|Lˉm,Aˉm−1. The efficient estimating function, known as the efficient score [11], under general settings involves complex recursive expressions, but under the assumption that

for k>m,j>m,m=0,...,K−1, the efficient score can be simplified as:

By setting eq. (6) to 0 and solving the estimating equation, one can obtain an efficient estimator for the causal parameter ψ. If the assumption in eq. (5) does not hold, the estimator from eq. (6) will not achieve semiparametric efficiency but will nonetheless be consistent. Equivalently, we can express eq. (6) in terms of the partially blipped down outcome Um with removal of the effect of treatments from time m onwards [4]. Specifically, Um is a vector with elements Um,k=Yk−∑l=mk−1γl,k∗Lˉl,Aˉl, for k=m+1,...,K.Um mimics the potential outcomes in the sense that E(Um|L¯m,A¯m)=E(Y_m+1a¯m−1,0|L¯m,A¯m). The alternative expression is

The theory for SNMMs can be directly applied to situations with clustered observations. We first construct blipped down outcomes U. For a single time point example with bivariate outcomes, they can be written as

The following is a parameterized example for blip function under the symmetry assumption:

Based on eq. (7), the efficient score for a single time point example can be written as

For the example in eq. (8), ∂U∂ψT has the following expression:

E∂U∂ψT|L is based on a model for treatment. If treatment probabilities are unknown, a multinomial logistic regression model can be used to estimate the treatment probability in each category of the joint treatment. The estimating equation in eq. (9) is doubly robust, because the estimator for ψ will be consistent if either the model for E∂U∂ψT|L or the model for EU|L is correctly specified. When both are correctly specified and the conditional variance VarU|L is also correctly specified, the estimator achieves the semiparametric efficiency bound for model (2) if assumption (5) is true. However, correct specification of VarU|L is not required for the consistency of the estimator.

4.1 Structural nested mean models

In the SNMM below, the blip function vector γm∗(.) models the effect of removing joint treatment am at time m from both outcomes, after having removed the effect of all subsequent treatments.

For repeatedly measured outcomes, eq. (10) can be further generalized to

where γ_m∗ is a vector with components γm,k∗(aˉm,lˉm;ψ) that models the effect of removing joint treatment at time m from both outcomes at time k, for k=m+1,...,K. A parameterized example for eq. (10) is shown below:

The ψ parameters in eq. (12) have similar interpretation as those in example (3), except that they contain subscript m to indicate the effect of joint treatment at time m is time-specific. A more parsimonious example is shown below:

Model (13) makes the symmetry assumption and also assumes the effect of treatments is the same for each time point. Notice that subscript m is omitted for ψ. The ψ parameters can be interpreted as effect of treatment per unit time. For example, ψ1 represents the effect of unit-specific treatment on the outcome in the same unit per unit time.

The sequential ignorability assumption is required to identify the contrast in eqs (10) and (11). It is a generalization of assumption A2 to a series of time-varying treatments as shown below:

It states that there is no unmeasured confounding between the partially blipped down potential outcomes and the joint treatment to a cluster at time m, conditional on the covariate and treatment history within the cluster.

4.2 Estimation

Similar to the single time point case, we first construct a series of partially blipped down outcomes Um for each time m:

The partially blipped down outcomes mimic the potential outcomes Yaˉm−1,0 in the sense that EUm|Lˉm,Aˉm=EYaˉm−1,0|Lˉm,Aˉm. For repeatedly measured outcomes, Um is a vector with components

for k=m+1,...,K.

For time-varying treatments, we suggest using the following score:

The efficient score in eq. (6) is usually too complicated to be used for situations with time-varying treatments. The estimator from the score in eq. (14) is not optimal but is reasonably efficient and doubly robust. It is straightforward to show that the estimator is consistent, if for each time m, either the model for treatment, E∂γm∗Lˉm,Aˉm;ψ∂ψT|Lˉm,Aˉm−1, or the model for outcome, EUm|Lˉm,Aˉm−1, is correctly specified. A simpler non-doubly robust estimating equation is

It does not involve modeling EUm|Lˉm,Aˉm−1 and assumes it to be 0. Therefore, the estimator from eq. (15) is consistent if and only if the model for treatment is correct.

For the parameterized example in eq. (13), we have the following expression for ∂γm∗Lˉm,Aˉm;ψ∂ψT:

Similarly, E∂γm∗Lˉm,Aˉm;ψ∂ψT|Lˉm,Aˉm−1 is based on a model for treatment. If treatment probabilities are unknown, a multinomial logistic regression model can be used to estimate treatment probability in each category of the combined treatment {Am,1,Am,2} received at time m. We can either fit a separate model for each time m, or fit a single model to pooled data from all person times. For the pooled treatment model, we may allow intercepts to be time-varying with the use of a spline.

4.3 Closed form estimator

If VarUm|Lˉm,Aˉm−1 is given, and a linear working model Dmδ=dmLˉm,Aˉm−1δ is assumed for EUm|Lˉm,Aˉm−1, where δ is the q−dimensional nuisance parameter in the outcome model and dm is a J×q matrix-valued function of Lˉm,Aˉm−1, and the blip function is linear in ψ with γm∗Lˉm,Aˉm=Rmψ=rmLˉm,Aˉm, where ψ is p−dimensional and rm is a J×p matrix-valued function of Lˉm,Aˉm, the estimator for ψ will have the following closed form (based on Robins and Hernán [12]):

where Cm=∑l=mK−1Rl−ERl|Lˉm,Aˉm−1.

A simple working model for VarUm|Lˉm,Aˉm−1 is σ2IJ×J for some constant variance σ2 at all times, where IJ×J is an identity matrix of dimension J. However, the independence structure and constant variance over time are often far from truth. A simulation study described in the next section will examine the impact of such misspecification on the efficiency of ψˆ. If one assumes a working model other than independence structure and constant variance, estimation will then require iterative procedures similar to those in generalized estimating equations (GEEs). We can start by assuming working model σ2IJ × J to obtain an initial estimate for ψ, and then iterate between estimating parameters in VarUm|Lˉm,Aˉm−1 and estimating ψ until convergence is achieved. We adapted the method of moments approach by Liang and Zeger [13] in estimating variance parameters [13]. An asymptotic variance estimator for ψˆ may be derived using the sandwich method, with the exact form found in Appendix A. Like GEE, VarUm|Lˉm,Aˉm−1 does not need to be correctly specified for the estimator for ψ to be consistent.

5.1 Algorithm

The data generating process is based on the likelihood for SNMMs described in Robins et al. [9]. Details about the likelihood can be found in Appendix B. A simulation study is conducted based on the setting where treatments and covariates are measured repeatedly at times m=0,...,K−1, and outcomes are measured at the last time point K only. The number of time points is K=5, and the cluster size is J=2. The total number of clusters is varied as N=100 or 1,000. The algorithm for generating data for one cluster is shown below. It includes two examples with different covariate effects in the conditional mean of outcomes. Variables with subscript m=−1 assume value 0.

For m=0,...,K−1, we generated covariates Lm from a multivariate normal distribution:

We generated joint treatment Am from a multinomial logistic regression model with probabilities pm,00, pm,10, pm,01, pm,11:

where

Last, observed outcomes YK were generated. We first generated the conditional mean of YK under the sequential ignorability assumption: