History-adjusted Marginal Structural Models for Estimating Time-varying Effect Modification

Marginal structural models

Marginal structural models (MSMs), introduced by Robins et al. (1–3), represent a major advance in statistical methodology for causal inference. In epidemiologic and clinical research, longitudinal data on treatment and covariates are often collected from cohorts of subjects. When treatment status can change over time, conventional analytic approaches (such as multivariable regression) often fail to provide valid causal inferences about treatment effects. MSMs address this well-recognized problem of time-dependent confounding. As a result, in many research contexts, MSMs are the best available analytic method.

MSMs have traditionally been restricted to the estimation of the causal effects of treatment or exposure conditional on baseline covariates (2). Thus, it has been possible to use this method to address questions such as “What is the effect of a treatment and how does it differ between study members with different covariate values at a given time point, such as entry to the study?” but not “How does the effect of a treatment differ as a result of changing values of a covariate over the course of the study?”. In other words, it has not been possible to use MSMs to estimate modification of causal effects by time-dependent covariates. In many research settings, estimation of time-dependent effect modification is of major interest. In clinical settings in particular, treatment decisions are often modified over time as a result of changing values for a patient’s covariates. Estimation of how past values of these covariates modify the future causal effect of a treatment has important implications for understanding the mechanistic action of the treatment, as well as for guiding clinical decision-making.

In this paper, we provide a practical introduction to history-adjusted marginal structural models (HA-MSMs), a type of MSM that allows estimation of effect modification by time-varying covariates. HA-MSMs are based on the concept, first introduced by Joffe et al. (4, 5) and implemented by Feldman et al. (6), of assuming a standard MSM at each time point during a study. Covariates and treatment history up to that time point are treated as baseline data, and the causal effect of treatment after that time point is estimated. Notably, HA-MSMs can assume a common model across time points, rather than simply assuming a separate MSM at each time point. HA-MSMs provide an alternative to structural nested models, an approach proposed by Robins (7) for the estimation of effect modification by time-dependent covariates.

Van der Laan et al. (8) recently established the formal statistical theory behind HA-MSMs. In this article, we use a practical application to provide a nontechnical introduction to HA-MSMs, including the underlying assumptions, their implementation using standard software, and the interpretation of results. Using observational cohort data from human immunodeficiency virus (HIV)-infected persons with detectable viremia on antiretroviral therapy, we estimate the effect of time to treatment modification on future CD4-positive T-lymphocyte (CD4 T-cell) count and how this effect differs given current CD4 T-cell count.

Antiretroviral therapy: when to switch?

HIV evolves rapidly in the presence of selective pressure. This leads to the accumulation of mutations that confer “resistance” to antiretroviral drugs. Therefore, the optimal manner in which to avoid the rapid emergence of resistance-associated mutations is to completely suppress viral replication. This can be achieved in many patients with standard three-drug combination regimens (9). However, a substantial proportion of treated patients fail to achieve complete viral suppression. Such patients are often switched to a new treatment regimen, but this can lead to the use of all available therapeutic options. Since many patients with incomplete viral suppression continue to do well immunologically (and therefore clinically) (10–12), some clinicians choose not to switch them to a new regimen as long as their CD4 T-cell counts remain elevated. Hence, clinicians are often faced with a dilemma in patients with detectable viremia who are on antiretroviral therapy: Should therapy be switched as soon as possible, thereby risking using up all available drugs quickly and exposing patients to increasingly complicated and potentially toxic regimens? Or should patients be maintained on a partially suppressive regimen as long as they are doing well immunologically and clinically, even though this approach will allow the ongoing accumulation of drug-resistance mutations that can limit future therapeutic options? (For a review of this issue, see the paper by Deeks (13).)

In this article, we estimate the effect of nonsuppressive therapy on future CD4 T-cell count and estimate how this outcome differs depending on a patient’s current CD4 T-cell count and time since virologic failure. Data for these analyses were drawn from the Study on the Consequences of the Protease Inhibitor Era (SCOPE), an observational study of a cohort of HIV-infected patients in San Francisco, California. Participants in SCOPE are seen at 4-month intervals. At each study visit, they complete interviewer- and self-administered questionnaires examining numerous domains, including socioeconomic status (housing, income, employment), antiretroviral medication use and adherence, occurrence of opportunistic infection or malignancy, and recreational drug use. Plasma HIV RNA levels and CD4/CD8 T-cell counts are measured at each visit, as well as between visits, according to clinician discretion. Importantly, decisions as to when and how to modify therapy are made by primary care providers.

We retrospectively identified all subjects who experienced virologic failure while being observed in SCOPE. Subjects became eligible for our analyses (t = 0) if they failed to achieve an undetectable HIV RNA level (<75 copies/ml) by week 24 on a new regimen or if they rebounded from an undetectable level. The exposure of interest was time to modification of the antiretroviral regimen, defined as switching or interrupting the use of at least one drug. This exposure was summarized as a binary variable at each time point, indicating whether or not a subject had switched from his original non-suppressive antiretroviral therapy regimen. Subjects were only allowed to switch once in our analyses. The method can be easily extended, however, to encompass more complex treatment patterns.

Below, we rely on this data structure to illustrate our method. We then present the results of our analyses and discuss their clinical significance.

The counterfactual framework

The causal effect of a treatment on an individual can be defined as the difference between that person’s outcomes with and without the treatment. Such outcomes are termed “counterfactual,” because only one outcome is observed for each individual. MSMs are models of how the population distribution of these counterfactual outcomes changes as a result of changes in treatment.

We begin by introducing some standard notation. Treatment over the course of the study (t = 0, …, K) is denoted Ā(K) = (A(0), …, A(K)), and covariates are denoted L̄ (K + 1), where treatment occurs after covariates are measured at a given time point and K + 1 is the end of follow-up. In our HIV example, Ā(K) is a vector of binary variables, equal to 1 for each time point up to the point at which a subject switches therapy and 0 thereafter. L̄_ā (K + 1) denotes the counterfactual CD4 T-cell count and other covariates that would have been observed over time if the subject had switched therapy at the time implied by Ā = ā. The set of these counterfactuals under each possible treatment history (switch time) is the full data set, X = (L̄_ā (K + 1), ā ∈ A), where A denotes the set of possible switch times. The outcome for a given time point t is the counterfactual CD4 T-cell count measured 8 months in the future (m = 8) under the switching time indicated by ā, denoted Y_ā (t + m).

If we observed the counterfactual CD4 T-cell counts for each subject under each possible switch time, we could estimate the causal effect of waiting to switch therapy by simply comparing the counterfactual outcomes under different switch times. However, we only observe the CD4 T-cell counts for each person under a single (nonrandom) switch time. As a result, in order to estimate the effect of time to switching treatments on CD4 T-cell count using the observed data, we must assume that the covariates we measured are sufficient to control for confounding. For example, within strata defined by our measured confounders, there must be no unmeasured variables that predict, at any time point, both the probability of switching treatment and CD4 T-cell count 8 months in the future. More formally, this assumption is stated as the sequential randomization assumption, which says that the treatment assignment at each time point depends only on the observed past:

where g(A(t)|X) = P(A(t)|X). Under the sequential randomization assumption, the treatment mechanism, which describes how treatment is assigned at each time point, can be written

In addition, we assume experimental treatment assignment. Experimental treatment assignment states that treatment is not assigned deterministically on the basis of a patient’s covariates. For instance, regardless of her covariate values, a patient who has not yet switched treatments must have a positive probability of both switching and not switching treatments at each time point.

Under these assumptions, we could use standard MSMs to ask, for example, “At baseline (when virologic failure occurs), how does time until switching to a new regimen affect CD4 T-cell count 8 months later? How does this effect differ depending on a patient’s CD4 T-cell count at baseline?” However, an MSM assumed at a single time point does not allow us to estimate how the effect of future non-suppressive therapy may change as a result of changes in a patient’s CD4 T-cell count. HA-MSMs directly address this question.

History-adjusted MSMs

HA-MSMs rely on the identical causal framework as standard MSMs but estimate a different parameter of interest. HA-MSMs assume a standard MSM at each time point during the study and model the counterfactual outcome indexed by treatment that occurs after that time point, conditional on some subset of the observed history up to that time point. Importantly, HA-MSMs allow the use of a common model across time points. In other words, HA-MSMs model some parameter of the counterfactual outcome if the study population were to follow their observed treatment history up to time j, followed by a specified counterfactual future treatment history until the outcome is measured, conditional on a subset of (possibly time-varying) covariates and/or treatment history measured before time j. In this article, we focus on HA-MSMs concerned with the mean of these counterfactual outcomes; however, the same framework can be readily adopted to model any other parameter.

We denote a future longitudinal treatment regimen, beginning at time j and continuing until the outcome is measured m time points later, as a(j, j + m − 1) ≡ (a(j), a(j + 1), …a(j + m − 1)), for j = 0, …, K + 1 − m. The effect modifiers of interest are denoted V(j) ⊂ (L̄ (j), Ā (j − 1)), a subset of a subject’s treatment and covariate history up to time j. For each time point in the study for which the outcome m time points later is defined, HA-MSMs model the expectation of the counterfactual outcome Y_{Ā(j−1),a(j, j+m−1)}(j + m), conditional on V(j), under each possible future treatment regimen. Thus, HA-MSMs are concerned with estimation of the following parameter: E(Y_{Ā(j−1),a(j, j+m−1)}(j + m)|V(j)).

Applied to our example, future antiretroviral treatment from time j until the outcome is measured, denoted a(j, j + m − 1), consists of a vector of counterfactual treatment decisions a(j), …a(j + m − 1), where a(t) = 0 if a subject is assigned to switch treatment at or before time t; otherwise, a(t) = 1. This vector of future treatment decisions exists for each subject beginning at each time point j = 0, …, K + 1 − m. We summarize a(j, j + m − 1) as $c(j)={\sum }_{l=j}^{j+m-1}a(l)$, which represents the future time (beginning at time point j) that the subject will spend on his/her original failing therapy before either treatment is switched or the outcome is measured. The current CD4 T-cell count at time j is denoted CD4(j), a subset of the full covariate history measured over time, L̄ (j). For each time point j, we are interested in the mean counterfactual CD4 T-cell count that would exist 8 months later (m = 8) among persons who had not yet switched therapy if they were to switch therapy at a specified counterfactual time after j.

To address this question, we might assume the following model:

In other words, we might assume that, among persons who have not yet switched treatment (Ā (j − 1) = 1), the counterfactual CD4 T-cell count 8 months later depends on additional time until switching (c(j)), but the magnitude of this effect differs depending on the amount of time a patient has already spent on nonsuppressive therapy (j) and current CD4 T-cell count (CD4(j)).

This model allows us to estimate the effect of each additional month until treatment modification on future CD4 T-cell count, among patients who have been on their current nonsuppressive therapy for different durations and have different current CD4 T-cell counts. For example, by testing whether β₄ = 0, we are testing the hypothesis that a subject’s current CD4 T-cell count modifies the effect of future time until switching.

Inverse-probability-of-treatment-weighted estimation

Several HA-MSM estimators are available (8); here, we focus on the inverse-probability-of-treatment-weighted (IPTW) estimator, which can be implemented with standard software using weighted least squares. For each time point j in the study, each subject receives a weight which is, informally, the inverse of that subject’s probability of receiving the treatment that she actually received, from time point j until the outcome is measured. If a subject has a longitudinal treatment regimen beginning at time point j that occurs frequently in the data among subjects with her covariate and treatment history, she receives a small j-specific weight. In contrast, if the subject has an unusual longitudinal treatment regimen given her covariates, the subject will receive a large weight. For example, patients whose CD4 T-cell counts have recently declined may be more likely to switch therapy. In this case, a subject who did not switch therapy despite a recent decline in CD4 T-cell count would receive a large weight.

The first step in implementing the IPTW estimator is to model the treatment mechanism or, in other words, to fit a predictive model of treatment at each time point t, given the observed past up to that time point: g(A(t)|Ā (t − 1), L̄ (t)), t = 0, …, K. For example, we can model the treatment decision (to switch therapy or not) made at every time point using logistic regression. Recall that once a subject switches therapy, he is no longer at risk of switching in the future. Thus, when fitting our model of the probability of staying on a certain regimen at a given time point (A(t) = 1), we fit the model only among subjects who have not already switched before that time point (A(t − 1) = 1). Note that, for the IPTW estimator to be consistent, the estimate of the treatment mechanism must be consistent and the covariates included in the model must be sufficient to control for confounding.

For each time point j = 0, …, K + 1 − m, the model of the treatment mechanism is used to estimate the denominator of the j-specific weight:

For subjects who do not switch therapy before the outcome is measured m months later, the denominator of the j-specific weight is

For subjects who have not switched therapy by time j but switch at some point T = j + C(j) before the outcome is measured (C(j) < m), the denominator of the j-specific weight is

Recall that subjects who have already switched therapy by time j do not contribute to our counterfactuals of interest.

The choice of a numerator for the weights will not affect the consistency of the IPTW estimator, as long as the numerator is only a function of treatment history and baseline covariates at time j. One common approach is to use “stabilized” weights (2). As applied to HA-MSMs, this approach involves estimating the following numerator for each j-specific weight:

where Ā(j,l − 1) = (A(j), …,A(l − 1)). The numerator can be estimated by simply fitting a model of the treatment decision at each time point, analogous to the model of the treatment mechanism fit for the denominator, but now including only those covariates contained in the effect modifiers of interest at time j.

Note that the same subject will have a separate weight for each time point j in the study, with denominators corresponding to the probability that the subject received her observed treatment from that time point j until the outcome is measured. Once each subject has been assigned a set of K + 1 − m weights, a weighted least-squares regression analysis is conducted using standard software, with each subject contributing K + 1 − m weighted lines of data. Because each subject contributes multiple lines of data, an approach such as the nonparametric bootstrap or robust standard error estimation using generalized estimating equations is needed to ensure accurate inference.

Modification of nonsuppressive antiretroviral therapy

Our results support a clinical scenario in which the major immunologic effects of exposure to nonsuppressive therapy accrue during the early months following virologic failure. One of the primary detrimental effects of delayed switching is hypothesized to arise because of continuing evolution of the virus in the presence of the drug, resulting in the accumulation of additional mutations. It may prove to be the case that the majority of such evolution has occurred within the first 5 months of nonsuppressive therapy. The changing role of CD4 T-cell count may also be due in part to the fact that the population remaining on nonsuppressive therapy for at least 5 months may consist disproportionately of persons better able to tolerate it.

Finally, we emphasize that the clinical relevance of the results presented here is limited by several factors: 1) the sample size was small; 2) due in part to sample size, we used a broad definition of treatment modification, which included treatment interruptions as well as initiation of new regimens with the goal of viral resuppression; 3) 19 of the persons included in the baseline MSM analysis achieved viral resuppression without treatment modification, suggesting that these persons did not fail therapy because of resistance; and 4) we made the assumption that time until switching had a linear effect on the outcome, which may or may not be an interesting summary measure of the true causal relation. Because of these concerns, the results presented are intended as an illustration of how HA-MSMs can be applied to a real question, rather than as a guide to clinical practice.

History-adjusted MSMs

HA-MSMs generalize MSM methodology. Like MSMs, HA-MSMs control for time-dependent confounding; in addition, HA-MSMs can be used to estimate time-dependent causal effect modification. As we have demonstrated, such an approach can be used to address practical clinical questions.

In the data example presented, loss of virologic failure represented a clear baseline. In many analyses, however, the assignment of a given time point as baseline is essentially arbitrary (for example, when baseline simply corresponds to entry into a study). In such cases, HA-MSMs offer an additional advantage, allowing gains in the precision of coefficient estimates by assuming a common model across time points. As a result of both this gain in precision and the ability of HA-MSMs to estimate effect modification in longitudinal data by time-varying covariates, we anticipate that these models will prove useful in multiple fields of applied research.