---
title: "Longitudinal Treatments"
author: "Noah Greifer"
date: "`r Sys.Date()`"
output:
    html_vignette:
        df_print: kable
        toc: true
vignette: >
  %\VignetteIndexEntry{Longitudinal Treatments}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown_notangle}
bibliography: references.bib
link-citations: true
editor_options:
  chunk_output_type: console
---

```{=html}
<style>
pre {
  overflow-x: auto;
}
pre code {
  word-wrap: normal;
  white-space: pre;
}
</style>
```

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, message = FALSE,
                      fig.width = 7, fig.height = 4, fig.align = "center")
options(width = 200, digits = 4)

me_ok <- rlang::is_installed("marginaleffects (>= 0.29.0)")

fmt <- function(x, digits = 2) {
  formatC(x, digits = digits, format = "f")
}
```

## Introduction

A longitudinal treatment is one that can be received, or not, at each of several time points, so that each unit has a treatment history rather than a single treatment status. Examples include a medication that can be started, stopped, and restarted at each clinic visit, a program a student can enroll in each semester, and a policy a jurisdiction can adopt in any year. The causal questions in these settings concern the joint effect of the treatments received across time: what would happen if units were treated at every time point compared with none, or at the first time point only compared with the last. Answering them requires accounting for *time-varying confounding*, in which covariates that are affected by earlier treatments go on to affect later treatments and the outcome. Standard regression adjustment cannot handle this [@danielMethodsDealingTimedependent2013; @mansourniaHandlingTimeVarying2017], and weighting is the usual remedy: weights are estimated for each unit's treatment history and used to fit a *marginal structural model* (MSM) for the outcome as a function of the treatment history [@robinsMarginalStructuralModels2000; @hernanMarginalStructuralModels2000]. For introductions to MSMs, see @thoemmesPrimerInverseProbability2016, @vanderweeleCausalInferenceLongitudinal2016, @coleConstructingInverseProbability2008, @williamsonMarginalStructuralModels2017, @stallworthyInvestigatingCausalQuestions2026, and @hernanCausalInferenceWhat2020.

In *WeightIt*, weights for longitudinal treatments are estimated with `weightitMSM()`, which fits a model for the treatment at each time point and multiplies the resulting weights together. This vignette explains the estimand and assumptions that are specific to longitudinal treatments and then walks through an analysis: estimating the weights, assessing balance at each time point, and estimating the joint effect of the treatments with an MSM. A second, briefer analysis shows how to account for loss to follow-up by including censoring weights in the same workflow. The two analyses use the simulated `msmdata` dataset included with *WeightIt*. For an introduction to the functions in *WeightIt*, see `vignette("WeightIt")`; for the weighting methods available, see `vignette("weighting-methods")`; for effect estimation in general, see `vignette("estimating-effects")`; and for a fuller treatment of balance assessment with longitudinal treatments, see `vignette("longitudinal-treat", package = "cobalt")`.

## Marginal Structural Models for Longitudinal Treatments

### The Estimand

Consider a study with $K$ time points. At each time point $k$, covariates $L_k$ are measured and then a treatment $A_k$ is received; $L_1$ contains the baseline covariates, and for $k > 1$, $L_k$ contains the covariates measured after treatment $A_{k-1}$ and before treatment $A_k$. The outcome $Y$ is measured after the last treatment. We write $\bar{A}_k = (A_1, \ldots, A_k)$ for the treatment history through time $k$ and $\bar{L}_k = (L_1, \ldots, L_k)$ for the covariate history, with $\bar{A} = \bar{A}_K$ the full treatment history. A *treatment regime* $\bar{a} = (a_1, \ldots, a_K)$ is a particular sequence of treatment values, and $Y(\bar{a})$ is the potential outcome a unit would have had under regime $\bar{a}$. With a binary treatment and three time points, there are eight regimes, including never treated, $(0, 0, 0)$, and always treated, $(1, 1, 1)$.

The estimands of interest are the expected potential outcomes under each regime, $E[Y(\bar{a})]$, and contrasts between them. The contrast between always treated and never treated, $E[Y(1, 1, 1)] - E[Y(0, 0, 0)]$, is a common one, but any pair of regimes can be compared, and the pattern across regimes answers questions about timing and duration, such as whether treatment at the last time point matters more than treatment at the first [@stallworthyInvestigatingCausalQuestions2026]. These are *joint* effects of the entire regime, not the effect of any one treatment holding the others fixed, and they are defined for the population from which the sample was drawn, so the estimand corresponds to the average treatment effect in the population (ATE); the other estimands available for point treatments do not apply.

A marginal structural model is a model for $E[Y(\bar{a})]$ as a function of the regime [@robinsMarginalStructuralModels2000; @robinsMarginalStructuralModels2000a]. It is *marginal* because it concerns the marginal—rather than joint—distribution of the potential outcomes [@breskinExploringSubtletiesInverse2018], and *structural* because it describes potential rather than observed outcomes. For a binary outcome and three time points, a *saturated* MSM has a parameter for every regime,

$$
\begin{aligned}
\text{logit}\, E[Y(a_1, a_2, a_3)] &= \beta_0 + \beta_1 a_1 + \beta_2 a_2 + \beta_3 a_3 \\
&\quad + \beta_4 a_1 a_2 + \beta_5 a_1 a_3 + \beta_6 a_2 a_3 + \beta_7 a_1 a_2 a_3
\end{aligned}
$$ and imposes no assumptions about how the regimes relate to each other. A parsimonious MSM, such as one with only the main effects of the three treatments or one in which the outcome depends only on the number of time points treated, has fewer parameters and is more precise, but it is a modeling assumption that can be wrong. With few time points, the saturated model is preferable; with many time points, a parsimonious model may be a practical necessity.

### Time-Varying Confounding

With a point treatment, confounding is addressed by adjusting for the covariates that affect both the treatment and the outcome. With a longitudinal treatment, the covariates measured between treatments play a double role. A covariate $L_2$ measured after $A_1$ may be affected by $A_1$, affect $A_2$, and affect $Y$: it is a confounder of the effect of $A_2$ and, at the same time, a mediator of the effect of $A_1$. Adjusting for $L_2$ in a regression of $Y$ on the treatment history removes the confounding of $A_2$ but blocks the part of the effect of $A_1$ that runs through $L_2$ and, if $L_2$ and $Y$ share an unmeasured cause, opens a non-causal path between $A_1$ and $Y$ (collider stratification). Not adjusting for $L_2$ leaves the effect of $A_2$ confounded. No choice of covariates for the regression yields the joint effect of the regime [@robinsNewApproachCausal1986; @robinsMarginalStructuralModels2000; @danielMethodsDealingTimedependent2013; @hernanCausalInferenceWhat2020].

Weighting escapes this dilemma because it uses the covariates to construct the weights rather than as terms in the outcome model. In the weighted sample, each treatment is independent of the covariate and treatment history that precedes it, as it would be in a sequentially randomized trial in which treatment at each time point is assigned at random, possibly with probabilities that depend on the past. The MSM is then fit to the weighted sample with only the treatment history (and, optionally, baseline covariates) on the right hand side, so the covariates affected by earlier treatments never enter the outcome model.

### Assumptions

The weighting estimator of an MSM identifies the joint effects under the longitudinal versions of the usual assumptions [@robinsMarginalStructuralModels2000; @hernanCausalInferenceWhat2020].

**Sequential exchangeability.** At each time point, treatment is independent of the potential outcomes given the observed history: $Y(\bar{a}) \perp A_k \mid \bar{A}_{k-1} = \bar{a}_{k-1}, \bar{L}_k$ for every $k$ and every regime $\bar{a}$. This is the assumption of no unmeasured confounding applied one time point at a time, and it requires that the covariates measured before each treatment include all the confounders of that treatment and the outcome, including the earlier treatments themselves. It cannot be verified from the data; its plausibility depends on what was measured and when.

**Positivity.** Every regime must have a positive probability of being received at every history that occurs in the population: $P(A_k = a_k \mid \bar{A}_{k-1} = \bar{a}_{k-1}, \bar{L}_k = \bar{l}_k) > 0$ for all $a_k$ and all histories $(\bar{a}_{k-1}, \bar{l}_k)$ with positive probability. Positivity is harder to satisfy with longitudinal treatments than with point treatments because it must hold at every time point, and because some treatment histories may be rare or absent by design (e.g., when a treatment, once started, is never stopped). Near-violations show up as very large weights [@westreichInvitedCommentaryPositivity2010; @plattPositivityAssumptionMarginal2012; @coleConstructingInverseProbability2008].

**Consistency and no interference.** A unit's observed outcome is its potential outcome under the regime it actually received, and a unit's outcome does not depend on the treatments other units receive.

**Correct specification.** The models for treatment at each time point must be correctly specified, or at least must produce weights that balance the covariates at each time point, and the MSM must be correctly specified as a function of the regime. Misspecification of the treatment models biases the weights and the estimate [@lefebvreImpactMisspecificationTreatment2008], which is why balance is assessed at every time point.

### Inverse Probability Weights

The weights for an MSM are the product over time points of the inverse of the probability of the treatment actually received, given the history to that point [@robinsMarginalStructuralModels2000]: $$w_i = \prod_{k=1}^{K} \frac{1}{P(A_k = A_{ik} \mid \bar{A}_{k-1} = \bar{A}_{i,k-1}, \bar{L}_k = \bar{L}_{ik})}.$$ Each factor is the weight that would be estimated for a point treatment at time $k$ with the history through time $k$ as the covariates, and this is how `weightitMSM()` computes them: it fits a treatment model at each time point using the formula supplied for that time point and multiplies the weights across time points. The factors can be estimated with any of the methods described in `vignette("weighting-methods")` that support longitudinal treatments, which are the methods that estimate a propensity score (`"glm"`, `"gbm"`, `"super"`, `"bart"`, `"cbps"`, and `"ipt"`); the optimization-based methods, which do not, cannot be used this way. The covariate balancing propensity score can alternatively estimate all the time points' models at once so that the product of the weights balances the covariates at every time point [@huffmanCovariateBalancingInverse2018].

The product of many inverse probabilities can be very variable, and the weights for a longitudinal treatment are often much more extreme than those for a point treatment. *Stabilized* weights replace the numerator of each factor with the probability of the observed treatment given the treatment history alone: $$sw_i = \prod_{k=1}^{K} \frac{P(A_k = A_{ik} \mid \bar{A}_{k-1} = \bar{A}_{i,k-1})}{P(A_k = A_{ik} \mid \bar{A}_{k-1} = \bar{A}_{i,k-1}, \bar{L}_k = \bar{L}_{ik})}.$$ Stabilized weights have a mean of 1 (when correctly specified), are far less variable than unstabilized weights, and yield more precise estimates [@robinsMarginalStructuralModels2000; @coleConstructingInverseProbability2008]. They are requested with `stabilize = TRUE`, which fits a saturated model of each treatment on the prior treatments for the numerator. The numerator does not change the target population, but it changes what the weights balance: in the stabilized weighted sample, each treatment is independent of the covariate history *given* the treatment history, and the association between treatments at different time points is preserved. This has two consequences. The MSM must include the treatment history flexibly enough to absorb that association, which the saturated MSM does, and balance on the prior treatments is not expected when assessing balance with stabilized weights, as we will see below. Baseline covariates can also be included in the numerator, which further reduces the variability of the weights at the cost of changing the interpretation of the MSM unless the same covariates are included in it; see @coleConstructingInverseProbability2008 for details and the `num.formula` argument of `weightitMSM()` for how to request this.

Balance is assessed one time point at a time: at time $k$, the weighted sample should show no association between $A_k$ and any variable in $\bar{A}_{k-1}$ or $\bar{L}_k$. The *cobalt* package computes the balance statistics for each time point and summarizes them across time points, and the effective sample size at each time point indicates the precision retained. These are the same criteria used to choose among weighting specifications for a point treatment, applied at each time point; see `vignette("longitudinal-treat", package = "cobalt")` for the tools and @jacksonDiagnosticsConfoundingTimevarying2016 for an alternative set of diagnostics.

## Estimating Weights for a Longitudinal Treatment

Below, we demonstrate estimating weights and the treatment effect for a longitudinal treatment using *WeightIt*.

### The Data

We use `msmdata`, a simulated dataset of 7500 units from a hypothetical study with three treatment periods. At each period, two covariates are measured and then a binary treatment is received, and a binary adverse outcome is measured after the third period. The treatments and later covariates were generated to depend on the earlier covariates and treatments, so the data exhibit time-varying confounding. The dataset is in "wide" format, with one row per unit and a separate column for each variable at each time point, which is the format `weightitMSM()` requires; `reshape()` can be used to convert data from long to wide format.

```{r}
library("WeightIt")
library("cobalt")

data("msmdata")

head(msmdata)
```

The baseline covariates are `X1_0` (a count) and `X2_0` (binary); `A_1` is the first treatment; `X1_1` and `X2_1` are measured after `A_1` and before `A_2`; `X1_2` and `X2_2` are measured after `A_2` and before `A_3`; and `Y_B` is the outcome. In the notation above, $L_1 = (X1_0, X2_0)$, $L_2 = (X1_1, X2_1)$, and $L_3 = (X1_2, X2_2)$.

### Initial Imbalance

The treatment model at each time point includes everything measured before that treatment: at the first time point, the baseline covariates; at the second, the baseline covariates, the first treatment, and the covariates measured after it; and at the third, all of these plus the second treatment and the covariates measured after it. We specify these as a list of formulas in temporal order and supply the list to `bal.tab()` to examine balance before weighting. Setting `which.time = .all` displays a balance table for each time point; the default, `which.time = .none`, displays only the summary across time points.

```{r}
bal.tab(list(A_1 ~ X1_0 + X2_0,
             A_2 ~ X1_1 + X2_1 + A_1 + X1_0 + X2_0,
             A_3 ~ X1_2 + X2_2 + A_2 + X1_1 + X2_1 + A_1 + X1_0 + X2_0),
        data = msmdata, stats = c("m", "ks"),
        which.time = .all)
```

Every covariate is imbalanced at every time point, with standardized mean differences well above .1, and the treatments at later time points are associated with the treatments before them. A model for `Y_B` that adjusted for all of these would block the effects of the earlier treatments that run through the later covariates, so we instead use them to estimate weights.

### Estimating the Weights

The same list of formulas is supplied to `weightitMSM()`. With the default `method = "glm"`, a logistic regression model is fit for each treatment on the variables in its formula, and the weights are the product of the inverse predicted probabilities. We start with unstabilized weights.

```{r}
W_un <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
                         A_2 ~ X1_1 + X2_1 + A_1 + X1_0 + X2_0,
                         A_3 ~ X1_2 + X2_2 + A_2 + X1_1 + X2_1 +
                           A_1 + X1_0 + X2_0),
                    data = msmdata, method = "glm")

W_un
```

Printing the object displays the treatment and covariates at each time point. `summary()` describes the distribution of the weights, once for each time point, since the same weights are summarized within the treatment groups defined at each time point.

```{r}
summary(W_un)
```

The weights are extremely variable, with a maximum of about `r round(max(W_un$weights))` and effective sample sizes around a quarter of the group sizes at the first two time points. This is typical of unstabilized weights for a longitudinal treatment, because the product of three inverse probabilities, each of which can be large, is larger still. We refit the weights with `stabilize = TRUE`, which multiplies each factor by the probability of the observed treatment given the prior treatments.

```{r}
W <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
                      A_2 ~ X1_1 + X2_1 + A_1 + X1_0 + X2_0,
                      A_3 ~ X1_2 + X2_2 + A_2 + X1_1 + X2_1 +
                        A_1 + X1_0 + X2_0),
                 data = msmdata, method = "glm",
                 stabilize = TRUE)

W
```

The printout now lists the stabilization factors, which are the prior treatments and their interactions at each time point.

```{r}
summary(W)
```

The largest weight has dropped to about `r round(max(W$weights))`, and the effective sample sizes have risen at every time point. The weights have a mean close to 1, as stabilized weights should.

### Assessing Balance

Supplying the `weightitMSM` object to `bal.tab()` assesses balance in the weighted sample at each time point.

```{r}
bal.tab(W, stats = c("m", "ks"), which.time = .all)
```

The covariates are balanced at every time point, with all standardized mean differences and Kolmogorov-Smirnov statistics near or below .1. The prior treatments are a different matter: `A_1` remains associated with `A_2` at the second time point, and `A_2` with `A_3` at the third, with standardized mean differences about as large as before weighting. This is expected with stabilized weights, whose numerator preserves the association among the treatments, and it is not a deficiency of the weights as long as the MSM includes the treatment history, which the saturated MSM we fit below does. Because the covariates are themselves associated with the prior treatments, small imbalances in the covariates can also appear in these marginal comparisons with stabilized weights; what the weights guarantee is balance within levels of the treatment history. With unstabilized weights, the prior treatments would be balanced as well, which can be verified by running `bal.tab(W_un, which.time = .all)`.

The summary across time points reports the largest imbalance for each variable across the time points at which it appears, which is a convenient way to check the whole specification at once; it is displayed by default when `which.time` is not set.

```{r}
bal.tab(W, stats = c("m", "ks"))
```

`love.plot()` displays the same information graphically, with one panel per time point.

```{r}
love.plot(W, stats = "m", binary = "std", abs = TRUE,
          thresholds = .1, which.time = .all)
```

The one covariate that remains slightly above the .1 threshold is `X1_2` at the third time point. Had this been larger, we would try another specification, as for a point treatment: adding squared terms or interactions to the treatment models, changing the method, or both. One option specific to longitudinal treatments is the covariate balancing propensity score, which estimates all the treatment models at once so that the product of the weights exactly balances the covariate means at every time point [@huffmanCovariateBalancingInverse2018].

```{r}
W_cbps <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
                           A_2 ~ X1_1 + X2_1 + A_1 + X1_0 + X2_0,
                           A_3 ~ X1_2 + X2_2 + A_2 + X1_1 + X2_1 +
                             A_1 + X1_0 + X2_0),
                      data = msmdata, method = "cbps")

bal.tab(W_cbps, stats = c("m", "ks"))
```

Mean balance is now exact at every time point, including on the prior treatments, at the cost of a lower effective sample size than the stabilized logistic regression weights achieve. Because M-estimation is not available for this version of CBPS, standard errors after weighting would have to be bootstrapped. The two specifications represent the usual trade-off between balance and precision, and either could be defended here; we proceed with the stabilized logistic regression weights, whose balance is adequate and which support standard errors that account for the estimation of the weights.

### Estimating the Treatment Effect

The MSM is fit with `glm_weightit()`, which incorporates the weights and, for `method = "glm"`, computes standard errors that account for their estimation using M-estimation, including the estimation of the stabilization factors. The outcome model includes the three treatments and all their interactions, so that it is saturated in the regime, along with the baseline covariates interacted with the treatments. Baseline covariates are measured before any treatment and are the only covariates that may be included; including any of the later covariates would reintroduce the problem that weighting was used to avoid. See `vignette("estimating-effects")` for the general procedure and for bootstrap standard errors, which are requested with the `vcov` argument.

```{r}
fit <- glm_weightit(Y_B ~ A_1 * A_2 * A_3 * (X1_0 + X2_0),
                    data = msmdata, weightit = W,
                    family = binomial)
```

The coefficients of this model are not themselves of interest. Instead, we use g-computation through the *marginaleffects* package to compute the expected potential outcome under each regime, which is the average over the sample of the predicted probability of the outcome with the treatments set to the values of that regime. Supplying the three treatments to `variables` in `avg_predictions()` produces one estimate per regime.

```{r, eval = me_ok}
library("marginaleffects")

p <- avg_predictions(fit, variables = c("A_1", "A_2", "A_3"))

p
```

Each row is a regime. The first row, with all three treatments set to 0, estimates that `r if (me_ok) fmt(100 * p$estimate[1], 1)`% of units would experience the adverse event if no one were treated at any time point, and the last row, with all three set to 1, estimates that `r if (me_ok) fmt(100 * p$estimate[8], 1)`% would if everyone were treated at every time point. The rows in between describe the effects of partial regimes; for example, treatment at the third time point alone (the second row) reduces the risk by more than treatment at the first time point alone (the fifth row).

To compare the regimes, we supply the predictions to `hypotheses()`. Setting the hypothesis to `~reference` contrasts each regime with the first, the never-treated regime.

```{r, eval = me_ok}
hypotheses(p, ~reference)
```

A specific contrast is requested by naming the rows. The joint effect of always being treated relative to never being treated is the difference between the eighth and first rows, and the corresponding log risk ratio is the log of their quotient.

```{r, eval = me_ok}
hypotheses(p, "b8 - b1 = 0")

hypotheses(p, "log(b8 / b1) = 0")
```

```{r, include = FALSE, eval = me_ok}
h <- hypotheses(p, "b8 - b1 = 0")
rr <- hypotheses(p, "log(b8 / b1) = 0")
```

Being treated at all three time points rather than none reduces the risk of the adverse event by `r if (me_ok) fmt(-100 * h$estimate, 1)` percentage points, with a 95% confidence interval from `r if (me_ok) fmt(-100 * h$conf.high, 1)` to `r if (me_ok) fmt(-100 * h$conf.low, 1)` points, which corresponds to a risk ratio of `r if (me_ok) fmt(exp(rr$estimate), 2)`. The standard errors here account for the estimation of the weights and the stabilization factors; because they were computed with M-estimation, no bootstrapping was needed, though bootstrapping remains an option by setting `vcov = "FWB"` or `vcov = "BS"` in `glm_weightit()`.

## Accounting for Loss to Follow-Up

In most longitudinal studies, some units stop being observed before the end of the study, so their later treatments and their outcome are unknown. If the units who drop out differ from those who remain in ways related to the outcome, an analysis restricted to the units who remain estimates the effect in a different population from the one the sample represents, and if dropout depends on a common cause of the treatment and the outcome, restricting to the uncensored units induces a non-causal association between them. Both are forms of selection bias [@hernanStructuralApproachSelection2004; @luClearerDefinitionSelection2022; @hernanCausalInferenceWhat2020]. The remedy is to weight the units who remain so that they resemble the full sample that was at risk of dropping out, using inverse probability of censoring weights (IPCW) [@robinsAnalysisSemiparametricRegression1995; @robinsCorrectingNoncomplianceDependent2000; @hernanMarginalStructuralModels2000; @coleConstructingInverseProbability2008; @seamanReviewInverseProbability2013].

### The Estimand with Censoring

With censoring, the potential outcome is $Y(\bar{a}, \bar{c} = 0)$, the outcome under regime $\bar{a}$ had the unit also remained under observation throughout, and the estimand is the contrast of $E[Y(\bar{a}, \bar{c} = 0)]$ across regimes. Identification requires an exchangeability assumption for censoring parallel to the one for treatment: at each time point, remaining under observation must be independent of the potential outcomes given the observed treatment and covariate history. It also requires positivity for censoring, meaning that every history must have a positive probability of remaining under observation. The censoring weights are the product over time points of the inverse of the probability of remaining under observation given the history, $$w^C_i = \prod_{k} \frac{1}{P(C_k = 0 \mid \bar{A}_k = \bar{A}_{ik}, \bar{L}_k = \bar{L}_{ik}, \bar{C}_{k-1} = 0)},$$ where $C_k = 1$ indicates that the unit was censored at time $k$. The final weight for each unit is the product of its treatment weights and its censoring weights, and a censored unit receives a weight of 0. Censoring weights can be stabilized in the same way as treatment weights, with the probability of remaining under observation given the treatment history in the numerator.

In `weightitMSM()`, a censoring model is a formula whose left side is wrapped in `.cens()`, placed in the list of formulas at the time point at which the censoring occurs. The censoring indicator must be 1 for units censored at that time point and 0 for units still under observation, following the survival convention. Each model in the list is fit only among the units still under observation when it is reached, so the treatments and covariates that are missing for censored units cause no problems. See `?.cens` and `?weightitMSM` for details.

### Simulating Loss to Follow-Up

`msmdata` has no censoring, so we create some. Below, we simulate dropout after the covariates of the third period are measured but before the third treatment is received, with dropout more likely for units with higher values of `X1_2` and for units treated at the second period. Units who drop out have no third treatment and no outcome, so we set those variables to missing for them, as they would be in real data.

```{r}
set.seed(7)

msmdata$C_2 <- rbinom(nrow(msmdata), 1,
                      prob = plogis(-4 + .35 * msmdata$X1_2 + .8 * msmdata$A_2))

# A_3 and Y_B are unobserved for units lost to follow-up
is.na(msmdata[msmdata$C_2 == 1, c("A_3", "Y_B")]) <- TRUE

table(msmdata$C_2)
```

About `r round(100 * mean(msmdata$C_2))`% of the sample is lost to follow-up.

### Estimating the Weights

The censoring model goes in the list of formulas after the second treatment model and before the third, and it includes the same history as the treatment model that follows it, since anything that affects both dropout and the outcome must be included. We request stabilized weights as before.

```{r}
Wc <- weightitMSM(list(A_1 ~ X1_0 + X2_0,
                       A_2 ~ X1_1 + X2_1 + A_1 + X1_0 + X2_0,
                       .cens(C_2) ~ X1_2 + X2_2 + A_2 + X1_1 + X2_1 +
                         A_1 + X1_0 + X2_0,
                       A_3 ~ X1_2 + X2_2 + A_2 + X1_1 + X2_1 +
                         A_1 + X1_0 + X2_0),
                  data = msmdata, method = "glm",
                  stabilize = TRUE)

Wc
```

The printout identifies the third entry as a censoring model and reports how many units were censored. The weights are the product of all four models' weights, and the censored units have weights of exactly 0, which `summary()` reports in its `# Zeros` column; we display the summary for the censoring model and the final treatment model only.

```{r}
summary(Wc, which.time = 3:4)
```

The censoring model's entry summarizes the weights among the units still under observation, and its effective sample size is measured against that sample. The final treatment model was fit among the `r sum(Wc$at.risk[, "A_3"])` units who remained, which the `at.risk` component of the output records.

### Assessing Balance

Balance is assessed as before, with one table per model. For the censoring model, the comparison is between the weighted units who remained under observation and the full sample at risk, rather than between two treatment groups; the weights are designed to make the first resemble the second.

```{r}
bal.tab(Wc, stats = c("m", "ks"), which.time = 3:4)
```

In the censoring table, the standardized mean differences between the weighted uncensored sample and the full sample are small for the covariates, with `X1_2`, the strongest predictor of dropout, the largest, and the association with `A_2` is preserved by the stabilization as with the treatments. At the third treatment, assessed among the units who remained, the covariates are balanced. As with any weighting specification, the balance tables and the effective sample sizes should guide the choice of the censoring model's specification and the weighting method.

### Estimating the Treatment Effect

The outcome model is fit exactly as before. Censored units have a weight of 0 and a missing outcome, and `glm_weightit()` tolerates missing values in the outcome and treatments for units whose weight is 0, so they are dropped from the fit without any change to the call.

```{r}
fit_c <- glm_weightit(Y_B ~ A_1 * A_2 * A_3 * (X1_0 + X2_0),
                      data = msmdata, weightit = Wc,
                      family = binomial)
```

```{r, eval = me_ok}
p_c <- avg_predictions(fit_c, variables = c("A_1", "A_2", "A_3"))

p_c

hypotheses(p_c, "b8 - b1 = 0")
```

```{r, include = FALSE, eval = me_ok}
h_c <- hypotheses(p_c, "b8 - b1 = 0")
```

The estimated effect of always being treated relative to never being treated is a reduction in risk of `r if (me_ok) fmt(-100 * h_c$estimate, 1)` percentage points. Because we introduced the censoring ourselves, we can compare this with the estimate from the full data above, a reduction of `r if (me_ok) fmt(-100 * h$estimate, 1)` points. An analysis that simply dropped the censored units would not recover this estimate in general: the units who dropped out had higher values of `X1_2`, which raises the risk of the outcome, so the units who remained are a lower-risk population in which the treatment contrast differs, and because dropout also depended on `A_2`, restricting to the units who remained induces an association between `A_2` and `X1_2` that treatment weights estimated among those units would not remove. The censoring weights restore the composition of the original sample and break that association. How large the bias from ignoring censoring would be in real data is unknown, which is why censoring should be modeled whenever it is plausibly related to the covariates or treatments.

Note that because the censoring weights multiply the treatment weights, the final weights can be more variable than either set alone, and balance at the time points after the censoring should be examined with particular care. When dropout depends strongly on the history, the censoring weights themselves become extreme, and weighting can perform poorly even when the censoring model is correct [@howeLimitationInverseProbabilityofcensoring2011].

## Other Considerations

The treatments at different time points need not be of the same type: `weightitMSM()` accepts binary, multi-category, and continuous treatments in any combination, with the weights at each time point computed as they would be for a point treatment of that type. The weighting method applies to every time point; to use different methods or options at different time points, the weights can be estimated separately with `weightit()` and multiplied by hand.

As the number of time points grows, the product of the weights becomes more variable, the saturated stabilization model and the saturated MSM gain parameters quickly, and positivity becomes harder to satisfy. Trimming the weights with `trim()`, using a parsimonious stabilization model through `num.formula`, and fitting a parsimonious MSM are the usual responses, each with the costs described above. Whatever the number of time points, the data must have one row per unit, so a dataset in long format must be reshaped to wide format before calling `weightitMSM()`.

Sampling weights can be supplied through `s.weights`, and missing covariate values among the units still under observation are handled as described on the help page for each method. For questions about the dosage and timing of exposures in developmental research, @stallworthyInvestigatingCausalQuestions2026 describe the *devMSMs* package, which provides a structured workflow for specifying, fitting, and interpreting MSMs for such questions. When reporting the analysis, the treatment and censoring models at each time point, the stabilization, the balance achieved at each time point, the effective sample sizes, and the form of the MSM should all be described; see `vignette("estimating-effects")` for reporting the effect estimates.

## References