Time-Varying Treatment Effects in Cox MSMs

Why Let the Treatment Effect Change Over Time?

A Cox marginal structural model with one treatment coefficient assumes that the treatment hazard ratio is constant throughout follow-up. That restriction can be scientifically implausible: a treatment may have an immediate effect that fades, a delayed effect that grows, or effects in opposite directions at different times. A treatment-by-time interaction relaxes the constant-effect restriction while retaining a Cox model for the marginal hazard under each treatment regime.

The treatment variable here is a binary indicator $A$ for one of two regimes, with $A=1$ denoting the treatment regime and $A=0$ the reference regime. Follow-up time is $t$. The model describes the hazard under each regime, not an effect conditional on adjusting for measured confounders in the outcome model. As in other weighted survival MSMs, its causal interpretation depends on the relevant consistency, exchangeability, and positivity conditions, appropriate treatment and censoring weights, and an adequate marginal structural model. Weighting does not make a misspecified outcome model correct.

Definition: A treatment-by-time interaction in a Cox MSM allows the log hazard ratio for treatment to depend on follow-up time. For a specified time function $g(t)$, a linear interaction model is $$ \lambda^a(t)=\lambda_0(t)\exp\{a[\beta_0+\beta_1 g(t)]\},\qquad a\in\{0,1\}. $$ Here $\lambda^a(t)$ is the marginal hazard at time $t$ under regime $a$, $\lambda_0(t)$ is the reference-regime baseline hazard, and $\beta_0,\beta_1$ are coefficients. The function $g(t)$ is specified before fitting; for example, it may be elapsed time or time centered at a chosen value.

In a weighted Cox fit, the treatment-by-time covariate is $A g(t)$. The weighting and risk-set construction are those of the weighted Cox MSM; the interaction changes the model’s linear predictor, not the causal weighting target. For a fixed treatment regime, its value at an event time is determined by that regime and the event time. The model’s central quantity is the coefficient function $\beta(t)=\beta_0+\beta_1g(t)$, whose exponentiation gives the time-specific hazard ratio.

Deriving the Time-Specific Hazard Ratio

The ratio follows directly from the two regime-specific hazards. Under the reference regime, $a=0$, the exponential term equals one. Under treatment, $a=1$, it equals $\exp\{\beta_0+\beta_1g(t)\}$. Thus the baseline hazard cancels:

\operatorname{HR}(t) =\frac{\lambda^1(t)}{\lambda^0(t)} =\frac{\lambda_0(t)\exp\{\beta_0+\beta_1g(t)\}}{\lambda_0(t)} =\exp\{\beta_0+\beta_1g(t)\}.

Result: In the model above, the treatment hazard ratio at time $t$ is $\exp\{\beta_0+\beta_1g(t)\}$. The coefficient $\beta_1$ is a change in the log hazard ratio per unit increase in $g(t)$, not itself a hazard ratio. If $g(t)=0$, the treatment hazard ratio is $\exp(\beta_0)$.

This result is a comparison of hazards at time $t$ among individuals still at risk at that time under the respective regimes. It is not a risk ratio or a ratio of survival probabilities. In particular, a hazard ratio above one indicates a higher instantaneous event rate at that time among those still event-free; it does not by itself quantify the difference in cumulative event probability.

Worked Example: A Linear Increase in the Log Hazard Ratio

Suppose time is measured in years and the fitted MSM is $\lambda^a(t)=\lambda_0(t)\exp\{a[-0.40+0.15t]\}$. Calculate the treatment hazard ratios at years 1, 3, and 5.

Here $g(t)=t$, $\hat{\beta}_0=-0.40$, and $\hat{\beta}_1=0.15$. The estimated log hazard ratio at year $t$ is $-0.40+0.15t$.

At year 1, this is $-0.40+0.15(1)=-0.25$, so $\widehat{\operatorname{HR}}(1)=\exp(-0.25)\approx0.779$. At year 3, it is $-0.40+0.45=0.05$, so $\widehat{\operatorname{HR}}(3)=\exp(0.05)\approx1.051$. At year 5, it is $-0.40+0.75=0.35$, so $\widehat{\operatorname{HR}}(5)=\exp(0.35)\approx1.419$.

The estimated association changes from a lower treatment hazard at year 1 to a slightly higher one at year 3 and a higher one at year 5. The crossing time solves $-0.40+0.15t=0$, giving $t=0.40/0.15\approx2.67$ years; exponentiating zero gives a hazard ratio of one. This is a model-based crossing, not evidence that a true treatment effect necessarily changes at precisely that time.

Centering Time and Choosing the Reference Point

The treatment main-effect coefficient depends on how time is coded. If $g(t)=t$, then $\beta_0$ is the log hazard ratio at time zero. If time zero is outside the observed follow-up or has no scientific relevance, that coefficient may be awkward to interpret. Centering time sets the reference point at a more useful time without changing the model’s fitted hazard ratios.

Result: For a linear interaction $\beta_0+\beta_1t$, choose a reference time $c$ and write $t-c$ instead. The equivalent coefficients are $\gamma_0=\beta_0+\beta_1c$ and $\gamma_1=\beta_1$, because $$ \beta_0+\beta_1t=(\beta_0+\beta_1c)+\beta_1(t-c). $$ Therefore the hazard ratio at every time is unchanged, while $\gamma_0$ is the log hazard ratio at $t=c$.

This equivalence is an algebraic reparameterization, not a different causal model. Substituting the stated coefficients into the centered expression gives $\beta_0+\beta_1c+\beta_1t-\beta_1c=\beta_0+\beta_1t$ for every $t$. Exponentiating either expression therefore gives the same hazard ratio at every time. $\square$

Worked Example: Centering at a Clinically Relevant Time

Suppose a fitted model uses $g(t)=t$ and has $\hat{\beta}_0=-0.30$ and $\hat{\beta}_1=0.12$ per year. Re-express it using time centered at year 2, and interpret the new main-effect coefficient.

Set $c=2$. Then $\hat{\gamma}_0=\hat{\beta}_0+2\hat{\beta}_1=-0.30+2(0.12)=-0.06$, while $\hat{\gamma}_1=0.12$. The equivalent model’s log hazard ratio is $-0.06+0.12(t-2)$.

At year 2, $t-2=0$, so the hazard ratio is $\exp(-0.06)\approx0.942$. At year 4, it is $\exp[-0.06+0.12(2)]=\exp(0.18)\approx1.197$. Using the original coefficients at year 4 gives $\exp[-0.30+0.12(4)]=\exp(0.18)$, the same result. Centering makes the main effect interpretable at year 2; it does not alter the model’s predictions.

Piecewise Interactions and Contrasts

A single linear interaction assumes the log hazard ratio changes at a constant rate on the chosen time scale. A piecewise specification instead allows the log hazard ratio to differ across prespecified time periods. For example, with a cut point $\tau$, one may use $\beta_0A+\delta A I(t>\tau)$. Before or at $\tau$, the log hazard ratio is $\beta_0$; after $\tau$, it is $\beta_0+\delta$. This model allows a step change at the cut point, rather than a smooth change. The cut point and functional form should be selected for scientific reasons or prespecified analysis plans, not chosen simply because they produce a preferred result.

Worked Example: A Prespecified Step Change

Suppose a model uses $A$ and $A I(t>2)$, with fitted coefficients $\hat{\beta}_0=0.30$ and $\hat{\delta}=-0.50$. Calculate the hazard ratio through year 2 and after year 2.

For $t\leq2$, $I(t>2)=0$. The log hazard ratio is $0.30$, so $\widehat{\operatorname{HR}}(t)=\exp(0.30)\approx1.350$. For $t>2$, the indicator is one. The log hazard ratio is $0.30-0.50=-0.20$, so $\widehat{\operatorname{HR}}(t)=\exp(-0.20)\approx0.819$.

The coefficient $\hat{\delta}$ is the change in the log hazard ratio after the cut point. It is not the post-cut-point log hazard ratio by itself. This parameterization imposes a discontinuity immediately after year 2; it does not describe a gradual transition.

For a general model with treatment terms $A g_1(t),\ldots,A g_K(t)$, write the log hazard ratio as $\beta_0+\sum_{k=1}^K\beta_k g_k(t)$. At a selected time $t$, the log hazard ratio is a linear contrast of the fitted coefficients, and exponentiating that contrast gives the hazard ratio. This makes clear why an interaction coefficient alone is usually not the effect of interest: the relevant quantity combines the coefficients according to the chosen time.

For example, in the linear model the contrast at time $t$ is $\beta_0+t\beta_1$. If $\widehat{\operatorname{Var}}(\hat{\beta})$ is an appropriate estimated covariance matrix for $\hat{\beta}=(\hat{\beta}_0,\hat{\beta}_1)^\mathsf{T}$, then the estimated variance of this log hazard ratio is $c(t)^\mathsf{T}\widehat{\operatorname{Var}}(\hat{\beta})c(t)$, where $c(t)=(1,t)^\mathsf{T}$. A confidence interval on the log scale can be exponentiated to obtain an interval for the hazard ratio. The covariance estimate must suit the weighted analysis; treating weighted records as if they were ordinary independent unweighted observations can give inappropriate uncertainty. The derivation of a suitable variance estimator is a separate issue from interpreting the interaction.

Model Choice, Diagnostics, and Common Mistakes

A constant treatment coefficient is the simpler alternative: it sets the interaction coefficient to zero and yields one hazard ratio for all follow-up times. It is attractive when the constant-effect assumption is scientifically plausible and adequately supported, but can conceal meaningful temporal variation. A prespecified piecewise model allows period-specific effects but assumes abrupt changes at its cut points. A continuous time interaction avoids that jump but imposes the chosen continuous shape. These are different restrictions on the effect function, not different ways to obtain the same answer.

Inspect estimated time-specific hazard ratios at clinically meaningful times, the uncertainty in those contrasts, and whether the fitted pattern is plausible over the observed follow-up. Compare the constant-effect restriction with the interaction model using an appropriate model comparison, while remembering that a failure to detect variation is not proof of constancy. Check that the relevant follow-up periods have adequate weighted risk-set information; sparse late risk sets can make the estimated interaction unstable. Also inspect the treatment and censoring weights and the chosen time coding. Unusually influential weights or extreme uncertainty can signal that the data provide weak support for the contrast, even if an optimizer returns coefficient estimates.

Mistake: calling $\exp(\beta_1)$ the hazard ratio. In a linear interaction, $\beta_1$ is a log-hazard-ratio slope per unit time. The hazard ratio at $t$ is $\exp\{\beta_0+\beta_1t\}$.
Mistake: interpreting the main effect without checking time coding. It is the log hazard ratio where the time function equals zero. Centering changes that reference point.
Mistake: assuming weighting guarantees a valid treatment effect. Weighting addresses measured treatment and censoring processes under its assumptions; it does not repair a wrong time-effect shape, unmeasured confounding, or practical positivity problems.
Mistake: treating a hazard ratio as a risk ratio. A time-specific hazard compares event rates among those still at risk, not cumulative event probabilities.

Finally, a time-varying treatment effect is not the same as a changing treatment history. This model lets the effect of the specified regime vary with follow-up time. It does not, by itself, encode when treatment starts, stops, or changes. Those require a different specification of the treatment history.

Key takeaway: State the time function and its reference point, calculate the log hazard ratio at the time of interest, and exponentiate it. Interpret the resulting hazard ratio as a time-specific marginal hazard contrast under the model and its causal assumptions.

Check Your Understanding

Answer each question, then check the concise solutions below.

A model has log hazard ratio $-0.2+0.1t$. What is its hazard ratio at $t=4$?
In a model using $t-3$, what does the treatment main-effect coefficient represent?
A model uses $0.4A-0.7A I(t>5)$. What are the hazard ratios for $t\leq5$ and $t>5$?
Why is $\exp(\beta_1)$ generally not the time-specific hazard ratio in a model with $A$ and $At$?

Solutions. (1) The log hazard ratio is $-0.2+0.1(4)=0.2$, so the hazard ratio is $\exp(0.2)\approx1.221$. (2) It is the log hazard ratio at $t=3$, where $t-3=0$. (3) Through year 5, the log hazard ratio is $0.4$, giving $\exp(0.4)\approx1.492$; after year 5, it is $0.4-0.7=-0.3$, giving $\exp(-0.3)\approx0.741$. (4) $\beta_1$ is the coefficient multiplying time in the log hazard ratio. At time $t$, the hazard ratio is $\exp(\beta_0+\beta_1t)$, which depends on both coefficients and the selected time.