Constructing Survival Curves from Weighted Cox Fits

From a Weighted Cox Fit to a Survival Curve

A weighted Cox fit estimates treatment effects on a hazard scale. To report the outcome that is often easier to interpret—a probability of surviving beyond each time—we also need a baseline cumulative hazard and a rule for applying the fitted treatment effects to a specified treatment regime. The resulting curve is a model-based estimate for the target population represented by the weighting procedure.

Let $g$ denote a treatment regime, and let $X_g(t)$ be the covariate vector that the marginal structural Cox model assigns to that regime at time $t$. For a time-constant treatment, this vector is constant; for a regime that changes treatment over time, it can change over time. Write the fitted model as

\lambda^g(t)=\lambda_0(t)\exp\{\beta^\mathsf{T}X_g(t)\}.

Here $\lambda^g(t)$ is the hazard under regime $g$, $\lambda_0(t)$ is the baseline hazard at covariate vector zero, and $\beta$ is the treatment-effect parameter. The corresponding survival curve is $S^g(t)=P(T^g>t)$, where $T^g$ is the event time under regime $g$.

Definition: The plug-in regime-specific survival estimate from a weighted Cox fit is $$ \widehat S^g(t)=\exp\left\{-\int_0^t \exp\bigl(\hat\beta^\mathsf{T}X_g(u)\bigr)\,d\widehat\Lambda_0(u)\right\}, $$ where $\hat\beta$ is the fitted coefficient vector and $\widehat\Lambda_0$ is the weighted Breslow estimate of the baseline cumulative hazard. The integral uses the regime’s covariate history, not each participant’s observed treatment history.

This construction relies on more than a fitted regression. For a causal interpretation, consistency, sequential exchangeability, and positivity must support identification of the regime-specific outcome, and the weights must be appropriate for the target population and risk sets. The marginal structural Cox model must adequately describe the regime-specific hazards over the time range of interest. The curve is also subject to the limitations of the fitted model and estimated weights; it is not a direct, assumption-free summary of the observed event times.

The Calculation at Event Times

In the weighted Breslow estimate, the baseline cumulative hazard increases at observed event times. If $t_j$ is an event time, let $Y_i(t_j)$ indicate that individual $i$ is at risk just before $t_j$, $dN_i(t_j)$ be that individual’s event count at $t_j$, and $W_i$ be the analysis weight for the risk interval. In the no-ties case, the increment has the form

\Delta\widehat\Lambda_0(t_j) = \frac{\sum_i W_i\,dN_i(t_j)} {\sum_i W_iY_i(t_j)\exp\{\hat\beta^\mathsf{T}X_i(t_j)\}}.

The numerator counts weighted events; the denominator is the weighted, model-adjusted risk-set total. For any regime, the predicted cumulative hazard through $t$ can then be evaluated as a sum:

\widehat H^g(t) = \sum_{t_j\leq t} \exp\{\hat\beta^\mathsf{T}X_g(t_j)\}\, \Delta\widehat\Lambda_0(t_j), \qquad \widehat S^g(t)=\exp\{-\widehat H^g(t)\}.

The regime-specific multiplier is evaluated at each increment. If the regime is constant, the multiplier is constant and can be taken outside the sum. If it changes, applying one multiplier to the entire baseline cumulative hazard generally gives the wrong result.

Worked Example: A Fixed Treatment Regime

Suppose treatment $A$ is coded 0 or 1, the fitted coefficient is $\hat\beta=\log(0.6)$, and the weighted Breslow baseline increments at times 1, 2, and 3 are $0.10$, $0.20$, and $0.15$. Calculate the fitted survival estimates under each fixed treatment level at time 2.

For the untreated regime, $X_0=0$, so its hazard multiplier is $\exp\{\hat\beta(0)\}=1$. The cumulative hazard through time 2 is

\widehat H^0(2)=1(0.10+0.20)=0.30, \qquad \widehat S^0(2)=\exp(-0.30)\approx 0.7408.

For the treated regime, $X_1=1$, so the multiplier is $\exp\{\log(0.6)\}=0.6$. Thus

\widehat H^1(2)=0.6(0.10+0.20)=0.18, \qquad \widehat S^1(2)=\exp(-0.18)\approx 0.8353.

The calculations use the same baseline increments but apply the fitted hazard multiplier for each regime. The estimates are model-based predictions; their numerical ordering here follows from the fitted coefficient and is not, by itself, evidence that treatment causally improves survival.

Worked Example: A Regime That Changes Treatment

Suppose $A(t)$ is coded 0 or 1, $\hat\beta=\log(2)$, and the baseline increments at times $0.5$, $1.5$, and $2.5$ are $0.08$, $0.12$, and $0.10$. Consider a regime that assigns $A=0$ through time 1 and $A=1$ after time 1. Calculate its survival estimate at each event time.

The treatment multiplier is 1 at time $0.5$, then 2 at times $1.5$ and $2.5$. At $0.5$,

\widehat H^g(0.5)=1(0.08)=0.08, \qquad \widehat S^g(0.5)=\exp(-0.08)\approx0.9231.

At $1.5$, add the increment weighted by the new regime value:

\widehat H^g(1.5)=0.08+2(0.12)=0.32, \qquad \widehat S^g(1.5)=\exp(-0.32)\approx0.7261.

At $2.5$, the total is

\widehat H^g(2.5)=0.08+2(0.12)+2(0.10)=0.52, \qquad \widehat S^g(2.5)=\exp(-0.52)\approx0.5945.

Using the treated multiplier of 2 on all three increments would instead give $2(0.08+0.12+0.10)=0.60$. That calculation imposes treatment from the beginning, so it does not represent the stated regime.

Why the Construction Works

For a regime with a hazard $\lambda^g(t)$, the cumulative hazard is $H^g(t)=\int_0^t\lambda^g(u)\,du$. When the hazard is continuous in time, survival solves $dS^g(t)/dt=-S^g(t)\lambda^g(t)$, with $S^g(0)=1$. Separating variables and integrating gives the exponential relationship used for prediction.

Derivation: Under an absolutely continuous hazard, $$ \frac{dS^g(t)}{S^g(t)}=-\lambda^g(t)\,dt. $$ Integrating from 0 to $t$, using $S^g(0)=1$, yields $$ \log S^g(t)=-\int_0^t\lambda^g(u)\,du, \qquad S^g(t)=\exp\left\{-\int_0^t\lambda^g(u)\,du\right\}. $$ Substituting the marginal structural Cox model gives $$ S^g(t)=\exp\left\{-\int_0^t\exp\bigl(\beta^\mathsf{T}X_g(u)\bigr)\,d\Lambda_0(u)\right\}. $$

The fitted curve replaces $\beta$ and $\Lambda_0$ with their estimates. Since the weighted Breslow estimate is a step function, its integral is evaluated as a sum over its jumps. This yields the event-time calculation above.

Two useful properties follow directly from that calculation. If the estimated baseline increments are nonnegative, each regime multiplier is positive, so the predicted cumulative hazard is nondecreasing. Exponentiating its negative therefore gives a survival estimate in $(0,1]$ that cannot increase over time.

Result: If $\Delta\widehat\Lambda_0(t_j)\geq0$ and $\hat\beta$ and $X_g(t_j)$ are finite, then $\widehat H^g(t)$ is nondecreasing and $\widehat S^g(t)$ is nonincreasing, with $0<\widehat S^g(t)\leq1$.

To verify the result, each summand $\exp\{\hat\beta^\mathsf{T}X_g(t_j)\}\Delta\widehat\Lambda_0(t_j)$ is nonnegative. Extending the sum to a later time can only add nonnegative terms, so $\widehat H^g(t)$ cannot decrease. The function $h\mapsto\exp(-h)$ is positive and decreasing, and $\widehat H^g(0)=0$, giving the stated bounds and monotonicity. $\square$

This is a structural check on the calculation, not a guarantee that the model is well specified or the causal assumptions hold. For example, a curve can be monotone and still be based on unstable weights, an incorrect regime coding, or a poor proportional-hazards fit.

Worked Example: Baseline Centering and Regime Prediction

Suppose software reports a baseline cumulative hazard $\widehat\Lambda_{x_{\mathrm{ref}}}$ for a reference covariate value $x_{\mathrm{ref}}=1$, rather than at zero. Let $\hat\beta=\log(2)$, and suppose the reported baseline increments are $0.04$ and $0.06$. Find the survival estimate at the second event time for the regime $x_g=0$.

The reported baseline is for linear predictor $\hat\beta x_{\mathrm{ref}}$. For $x_g=0$, the hazard relative to that reference is multiplied by $\exp\{\hat\beta(x_g-x_{\mathrm{ref}})\}=\exp\{-\log(2)\}=0.5$.

\widehat H^g=0.5(0.04+0.06)=0.05, \qquad \widehat S^g=\exp(-0.05)\approx0.9512.

Applying the reference curve without transforming it would give $\exp(-0.10)\approx0.9048$, which predicts for $x=1$, not $x=0$. Before using a software-produced baseline, check which covariate value it represents and whether its reported quantity is cumulative hazard or survival.

Diagnostics, Alternatives, and Limits

A practical check is to compute the curve from the event-time increments directly for at least one time point, then compare it with the software prediction. Check that the regime vector is coded in the same way as the fitted model, that each interval uses the treatment history assigned by the regime, and that the baseline estimate corresponds to the intended reference covariate value. Also inspect the risk sets and weights used to estimate the baseline: very large or highly variable weights can make estimates sensitive to a few participants. Weight stabilization or truncation changes the estimating procedure and target interpretation and should be reported, not treated as a harmless numerical adjustment.

A common mistake is to calculate $\exp\{-\widehat\Lambda_0(t)\}$ and call it the treatment-specific survival curve. That expression is the curve at $X=0$, the baseline covariate value; it omits the fitted treatment multiplier. A second mistake is to multiply the entire baseline cumulative hazard by one treatment effect for a regime whose treatment changes over time. The correction is to apply the regime-specific multiplier separately to every baseline increment.

A weighted Kaplan–Meier estimator can provide a nonparametric or less model-structured survival summary under suitable weighting and censoring conditions, but it does not use the Cox model’s treatment-effect structure to predict a specified regime in the same way. An unweighted Cox model estimates associations under its own model and does not, by itself, correct treatment-confounding through weighting. A parametric survival model can provide smoother curves, but depends on additional distributional assumptions. The weighted Cox construction is useful when the marginal structural proportional-hazards model is scientifically defensible; it does not remove the need to assess that assumption, censoring, positivity, weight behavior, or the time horizon supported by the data.

Check Your Understanding

Try each calculation or explanation, then check the solutions below.

For a fixed regime with hazard multiplier 1.5 and baseline increments 0.04 and 0.06, calculate the fitted survival estimate after the second increment.
A regime has multiplier 1 for the first baseline increment, then multiplier 0.5 for the next two. The increments are 0.10, 0.20, and 0.10. Find its cumulative hazard and survival estimate after the third increment.
Why is a predicted survival curve that increases between event times a warning sign for the calculation described here?
Software reports a baseline curve at $x_{\mathrm{ref}}$, but the target regime has a different fitted linear predictor. What adjustment is needed before interpreting the curve for that regime?

Solutions. (1) The cumulative hazard is $1.5(0.04+0.06)=0.15$, so survival is $\exp(-0.15)\approx0.8607$. (2) The cumulative hazard is $1(0.10)+0.5(0.20)+0.5(0.10)=0.25$, so survival is $\exp(-0.25)\approx0.7788$. (3) With nonnegative Breslow increments and positive regime multipliers, cumulative hazard cannot decrease, so its negative exponential cannot increase. Check for a sign error, incorrect ordering or subtraction of increments, or a plotting/interpolation issue. (4) Multiply the reported baseline increments by the exponential of the difference between the regime and reference linear predictors, $\exp\{\hat\beta^\mathsf{T}(X_g-x_{\mathrm{ref}})\}$, before accumulating the regime hazard; then exponentiate the negative cumulative hazard to obtain survival.

Constructing Survival Curves from Weighted Cox Fits

What You'll Learn