I am going to cut right to the chase, since this is the third of three posts related to confounding and weighting, and it’s kind of a long one. (If you want to catch up, the first two are here and here.) My aim with these three posts is to provide a basic explanation of the marginal structural model (MSM) and how we should interpret the estimates. This is obviously a very rich topic with a vast literature, so if you remain interested in the topic, I recommend checking out this (as of yet unpublished) text book by Hernán &amp; Robins for starters.
The DAG below is a simple version of how things can get complicated very fast if we have sequential treatments or exposures that both affect and are affected by intermediate factors or conditions.
\(A_0\) and \(A_1\) represent two treatment points and \(L_0\) and \(L_1\) represent measurements taken before and after treatments, respectively. Both treatments and at least \(L_1\) affect outcome \(Y\). (I am assuming that the \(A\)’s and \(L\)’s are binary and that \(Y\) is continuous. \(\epsilon\) is \(N(0, \sigma_\epsilon^2)\).)
An example of this might be a situation where we are interested in the effect of a drug treatment on mental health well-being for patients with prehypertension or hypertension. A physician’s decision to administer the drug at each visit is influenced by the patient’s level of hypertension. In turn, the treatment \(A_0\) potentially reduces the probability of hypertension - \(P(L_1=1)\). And finally, \(L_1\) influences the next treatment decision and ultimately the mental health outcome \(Y\).
The complicating factor is that the hypertension level following the first treatment (\(L_1\)) is both a mediator the effect of treatment \(A_0\) and confounder of the treatment effect \(A_1\) on \(Y\). To get an unbiased estimate the effect of the combined treatment regime (\(A_0\) and \(A_1\)) we need to both control for \(L_1\) and not control for \(L_1\). This is where MSMs and inverse probability weighting (IPW) come into play.
The MSM is marginal in the sense that we’ve been talking about in this series - the estimate will be a population-wide estimate that reflects the mixture of the covariates that influence the treatments and outcomes (in this case, the \(L\)’s). It is structural in the sense that we are modeling potential outcomes. Nothing has changed from the last post except for the fact that we are now defining the exposures as a sequence of different treatments (here \(A_0\) and \(A_1\), but could easily extend to \(n\) treatments - up to \(A_n\).)
Imagine an experiment …
To understand the MSM, it is actually helpful to think about how a single individual fits into the picture. The tree diagram below literally shows that. The MSM posits a weird experiment where measurements (of \(L\)) are collected and treatments (\(A\)) are assigned repeatedly until a final outcome is measured. In this experiment, the patient is not just assigned to one treatment arm, but to both! Impossible of course, but that is the world of potential outcomes.
At the start of the experiment, a measurement of \(L_0\) is collected. This sends the patient down the one of the branches of the tree. Since the patient is assigned to both \(A_0=0\) and \(A_0=1\), she actually heads down two different branches simultaneously. Following the completion of the first treatment period \(A_0\), the second measurement (\(L_1\)) is collected. But, two measurements are taken for the patient - one for each branch. The results need not be the same. In fact, if the treatment has an individual-level effect on \(L_1\), then the results will be different for this patient. In the example below, this is indeed the case. Following each of those measurements (in parallel universes), the patient is sent down the next treatment branches (\(A_1\)). At this point, the patient finds herself in four branches. At the end of each, the measurement of \(Y\) is taken, and we have four potential outcomes for individual {i}: \(Y^i_{00}\), \(Y^i_{10}\), \(Y^i_{01}\), and \(Y^i_{11}\).
A different patient will head down different branches based on his own values of \(L_0\) and \(L_1\), and will thus end up with different potential outcomes. (Note: the values represented in the figure are intended to be average values for that particular path.)
How do we define the causal effect?
With four potential outcomes rather than two, it is less obvious how to define the causal effect. We could, for example, consider three separate causal effects by comparing each of the treatment “regimes” that include at least one exposure to the intervention to the single regime that leaves the patient entirely unexposed. That is, we could be interested in (at the individual \(i\) level) \(E^i_1 = Y^i_{10}-Y^i_{00}\), \(E^i_2 = Y^i_{01}-Y^i_{00}\), and \(E^i_3 = Y^i_{11}-Y^i_{00}\). This is just one possibility; the effects of interest are driven entirely by the research question.
When we have three or four or more intervention periods, the potential outcomes can start to pile up rapidly (we will have \(2^n\) potential outcomes for a sequence of \(n\) treatments.) So, the researcher might want to be judicious in deciding which contrasts to be made. Maybe something like \(Y_{1111} - Y_{0000}\), \(Y_{0111} - Y_{0000}\), \(Y_{0011} - Y_{0000}\), and \(Y_{0001} - Y_{0000}\) for a four-period intervention. This would allow us to consider the effect of starting (and never stopping) the intervention in each period compared to never starting the intervention at all. By doing this, though, we would be using only 5 out of the 16 potential outcomes. If the remaining 11 paths are not so rare, we might be ignoring a lot of data.
The marginal effect
The tree below represents an aggregate set of branches for a sample of 5000 individuals. The sample is initially characterized only by the distribution of \(L_0\). Each individual will go down her own set of four paths, which depend on the starting value of \(L_0\) and how each value of \(L_1\) responds in the context of each treatment arm.
Each individual \(i\) (at least in theory) has four potential outcomes: \(Y^i_{00}\), \(Y^i_{10}\), \(Y^i_{01}\), and \(Y^i_{11}\). Averaging across the sample provides a marginal estimate of each of these potential outcomes. For example, \(E(Y_{00})=\sum_i{Y^i_{00}}/5000\). This can be calculated from the tree as \[(1742*53 + 1908*61 + 392*61 + 958*69)/5000 = 59.7\] Similarly, \(E(Y_{11}) = 40.1\) The sample average causal effects are estimated using the respective averages of the potential outcomes. For example, \(E_3\) at the sample level would be defined as \(E(Y_{11}) - E(Y_{00}) = 40.1 - 59.7 = -19.6\).
Back in the real world
In reality, there are no parallel universes. Maybe we could come up with an actual randomized experiment to mimic this, but it may be difficult. More likely, we’ll have observed data that looks like this:
Each individual heads down his or her own path, receiving a single treatment at each time point. Since this is not a randomized trial, the probability of treatment is different across different levels of \(L_0\) and \(L_1\) and that \(L_0\) and \(L_1\) are associated with the outcome (i.e. there is confounding).
Estimating the marginal effects
In the previous posts in this series, I provided some insight as to how we might justify using observed data only to estimate these sample-wide average potential outcomes. The most important assumption is that when we have measured all confounders, we may be able to say, for example, \(E(Y_{01}) = E(Y | A_0 = 0 \ \&amp; \ A_1 = 1 )\). The potential outcome for everyone in the sample is equal to the observed outcome for the subgroup who actually followed the particular path that represents that potential outcome. We will make the same assumption here.
At the start of this post, I argued that given the complex nature of the data generating process (in particular given that \(L_1\) is both a mediator and confounder), it is challenging to get unbiased estimates of the intervention effects. One way to do this with marginal structural models (another way is using g-computation, but I won’t talk about that here). Inverse probability weighting converts the observed tree graph from the real world to the marginal tree graph so that we can estimate sample-wide average (marginal) potential outcomes as an estimate for some population causal effects.
In this case, the inverse probability weight is calculated as \[IPW = \frac{1}{P(A_0=a_0 | L_0=l_0) \times P(A_1=a_1 | L_0=l_0, A_0=a_0, L_1=l_1)}\] In practice, we estimate both probabilities using logistic regression or some other modeling technique. But here, we can read the probabilities off the tree graph. For example, if we are interested in the weight associated with individuals observed with \(L_0=1, A_0=0, L_1=0, \textbf{and } A_1=1\), the probabilities are \[P(A_0 = 0 | L_0=1) = \frac{676}{1350}=0.50\] and \[P(A_1=1 | L_0=1, A_0=0, L_1=0) = \frac{59}{196} = 0.30\]
So, the inverse probability weight for these individuals is \[IPW = \frac{1}{0.50 \times 0.30} = 6.67\] For the 59 individuals that followed this pathway, the weighted number is \(59 \times 6.67 = 393\). In the marginal world of parallel universes, there were 394 individuals.
Simulating data from an MSM
Before I jump into the simulation, I do want to reference a paper by Havercroft and Didelez that describes in great detail how to generate data from a MSM with time-dependent confounding. It turns out that the data can’t be generated exactly using the intial DAG (presented above), but rather needs to come from something like this:
where \(U\) is an unmeasured, maybe latent, covariate. The observed data (that ignores \(U\)) will indeed have a DAG that looks like the one that we started with.
When doing simulations with potential outcomes, we can generate all the potential outcomes for each individual using a parallel universe approach. The observed data (treatment choices and observed outcomes) are generated separately. The advantage of this is that we can confirm the true causal effects because we have actually generated potential outcomes. The disadvantage is that the code is considerably more complicated and the quantity of data generated grows. The situation is not so bad with just two treatment periods, but the size of the data increases exponentially with the number of treatments: as I mentioned earlier, there will be \(2^n\) potential outcomes for each individual.
Alternatively, we can just generate the observed data directly. Since we know the true causal parameters we actually “know” the causal effects and can compare our estimates.
I will go through the convoluted approach because I think it clarifies (at least a bit) what is going on. As an addendum, I will show how all of this could be done in a few lines of code if we take the second approach …
library(broom)
library(simstudy)
# define U, e and L0
defA0