Marginal Structural Models and Causal Inference in Epidemiology

In observational studies with exposures or treatments that vary over time, standard approaches for adjustment of confounding are biased when there exist time-dependent confounders that are also affected by previous treatment. This paper introduces marginal structural models, a new class of causal models that allow for improved adjustment of confounding in those situations. The parameters of a marginal structural model can be consistently estimated using a new class of estimators, the inverse-probability-of-treatment weighted estimators.

Marginal structural models (MSMs) are a new class of causal models for the estimation, from observational data, of the causal effect of a time-dependent exposure in the presence of time-dependent covariates that may be simultaneously confounders and intermediate variables. ^1-3 The parameters of a MSM can be consistently estimated using a new class of estimators: the inverse-probability-of-treatment weighted (IPTW) estimators. MSMs are an alternative to structural nested models (SNMs), the parameters of which are estimated through the method of g-estimation. ^4-6

The usual approach to the estimation of the effect of a time-varying exposure or treatment has been to model the probability of disease as a function of past exposure and past confounder history, using analytic methods such as stratified analysis and its parametric analogs (for example, logistic or proportional hazards regression). We will show in sections 4 and 7.1 that these standard approaches may be biased, whether or not one further adjusts for past confounder history in the analysis, when (1) there exists a time-dependent covariate that is a risk factor for, or predictor of, the event of interest and also predicts subsequent exposure, and (2) past exposure history predicts subsequent level of the covariate. We refer to covariates satisfying condition 1 as time-dependent confounders. Conditions 1 and 2 will be true in many observational studies, particularly those in which there is confounding by indication. For example, in a study of the effect of zidovudine (AZT) treatment on mortality among human immunodeficiency virus (HIV)-infected subjects, the time-dependent covariate CD4 lymphocyte count is both an independent predictor of survival and initiation of therapy with AZT and is itself influenced by prior AZT treatment. In a study of the effect of obesity on mortality, the development of clinical cardiac or respiratory disease is an independent predictor of both mortality and subsequent weight loss and is influenced by prior weight gain. Conditions 1 and 2 will always hold when there are time-dependent covariates that are simultaneously confounders and intermediate variables. A more detailed description of the bias of standard methods, as well as several additional epidemiologic examples of time-dependent confounding, has been presented elsewhere. ^5,7

1. Time-Dependent Confounding

Consider a follow-up study of HIV-infected patients. Let A_k be the dose of the treatment or exposure of interest, say AZT, on the k^th day since start of follow-up. Let Y be a dichotomous outcome of interest (for example, Y = 1 if HIV RNA is not detectable in the blood and is 0 otherwise) measured at end of follow-up on day K + 1. Our goal is to estimate the causal effect of the time-dependent treatment A_k on the outcome Y.

Figure 1 is a causal graph that represents our study with K = 1. A causal graph is a directed acyclic graph in which the vertices (nodes) of the graph represent variables and the directed edges (arrows) represent direct causal effects. ⁸ In Figure 1, L_k represents the value on day k of the vector of all measured risk factors for the outcome, such as age, CD4 lymphocyte count, white blood count (WBC), hematocrit, diagnosis of acquired immunodeficiency syndrome (AIDS), and the presence or absence of various symptoms or opportunistic infections such as oral candidiasis. Similarly, U_k represents the value on day k of all unmeasured causal risk factors for Y. Figure 1, b, differs from Figure 1, a, only in that the arrows from the unmeasured causal risk factors into the treatment variables have been removed. When, as in Figure 1, b, there is no arrow from unmeasured causal risk factors into treatment variables, we say that there are no unmeasured confounders given data on the measured confounders L_k. ^9,10 Figure 1, c, differs from Figure 1, a and b, in that none of the risk factors for Y (measured or unmeasured) has arrows into any treatment variable. Note, however, that earlier treatment A₀ can causally affect later treatment A₁. When, as in Figure 1, c, there is no arrow from any (nontreatment) risk factor into any treatment variable, there is no confounding by either measured or unmeasured factors, in which case we say that treatment is unconfounded. ^9,10

Figure 1

The distinctions drawn above apply equally to more familiar point-treatment studies in which the treatment is not time-dependent. As indicated in Figure 2, a point-treatment study is a special case of the general set-up in which K = 0. Figure 2, a-c, contains the analogs of Figure 1, a-c, for a point-treatment study.

Figure 2

As in any observational study, we cannot determine from the observed data on the L_k, A_k, and Y whether there is confounding by unmeasured risk factors. We can only hope that whatever residual confounding there may be due to the U_k is small. Under the untestable assumption that there is no unmeasured confounding given the L_k, we can, however, empirically test from the data whether treatment is unconfounded. Specifically, a sufficient condition for treatment to be unconfounded is that, at each time k, among subjects with the same past treatment history A₀, …, A_k-1, the treatment A_k is unassociated with the past history of measured covariates L₀, …, L_k. ^9-11 For example, in our point-treatment study, treatment will be unconfounded if A₀ is unassociated with L₀.

2. Counterfactuals in Point-Treatment Studies

We begin with a review of how one would estimate the effect of A₀ on Y in the point-treatment study of Figure 2. Suppose treatment A₀ is dichotomous; suppose further that Figure 2, c, is the true causal graph, that is, that neither measured nor unmeasured covariates confound the relation between treatment and the outcome. Then the crude risk difference, risk ratio, and odds ratio each measure the causal effect of the treatment A₀ on the outcome Y, although on different scales. The crude risk difference is cRD = pr[Y = 1|A₀ = 1] - pr[Y = 1|A₀ = 0], the crude risk ratio is cRR = pr[Y = 1|A₀ = 1]/pr[Y = 1|A₀ = 0], the crude odds ratio is cOR = pr[Y = 1|A₀ = 1]pr[Y = 0|A₀ = 0]/{pr[Y = 1|A₀ = 0]pr[Y = 0|A₀ = 1]}, and, for example, pr[Y = 1|A₀ = 1] is the probability that Y = 1 among treated subjects (A₀ = 1). We assume that the study subjects are a random sample from a large, possibly hypothetical, source population. Probabilities refer to proportions in the source population.

The causal contrasts that correspond to these associational parameters involve counterfactual variables. Specifically, the variable Y_a0=1 denotes a subject's outcome if treated and Y_a0=0 denote a subject's outcome if left untreated. For a given subject, the causal effect of treatment, measured on a difference scale, is Y_a0=1 - Y_a0=0. For no subject are both Y_a0=0 and Y_a0=1 observed. If a subject is treated (A₀ = 1), the subject's observed outcome Y equals Y_a0=1, and Y_a0=0 is unobserved. If A₀ = 0, Y equals Y_a0=0, and Y_a0=1 is unobserved. Let pr(Y_a0=1 = 1) and pr(Y_a0=0 = 1), respectively, be the probability that Y_a0=1 is equal to 1 and Y_a0=0 is equal to 1. Then, if treatment A₀ is unconfounded, the crude RD equals the causal risk difference pr[Y_a0=1 = 1] - pr[Y_a0=0 = 1] in the source population. The causal risk difference is the average of the individual causal risk differences Y_a0=1 - Y_a0=0. Similarly, the crude RR equals the causal RR, pr(Y_a0=1 = 1)/pr(Y_a0=0 = 1), and the crude OR equals the causal OR, pr(Y_a0=1 = 1)pr(Y_a0=0 = 0)/{pr(Y_a0=1 = 0)pr(Y_a0=0 = 1)}. Because of the possibility of effect modification, the population causal parameter need not equal the causal parameter within a stratum of the measured risk factors L₀ even if treatment is unconfounded. Effect modification is considered in section 9.

3. Models for Point-Treatment Studies

The causal RD, RR, and OR can also be expressed in terms of the parameters of the following linear, log linear, and linear logistic models for the two counterfactual probabilities pr(Y_a0=1 = 1) and pr(Y_a0=0 = 1). MATH MATH MATH where Y_a0 is Y_a0=1 if a₀ = 1 and Y_a0 is Y_a0=0 if a₀ = 0. Specifically, the causal RD = ψ₁, the causal RR = e^θ1, and the causal OR = e^β1. Models 1-3 are saturated MSMs. They are marginal models, because they model the marginal distribution of the counterfactual random variables Y_a0=1 and Y_a0=0 rather than the joint distribution (that is, models 1-3 do not model the correlation of Y_a0=1 and Y_a0=0). They are structural models, because they model the probabilities of counterfactual variables and in the econometric and social science literature models for counterfactual variables are often referred to as structural. ^8,12 Finally, they are saturated, because each has two unknown parameters and thus each model places no restriction on the possible values of the two unknown probabilities pr(Y_a0=1 = 1) and pr(Y_a0=0 = 1). Note that these models do not include covariates, because they are, by definition, models for causal effects on the entire source population; they are not models for observed associations.

Equation U1

Equation U2

Equation U3

The crude RD, RR, and OR can also be expressed in terms of the parameters of the following saturated linear, log linear, and linear logistic models for the observed outcome Y. MATH MATH MATH These are models for associations observed when comparing subpopulations (defined by levels of treatment) of the source population. The crude RD equals ψ^′₁, the crude RR equals e^θ′1, and the crude OR equals e^β′1. The parameters of the associational models 4-6 will differ from the parameters of the MSMs 1-3, except when treatment is unconfounded. Because models 4-6 are models for the observed data, (asymptotically) unbiased estimates of the model parameters can be obtained using standard statistical software (assuming no selection bias or measurement error). When treatment is unconfounded, these same estimates will also be unbiased for the corresponding causal parameters of models 1-3. For example, to fit models 4-6, one could use the general-purpose SAS program Proc Genmod, using the model statement Y = A₀ with the outcome Y specified as a binomial variable. To estimate ψ^′₁, one would specify the identity link; to estimate θ^′₁, the log link; and for β^′₁, the logit link. Programs analogous to Proc Genmod also exist in other packages, such as S-Plus, Gauss, and Stata. ^13-16

Equation U4

Equation U5

Equation U6

4. No Unmeasured Confounders

Suppose now that treatment is confounded. Then the crude association parameter will not equal the corresponding causal parameter. Similarly, the parameters of the MSMs will fail to equal the parameters of the corresponding observed data models (for example, β₀ ≠ β^′₀ and β₁ ≠ β^′₁). Assuming we have no unmeasured confounders given data on measured confounders L₀, unbiased estimates of the causal parameters ψ₁, θ₁, and β₁ can, however, still be obtained using Proc Genmod by performing a weighted analysis. Specifically, using the weight statement (that is, option SCWGT) in Proc Genmod, each subject i is assigned a weight w_i equal to the inverse of the conditional probability of receiving his or her own treatment. That is, w_i = 1/pr[A₀ = a_0i|L₀ = l_0i], where, for example, l_0i is the observed value of the variable L₀ for subject i. The true weights w_i are unknown but can be estimated from the data in a preliminary logistic regression of A₀ on L₀. For example, we might specify the logistic regression model MATH where A₀ is AZT treatment, L₀ is, for example, the column vector of covariates with components age, CD4 count, WBC count, hematocrit, and presence of symptoms, and α₁ is a row vector of unknown parameters. We can then obtain estimates (α̂₀, α̂₁) of (α₀, α₁) using standard logistic regression software. Then, for a subject i with A₀ = 0 and L₀ = l_0i, we would estimate w_i = 1/pr[A₀ = 0|L₀ = l_0i] by 1/{1/[1 + exp(α̂₀ + α̂₁l_0i)]} = 1 + exp(α̂₀ + α̂₁l_0i). For a subject with A₀ = 1 and L₀ = l_0i, we would estimate w_i = 1/pr[A₀ = 1|L₀ = l_0i] by {1 + exp(α̂₀ + α̂₁l_0i)}/exp(α̂₀ + α̂₁l_0i) = 1 + exp(-α̂₀ - α̂₁l_0i).

Equation U7

In summary, if there are no unmeasured confounders given data on L₀, one can control confounding (due to L₀) by modifying the crude analysis by weighting each subject i by w_i. The denominator of w_i is informally the probability that a subject had his or her own observed treatment. Thus, we refer to these weighted estimators as IPTW estimators.

Why does this approach work? The effect of weighting in Proc Genmod is to create a pseudopopulation consisting of w_i copies of each subject i. That is, if, for a given subject, w_i = 4, the subject contributes four copies of him- or herself to the pseudopopulation. This new pseudopopulation has the following two important properties. First, in the pseudopopulation, unlike the actual population, A₀ is unconfounded by the measured covariates L₀. Second, pr(Y_a0=1 = 1) and pr(Y_a0=0 = 1) in the pseudopopulation are the same as in the true study population so that the causal RD, RR, and OR are the same in both populations. Hence, it follows that we can unbiasedly estimate the causal RD, RR, and OR by a standard crude analysis in the pseudopopulation. But this is exactly what our IPTW estimator does, because the weights w_i serve to create, as required, w_i copies of each subject. In the Appendix, we present a detailed numerical example to clarify further our IPTW methodology, and we compare our methodology with the propensity score methodology of Rosenbaum and Rubin. ¹⁷

5. Unmeasured Confounding

In the presence of unmeasured confounding factors U₀, one could still unbiasedly estimate the causal risk difference, risk ratio, and odds ratio as above if one used weights w_i = 1/pr[A₀ = a_0i|L₀ = l_0i, U₀ = u_0i] in implementing the analysis in Proc Genmod. Nevertheless, because data on U₀ are not observed, it is not possible to estimate these weights unbiasedly. Indeed, unbiased estimation by any method is impossible in the presence of unmeasured confounding factors without strong additional assumptions.

6. Multilevel Treatment and Unsaturated MSMs

Suppose again that treatment is unconfounded (as represented in Figure 2, c) but now A₀ is an ordinal variable representing a subject's daily dose in units of 100 mg of AZT. Possible values of A₀ are 0, 1, …, 14, 15. In that case, the number of potential outcomes associated with each subject will be 16. Specifically, we let Y_a0 be the value of Y that would have been observed had the subject received dose a₀ rather than the observed dose. Thus, in principle, a subject has a separate counterfactual variable for each of the 16 possible AZT doses a₀. The subject's observed outcome Y is the outcome Y_a0 corresponding to the dose a₀ equal to the subject's observed dose. For expositional convenience, we continue to refer to all of the Y_a0's as counterfactuals, even though for a₀ equal to the observed dose, Y_a0 is the factual variable Y. Because there are so many potential outcome variables Y_a0, we can no longer conveniently perform a saturated analysis. Rather, we would usually assume a parsimonious dose-response relationship by specifying a linear logistic MSM such as MATH This model says that the probability of success had all subjects been treated with dose a₀ is a linear logistic function of the dose with slope parameter β₁ and intercept β₀, so e^β1 is the causal OR associated with an increase in AZT dose of 100 mg.

Equation U8

We contrast the MSM model 8 with the following linear logistic association model for the observed data. MATH Assuming no selection bias or measurement error, we can unbiasedly estimate the associational parameters β^′₀ and β^′₁ by fitting the linear logistic model 9 using a standard logistic regression software package such as SAS Proc Logistic or Proc Genmod. If the treatment is unconfounded, then the parameters of models 8 and 9 are equal. As a consequence, our logistic regression estimate of β^′₁ is also an unbiased estimate of our causal parameter β₁.

Equation U9

If treatment is confounded by the measured variables L₀, then β₁ ≠ β^′₁ and our standard logistic regression estimate of β^′₁ is a biased estimate of the causal parameter β₁ owing to confounding by L₀. However, even when treatment is confounded, if there are no unmeasured confounders given L₀, then one can obtain unbiased estimates of the causal parameter β₁ of model 8 by fitting the logistic model 9 with Proc Genmod if one uses subject-specific weights w_i = 1/pr(A₀ = a_0i|L₀ = l_0i). Again, in practice, w_i is unknown and one must estimate it from the data by specifying a model.

For example, one might specify the following polytomous logistic model:MATH which can be fit in SAS using Proc Logistic or Genmod to obtain estimates of the parameters α₀₁, α₀₂, …, α₀₁₅, and α₁.

Equation U10

6.

1. Stabilized Weights

The probabilities pr[A₀ = a_0i|L₀ = l_0i] may vary greatly between subjects when components of L₀ are strongly associated with A₀. This variability can result in extremely large values of the weight w_i for a few subjects. These few subjects will contribute a very large number of copies of themselves to the pseudopopulation and thus will dominate the weighted analysis, with the result that our IPTW estimator will have a large variance and will fail to be approximately normally distributed. If the MSM model is saturated (for example, models 1-3), this variability is unavoidable, because it reflects a lack of information in the data as a result of the confounders L₀ being highly correlated with treatment A₀. For unsaturated MSMs, such as model 8, however, this variability can, to a considerable extent, be mitigated by replacing the weight w_i by the stabilized weight sw_i = pr[A₀ = a_0i]/pr[A₀ = a_0i|L₀ = L_0i]. To understand the stabilized weight, suppose A₀ was unconfounded so that A₀ and L₀ are unassociated and pr[A₀ = a_0i] = pr[A₀ = a_0i|L₀ = l_0i]. Then sw_i = 1, and each subject contributes the same weight. When A₀ is confounded, sw_i will not be constant but will vary around the number 1, depending on the subject's value of L₀. sw_i, however, will still tend to be much less variable than w_i. Furthermore, Robins ^1,2 shows that, when we use the weight sw_i rather than the weight w_i in Proc Genmod, the estimates of the parameters β of an MSM remain unbiased and, in the case of an unsaturated MSM, will generally be less variable. For saturated MSMs, the variability of our estimate of β will be the same whether we use the stabilized or unstabilized weights.

Of course, pr[A₀ = a_0i] and pr[A₀ = a_0i|L₀ = l_0i] are unknown and must be estimated. pr[A₀ = a_0i|L₀ = l_0i] can be estimated as described above; pr[A₀ = a_0i] can be estimated as the proportion of subjects in the study sample with A₀ equal to a_0i. This estimate is equivalent to that obtained by fitting the polytomous model MATH where we place an asterisk on the parameter α^*_0a0 to indicate that this parameter will differ from the parameter α_0a0 of model 10 when A₀ is confounded.

Equation U11

In Ref 2, Robins introduces augmented IPTW estimators. These estimators are even more efficient than the IPTW estimator that uses stabilized weights but are more difficult to compute.

6.

2. Continuous Treatment

Suppose we were able to measure a subject's daily intake of AZT to the nearest tenth of a milligram, so that now AZT is essentially a continuous treatment. Further, for expositional convenience, assume that no subject has AZT dose near 0, and that we can effectively model the distribution of A₀ as normal. Now each individual has an extremely large number of counterfactual outcomes Y_a0. One can still obtain unbiased estimates of the causal parameter β₁ of model 8 by fitting the logistic model 9 with Proc Genmod if one uses the stabilized weights sw_i = f(a_0i)/f(a_0i|l_0i), where f(a₀|l₀) is the conditional density of the continuous variable A₀ given L₀, and f(a₀) is the marginal density of the continuous variable A₀. To estimate f(a₀|l₀), one might specify that, given L₀, A₀ is normal with mean α₀ + α₁L₀ and variance ς². Then unbiased estimates (α̂₀, α̂₁, μ̂²) of (α₀, α₁, ς²) can be obtained by ordinary least-squares regression of A₀ on L₀ using, for example, Proc REG in SAS. Then f(a_0i|l_0i) would be estimated by the normal density (2πμ̂²)^-1/2 exp{-[a_0i - (α̂₀ + α̂₁l_0i)]²/2μ̂²}. To estimate the numerator f(a_0i) of the stabilized weight sw_i, one might specify that A₀ is normal with mean α^*₀ and variance ς*². f(a_0i) could be estimated by the normal density (2πμ̂*²)^-1/2 exp[-(a_0i - α̂^*₀)²/2μ̂*²] where α̂^*₀ is the average of the observed A₀s and μ̂*² is their empirical variance. When A₀ is continuous, estimates based on the unstabilized weights w_i = 1/f(a_0i|l_0i) have infinite variance and thus cannot be used. ^1,2

6.

3. Confidence Intervals

As described above, we shall estimate the parameters of the MSM 8 by fitting the association model 9 in Proc Genmod using estimates of the stabilized weights sw_i. If we choose the option repeated and specify an independence working correlation matrix, the Proc Genmod program will also output a 95% robust Wald confidence interval for β₁ given by β̂₁ ± 1.96 MATH, where var (β̂₁) is the so-called robust¹⁸ or sandwich estimator of the variance of β̂₁. Robins ^1,2 shows that the robust Wald interval will have coverage probability of at least 95%, although narrower valid intervals can be obtained with some additional programming. ^1,2 The ordinary nonrobust model-based Wald confidence interval outputted by most weighted logistic regression programs will not be guaranteed to provide at least 95% coverage and thus should be avoided. Other software packages such as S-Plus, Gauss, and STATA also offer robust variance estimators and could be used in place of Proc Genmod. ^13,19-21

Equation U12

7. Time-Dependent Treatments

We now return to the setting of section 1, in which A_k is the dose of treatment AZT on the k^th day from start of follow-up and Y is a dichotomous outcome variable measured at end of follow-up on day K + 1. Similarly, L_k represents the value on day k of the vector of all measured risk factors for the outcome. Let Ā_k = (A₀, A₁, …, A_k) be the treatment or exposure history through day k and let Ā = Ā_K. Define ¯L_k and ¯L similarly. Let Y_ā be the value of Y that would have been observed had all subjects received dose history ā = (a₀, a₁, …, a_K) rather than their observed dose history Ā. Note that, even if a_k is dichotomous on each day k (that is, on each day a subject is either on or off treatment), there will be 2^K dose histories ā and thus 2^K possible counterfactuals, only one of which is observed (that is, is factual). Thus, it may not be possible to estimate a saturated MSM. Therefore, we would generally assume some sort of parsimonious dose-response relationship by specifying a linear logistic MSM such as MATH where cum (ā) = ∑_k=0^K a_k is the cumulative dose through end-of-follow-up associated with the dose history ā.

Equation U13

We contrast the MSM 12 with the following linear logistic association model for the observed data:MATH Assuming no loss to follow-up selection bias or measurement error, we can unbiasedly estimate the parameters β^′₁ by fitting the linear logistic model 13 using a standard logistic regression software package. If the treatment is unconfounded, the parameters of models 12 and 13 are equal. As a consequence, our logistic regression estimate of β^′₁ is also an unbiased estimate of our causal parameter β₁. If the treatment is confounded, then β₁ ≠ β^′₁ and our standard logistic regression estimate of β^′₁ is a biased estimate of the causal parameter β₁ as a result of confounding by ¯L_k. When treatment is confounded, however, if there is no unmeasured confounder given the L_k, then one can still obtain unbiased estimates of the causal parameter β₁ of model 12 by fitting the logistic model 13 with the stabilized weights MATH where ∏_k=0^K b_k = b₀ × b₁ × b₂ × … × b_K and Ā_-1 is defined to be 0. Note in the special case in which K = 0 (that is, a point-treatment study), models 12 and 13 reduce to our previous models 8 and 9 and 14 reduces to our previous sw_i. The denominator of sw_i is informally the conditional probability that a subject had his or her own observed treatment history through time K. With time-dependent treatments the variation in the unstabilized weights will often be enormous, with the result that the resulting estimator of β can be highly variable with a markedly non-normal sampling distribution. We therefore strongly recommend the use of stabilized weights.

Equation U14

Equation U15

We emphasize that when treatment is confounded, it is the parameter β₁ of our MSM 12, as opposed to the parameter β^′₁ of the association model 13 that is of policy importance. To see why, consider a new subject from the source or target population exchangeable with the N study subjects. We would like to administer to the subject the treatment ā that minimized probability that he or she has detectable HIV in the serum at the end of follow-up, that is, pr(Y_ā = 1). Thus, for example, if the parameter β₁ of MSM 12 is positive (that is, the probability of having HIV in the blood increases with increasing duration of AZT treatment), we would withhold AZT treatment from our subject. In contrast, the parameter β^′₁ of model 13 may be confounded by the association of covariates with treatment. For example, suppose physicians preferentially started AZT on subjects who, as indicated by their prognostic factor history (for example, CD4 count), were doing poorly. Further, suppose that AZT had no causal effect on Y (that is, β₁ = 0). Then the parameter β^′₁ (and thus our estimate from the unweighted logistic regression) will be positive but will have no causal interpretation as the effect of AZT on Y.

7.

1. Bias Induced by Controlling for a Variable Affected by Treatment

One might suppose that an alternative approach to controlling confounding by measured covariates is an unweighted logistic model that adjusts for confounder history ¯L ≡ ¯L_K, such as MATH where, for simplicity, we here assume L_k consists of a single covariate CD4 count at time k. Nevertheless, even under our assumption of no unmeasured confounders, the parameter β^″₁ differs from the parameter β₁ of our MSM. What is worse is that the parameter β^″₁ will generally not have a causal interpretation, even if the model for pr[Y = 1|Ā = ā, ¯L = ¯l] is correctly specified. This is because cum (Ā) depends on a subject's entire treatment history, including A₀, and A₀ may affect the time-dependent covariates L_k and L_k-1. Fitting a logistic model that adjusts for a covariate that is both affected by treatment and is a risk factor for the outcome provides an unbiased estimate of the association parameter β^″₁ but a biased estimate of the causal parameter β₁. This is true even under the null hypothesis of no direct, indirect, or overall treatment effect (so that β₁ of model 12 equals 0) when, as in Figure 1, a component of L_k (for example, red blood count) and the outcome Y have an unmeasured common cause U₀ (for example, the baseline number of bone marrow stem cells). ^5,7,22-26

Equation U16

To summarize, standard regression methods adjust for covariates by including them in the model as regressors. These standard methods may fail to adjust appropriately for confounding due to measured confounders L_k when treatment is time varying, because (1) L_k may be a confounder for later treatment and thus must be adjusted for, but (2) may also be affected by earlier treatment and thus should not be adjusted for by standard methods. A solution to this conundrum is to adjust for the time-dependent covariates L_k by using them to calculate the weights sw_i rather than by adding the L_k to the regression model as regressors.

8. Estimation of the Weights

We now describe how to estimate the weights sw_i. For simplicity, we again assume that the treatment A_k at each time k is dichotomous. Consider first the denominator of model 14. We begin by estimating the unknown probability pr[A_k = 1|Ā_k = 1|Ā_k-1 = ā_k-1, ¯L_k = ¯l_k] using a pooled logistic model that treats each person-day as an observation. For example, we might fit the model MATH where, for example, l_k is the vector of CD4 count, WBC, hematocrit, and an indicator for symptoms at time k, and the α₄, α₅, α₆, and α₇ are row vectors. This model says that the probability of being treated on day k depends in a linear logistic fashion on the day k, the previous 2 days' treatment, the current and previous days' covariates, an interaction between yesterday's treatments and today's covariates, and the baseline covariates.

Equation 17

One can fit model 15 using any standard logistic regression program. The numerator probabilities can be estimated similarly, except that, in fitting model 15, we remove the last four terms that are functions of the covariates. That is, we fit the model as follows:MATH

Equation U18

For each subject i, we then have our logistic program output the estimated predicted values p̂_0i, …, p̂_Ki from the fit of model 15, which are maximum likelihood estimates of pr[A_k = 1|Ā_k-1 = ā_(k-1)i, ¯L_k = ¯l_ki]. Similarly, we have outputted the predicted values p̂^*_1i, …, p̂^*_Ki from model 16, which are estimates of the quantities pr[A_k = 1|Ā_k-1 = ā_(k-1)i]. Then we estimate sw_i by MATH For example, 1 - p̂^*_ki is an estimate of the probability pr[A_k = a_ki|Ā_k-1 = ā_(k-1)i] when a_ki = 0. The data analyst will need to write a small program to compute sw_i for each subject from the predicted values outputted from the fit of models such as 15 and 16.

Equation U19

Under our assumption of no unmeasured confounders, the resulting estimate of the causal parameter β₁ will be unbiased, provided the model 15 for pr[A_k = 1|Ā_k-1 = ā_k-1, ¯L_k = ¯l_k] is correctly specified. Furthermore, under these same conditions, the 95% robust Wald confidence interval will be guaranteed to cover β₁ at least 95% of the time. The estimate of β₁ will remain unbiased even if the model 16 for pr[A_k = 1|Ā_k-1 = ā_k-1] is misspecified. ^1,2 Indeed, if model 15 is correct and treatment is confounded, model 16 is guaranteed to be somewhat misspecified because of the noncollapsibility of logistic models. ²⁵

9. Effect Modification by Pretreatment Covariates

MSMs can be generalized to allow one to include pretreatment covariates. For example, model 12 could be generalized to MATH where V is a component of the vector of measured pretreatment covariates L₀, and β₃ denotes a treatment-covariate interaction. Note in model 18, β₁ + β₃v represents the effect of cumulative treatment on a linear logistic scale within level v of the baseline variable V. As our IPTW estimators already automatically adjust for any confounding due to V, the particular subset V of L₀ that an investigator chooses to include in model 18 should only reflect the investigator's substantive interest. For example, a variable V should be included in model 18 only if the investigator both believes that V may be an effect modifier and has greater substantive interest in the causal effect of treatment within levels of the covariate V than in the source population as a whole.

Equation U20

We obtain unbiased estimates of the parameters of model 18 under the assumption of no unmeasured confounders by fitting an association model such as MATH using Proc Genmod with the estimated weights sw_i of Eq 17, modified only in that p^*_k is now the estimated predicted value from the fit of a model such as MATH Elementary epidemiologic textbooks emphasize that effect modification is logically distinct from confounding. Nonetheless, many students have difficulty understanding the distinction, because the same statistical methods (stratification and regression adjustments) are used both for confounder control and detection of effect modification. Thus, there may be some advantage to teaching elementary epidemiologic methods using marginal structural models, because then methods for confounder control (inverse-probability-of-treatment weighting) are distinct from methods for detection of effect modification (adding treatment covariate interaction terms to an MSM).

Equation U21

Equation U22

Finally, an important caveat: MSMs cannot be used to model the interaction of treatment with a time-varying covariate. For this, structural nested models should be used. ^4-6 Therefore, it is not valid to include in the covariate V in model 18 any components of the time-dependent covariate L_k measured at any time k > 0.

10. Censoring by Loss to Follow-Up

Heretofore, we have assumed that each study subject is observed until end of follow-up at time K + 1. In this section, we allow for censoring by loss to follow-up. Specifically, let C_k = 1 if a subject was lost to follow-up by day k and C_k = 0 otherwise. We assume that once a subject is lost to follow-up, the subject does not later re-enter follow-up. No new idea is required to account for censoring, provided we conceptualize censoring as just another time-varying treatment. From this point of view, to want to adjust for censoring is only to say that our interest is in the causal effect of the treatment Ā if, contrary to fact, all subjects had remained uncensored, rather than having followed their observed censoring history. Our goal remains to estimate the parameter β₁ of the logistic MSM 12 except now Y_ā refers to a subject's outcome if, possibly contrary to fact, the subject has followed treatment history ā and has never been censored. Again, we can do so if there are no unmeasured confounders for both treatment and censoring. To formalize this idea, one adds at each time k the variable C_k to the graph in Figure 1 just before L_k and after A_k-1. Then, the assumption of no unmeasured confounders for treatment and censoring is that no arrow arising from the unmeasured causal risk factors U goes directly into either C_k or A_k for any k. In that case, the measured covariates L_k are sufficient to adjust both for confounding and selection bias due to loss to follow-up.

Again, we can obtain unbiased estimates of the causal parameters β₁ by fitting the linear logistic association model 13 with appropriate weights included. Because the outcome Y is unobserved unless the subject does not drop out, that is, ¯C = (C₀, …, C_K+1) = 0, our weighted logistic regression fit of model 13 is restricted to uncensored subjects. The required subject-specific weight is sw_i × sw_i^†, where MATH and, in addition, in defining and estimating sw_i, we now add to the right side of each conditioning event in models 14-16 the event ¯C_k = 0, because otherwise, A_k would not be observed. The unknown probabilities in sw_i^† can be estimated using a pooled logistic model that treats each person-day as an observation. Specifically, we fit analogs of models 15 and 16 for logit pr[C_k = 0|¯C_k-1 = 0, Ā_k-1 = ā_k-1, ¯L_k = ¯l_k] and for logit pr[C_k = 0|¯C_k-1 = 0, Ā_k-1 = ā_k-1]. Note that the denominator of the product sw_i × sw_i^† is informally the conditional probability that an uncensored subject had his or her observed treatment and censoring history through time K + 1. Thus, we refer to our weighted logistic estimator as an inverse-probability-of-treatment-and-censoring weighted estimator. If we view (A_k, C_k) as a joint treatment at time k, then one can informally interpret this denominator as simply the probability that a subject follows his or her own treatment history, which is exactly the interpretation that we had previously in the absence of censoring.

Equation U23

11. Limitations of Marginal Structural Models

It is shown in Ref 2 and Appendix 2 that our IPTW estimators will be biased and thus MSMs should not be used in studies in which at each time k there is a covariate level l_k such that all subjects with that level of the covariate are certain to receive the identical treatment a_k. For example, this circumstance implies that MSMs should not be used in occupational cohort studies. To see why, consider an occupational cohort study in which A_k is the level of exposure to an industrial chemical at time k and L_k = 1 if a subject is off work at time k and L_k = 0 otherwise. Then all subjects with L_k = 1 have A_k = 0, because all subjects off work are unexposed. Similarly, in a study of the effect of screening on mortality from cervical cancer, women who have had their cervix operatively removed by time k (which we denote by L_k = 0) cannot receive exposure (that is, screening) at that time, so MSMs should not be used. Nevertheless, g-estimation of structural nested models can always be used to estimate exposure effects, even in studies in which MSMs cannot be used. In many studies, such as the analysis of the Multicenter AIDS Cohort Study data described in our companion paper, ²⁷ we believe, based on substantive considerations, the above difficulty does not occur and MSMs are a practical method.

12. Conclusion

We have described how to use MSMs to estimate the causal effect of a time-varying exposure or treatment on a dichotomous outcome. In our companion paper, ²⁷ we extend our results to survival time outcomes and compare and contrast methods based on MSMs to alternative, previously proposed methods, based on g-estimation of structural nested models and on estimation of the g-computation algorithm formula.