Worked Example of a Weighted Cox Survival MSM

From Treatment Histories to a Marginal Cox Model

A weighted Cox survival marginal structural model (MSM) uses treatment-history weights to fit a hazard model for a target population in which measured treatment-confounder feedback has been addressed by weighting. The calculation is a Cox partial-likelihood calculation, but the risk sets are weighted: each person contributes according to their treatment history rather than as one unweighted unit. The resulting treatment coefficient describes a hazard contrast under the specified MSM; it is not itself an absolute survival difference.

For a simple example, suppose the two interventions of interest are sustained treatment $a=0$ and sustained treatment $a=1$. A proportional-hazards MSM is

\lambda^a(t)=\lambda_0(t)\exp(\beta a),\qquad a\in\{0,1\},

where $\lambda^a(t)$ is the counterfactual hazard under regime $a$, $\lambda_0(t)$ is the baseline hazard, and $\exp(\beta)$ is the regime-1 versus regime-0 hazard ratio. For a general longitudinal treatment history, the fitted covariate must encode the treatment contrast specified by the scientific question—for example, a regime indicator or a time-varying treatment term. A coefficient for current treatment should not be interpreted as the effect of a sustained regime unless the model and design support that interpretation.

Definition: Let $N_i(t)$ indicate whether person $i$ has experienced the event by time $t$, and let $Y_i(t)$ indicate whether they are at risk just before $t$. Let $A_i(t)$ be the treatment covariate in the chosen MSM and $W_i(t)$ the applicable treatment-history weight. At event time $t_j$, the weighted Cox risk-set denominator is $$ S_w^{(0)}(\beta,t_j)=\sum_{i=1}^{n}Y_i(t_j)W_i(t_j)\exp\{\beta A_i(t_j)\}. $$ For distinct event times, the weighted partial log likelihood is $$ \ell_w(\beta)=\sum_j W_{(j)}(t_j)\left[\beta A_{(j)}(t_j)-\log S_w^{(0)}(\beta,t_j)\right], $$ where $(j)$ denotes the person with the event at $t_j$. The event contribution uses that person's applicable weight. This tutorial assumes the weights are nonnegative and fixed for the calculation; estimated-weight uncertainty requires additional care.

The target interpretation relies on substantive and statistical conditions, not on the partial likelihood alone. In particular, causal interpretation requires consistency, sequential exchangeability given the measured history, and treatment positivity for the regimes being compared. The weighted pseudo-population must be adequately represented by the chosen MSM, including its proportional-hazards restriction if that is the model fitted. For the Cox calculations below, assume independent subjects, correctly constructed risk sets, no tied event times, and weights that are predictable with respect to the event process. Ties require an explicitly chosen approximation or an appropriate tied-event method.

A Small Weighted Partial-Likelihood Fit

The first example isolates the arithmetic of weighted risk sets. Think of two independent strata, each representing a separate baseline-hazard stratum. In each stratum, just before its one event, the risk set contains treatment $A=1$ with total weight 1 and treatment $A=0$ with total weight 1. The event in the first stratum occurs in an $A=1$ person; the event in the second occurs in an $A=0$ person. Event weights are 1. The common treatment coefficient is estimated, while each stratum has its own baseline hazard.

Worked Example: Estimating a Common Treatment Coefficient

For either event, the weighted denominator is $1+\exp(\beta)$. The two event contributions therefore give

\ell_w(\beta)=\beta-\log\{1+\exp(\beta)\}-\log\{1+\exp(\beta)\} =\beta-2\log\{1+\exp(\beta)\}.

Differentiating,

U_w(\beta)=\frac{d\ell_w}{d\beta} =1-\frac{2\exp(\beta)}{1+\exp(\beta)}.

Set the score to zero. Then $\exp(\hat\beta)/(1+\exp(\hat\beta))=1/2$, so $\exp(\hat\beta)=1$ and $\hat\beta=0$. The fitted hazard ratio is $\exp(\hat\beta)=1$: this small example estimates equal hazards under the two regimes. The information is

\mathcal I(\hat\beta) =\sum_{j=1}^{2}p_j(1-p_j) =2(1/2)(1/2)=1/2, \qquad \operatorname{SE}_{\rm model}(\hat\beta)=\sqrt{2}.

Here $p_j=\exp(\hat\beta)/(1+\exp(\hat\beta))=1/2$ is the fitted treatment share in the risk-set denominator. The estimate is finite because the event pattern includes both a treated and an untreated event, while both treatment levels are represented in each risk set.

This derivation illustrates a general result. Define $S_w^{(1)}(\beta,t_j)=\sum_iY_i(t_j)W_i(t_j)\exp\{\beta A_i(t_j)\}A_i(t_j)$ and $S_w^{(2)}$ by replacing $A_i(t_j)$ with $A_i(t_j)^2$. Differentiating the log denominator gives $S_w^{(1)}/S_w^{(0)}$, the weighted risk-set mean of treatment. Differentiating once more gives its weighted variance.

Result: For the stated distinct-event weighted partial log likelihood, $$ U_w(\beta)=\sum_j W_{(j)}(t_j)\left[A_{(j)}(t_j)-\frac{S_w^{(1)}(\beta,t_j)}{S_w^{(0)}(\beta,t_j)}\right], $$ and $$ -\ell_w''(\beta)=\sum_j W_{(j)}(t_j)\left[ \frac{S_w^{(2)}(\beta,t_j)}{S_w^{(0)}(\beta,t_j)} -\left\{\frac{S_w^{(1)}(\beta,t_j)}{S_w^{(0)}(\beta,t_j)}\right\}^{2} \right]. $$ Thus the negative second derivative is a sum of event-weighted, weighted risk-set variances. It is nonnegative, so the log partial likelihood is concave. It is strictly concave over a parameter range if the contributing risk sets have positive weighted variation in the modeled covariate there.

To verify the result, write $S_w^{(0)}=\sum_i q_i\exp(\beta A_i)$, where $q_i=Y_i(t_j)W_i(t_j)$ at the event time. Its derivative is $S_w^{(1)}$, so the first derivative of $\log S_w^{(0)}$ is $S_w^{(1)}/S_w^{(0)}$. The derivative of that ratio is $S_w^{(2)}/S_w^{(0)}-(S_w^{(1)}/S_w^{(0)})^2$. This is the variance of $A_i$ under probabilities proportional to $q_i\exp(\beta A_i)$, hence is nonnegative. Applying these derivatives to each event term proves both formulas. $\square$

From the Fitted Hazard to Survival

After estimating $\beta$, estimate the baseline cumulative hazard using a weighted Breslow increment. With distinct event times and event weights as in the partial likelihood, the increment at $t_j$ is

\Delta\widehat{\Lambda}_0(t_j) =\frac{W_{(j)}(t_j)}{S_w^{(0)}(\hat\beta,t_j)}.

The fitted survival curve under regime $a$ is then $\widehat S^a(t)=\exp\{-\exp(\hat\beta a)\widehat\Lambda_0(t)\}$, with the cumulative baseline estimate formed by summing the increments through time $t$. This step relies on the proportional-hazards MSM: the regime-specific hazard multiplies the baseline by a constant $\exp(\beta a)$. It is not a direct estimate of a separate hazard shape for each regime.

Worked Example: Calculating Absolute Survival Effects

Suppose a fitted model gives $\hat\beta=-0.4055$, so $\exp(\hat\beta)\approx0.6666$. Suppose the estimated baseline cumulative hazard through a chosen time is $0.12+0.08=0.20$. Under regime 0,

\widehat S^0(t)=\exp(-0.20)=0.8187, \qquad \widehat P(T^0\le t)=1-0.8187=0.1813.

Under regime 1,

\widehat S^1(t)=\exp\{-0.20\exp(-0.4055)\} \approx \exp(-0.1333)=0.8752, \qquad \widehat P(T^1\le t)=1-0.8752=0.1248.

The estimated survival difference is $0.8752-0.8187=0.0565$, about 5.65 percentage points in favor of regime 1 at this time. The hazard ratio is approximately 0.667, but the absolute survival difference depends on the baseline cumulative hazard and the time point; it is not determined by the hazard ratio alone.

The connection between the Breslow estimate and survival follows by integrating the fitted hazard. Under regime $a$, $\widehat\lambda^a(t)=\exp(\hat\beta a)\,d\widehat\Lambda_0(t)/dt$ in the continuous-time notation. Therefore the cumulative hazard through $t$ is $\exp(\hat\beta a)\widehat\Lambda_0(t)$. Since survival is the exponential of minus the cumulative hazard, the displayed survival formula follows. This argument presumes a valid fitted baseline cumulative hazard and does not remove the need for the causal assumptions behind the weighted model. $\square$

Uncertainty, Diagnostics, and Limits

The inverse observed information gives the usual model-based variance when the partial-likelihood model and its independent-observation variance assumptions are appropriate. In a weighted analysis, subjects' contributions can be unequal, and the model-based variance can understate uncertainty. A subject-level sandwich variance replaces the model-based variance with a “bread” term from the derivative of the score and a “meat” term from the outer products of subject-level score contributions. If a person contributes multiple records, those records must be combined into one subject-level contribution before forming the meat. This protects against some forms of variance misspecification; it does not correct confounding, positivity failure, or a wrong MSM.

For illustration, suppose a fitted model reports $\hat\beta=-0.4055$, with model-based information contributions $2(4/7)(3/7)$, $(2/11)(9/11)$, and $(4/7)(3/7)$ across its risk sets. Then

\mathcal I(\hat\beta) =2(4/7)(3/7)+(2/11)(9/11)+(4/7)(3/7) =\frac{36}{49}+\frac{18}{121} =\frac{5238}{5929} \approx0.88345.

The corresponding inverse-information standard error is approximately $1.064$. Using these displayed rounded inputs, a model-based 95% Wald interval is

-0.4055\pm1.96(1.064) =[-2.49094,\ 1.67994] \approx[-2.491,\ 1.680].

This interval is wide, reflecting limited information. It concerns the log hazard ratio; exponentiating its endpoints gives a hazard-ratio interval. A robust interval may differ because its variance estimate uses subject-level score variability rather than only the information. If treatment weights were estimated, treating them as fixed can omit additional uncertainty. A subject-level bootstrap that re-estimates the weights and refits the Cox model is one way to propagate that uncertainty, provided the resampling and fitting procedure is appropriate.

Useful checks include examining the distribution of weights and the effective representation of each treatment regime in each event risk set. Very large weights, nearly empty weighted risk sets, or extreme changes in weighted covariate composition can make estimates unstable and signal weak positivity. The weighted risk-set treatment means and event counts also help reveal whether the coefficient is informed by meaningful within-risk-set contrasts. Check proportional hazards with suitable residual or time-interaction diagnostics; a single coefficient may obscure effects that vary over follow-up.

A common mistake is to read $\exp(\hat\beta)$ as a risk ratio, or to say that a hazard ratio of 0.67 means 33% fewer events by a specified time. It is a ratio of instantaneous hazards among those still event-free at each time, under the model. Absolute risks or survival differences require the baseline hazard and a time horizon, as in the calculation above. Likewise, robust standard errors address a variance calculation, not a failure of treatment positivity or a misspecified causal model.

Compared with an unweighted Cox fit, the weighted fit targets a marginal contrast in the weighted pseudo-population, whereas an unweighted fit generally describes a conditional association given its included covariates. Compared with a pooled logistic survival MSM, the Cox model uses continuous event times and a partial likelihood; pooled logistic models organize follow-up into intervals and model discrete interval event probabilities. These alternatives have different time-scale and approximation considerations. Neither method makes poor weights, sparse events, or an incorrect causal model harmless.

Check Your Understanding

Try each question, then check the brief solution that follows it.

If a weighted risk set has treatment-1 weight total 3 and treatment-0 weight total 2, write its Cox denominator for a binary treatment coefficient $\beta$.
In the two-stratum example, why is the fitted coefficient finite rather than tending to positive or negative infinity?
For $\hat\beta=\log(0.5)$ and baseline cumulative hazard $0.10$, calculate fitted survival under regimes 0 and 1.
What does a subject-level sandwich variance change, and what does it not fix?

Solutions. (1) The denominator is $2+3\exp(\beta)$: the control total contributes $2\exp(0)$, and the treated total contributes $3\exp(\beta)$. (2) Both treatment levels occur in the risk sets, and the observed events include one event in each treatment group, so neither extreme coefficient can keep increasing the partial likelihood; the score is zero at $\beta=0$. (3) Regime 0 survival is $\exp(-0.10)=0.9048$. Regime 1 survival is $\exp\{-0.10\exp(\log(0.5))\}=\exp(-0.05)=0.9512$. (4) The sandwich changes the variance estimate to reflect subject-level score variability, including dependence among a person's repeated records when contributions are correctly clustered. It does not fix confounding, positivity failure, a wrong treatment-weight model, or a misspecified proportional-hazards MSM.

Worked Example of a Weighted Cox Survival MSM

What You'll Learn