A standard Cox model assumes each covariate is measured once and stays constant throughout follow-up. In many studies that assumption is wrong. A patient’s bilirubin rises over months, a worker’s cumulative exposure increases with each shift, a trial subject crosses over to the active arm mid-study. Treating a changing value as fixed introduces bias in the estimated hazard ratio. Time-varying covariates (TVC) solve this by letting the covariate take a different value in each interval of follow-up.
The counting-process form of the Cox model handles time-varying covariates naturally. Instead of one row per subject you create one row per interval: each row records when the subject entered the risk set (tstart), when they exited (tstop), the event indicator for that interval, and the covariate values that were constant during it. The split_episodes() function builds this dataset from the raw repeated-measurement data, and the formula response Surv(tstart, tstop, status) describes it, just as in R.
The data
The pbcseq dataset contains repeated lab measurements for subjects in the Mayo Clinic PBC trial. It is a companion to the cross-sectional pbc dataset: while pbc gives each subject a single row with final outcome, pbcseq records the values of bilirubin, albumin, prothrombin time, and several other markers at each clinical visit.
import greenwood as gw
# Load both datasets as Polars DataFrames
pbc = gw.load_dataset("pbc", backend="polars")
pbcseq = gw.load_dataset("pbcseq", backend="polars")
# Preview the repeated lab measurements
pbcseq.select(["id", "day", "bili", "albumin", "protime"])
PolarsRows1,945Columns5 |
|
|
|
|
|
|
| 0 |
1 |
0 |
14.5 |
2.6 |
12.2 |
| 1 |
1 |
192 |
21.3 |
2.94 |
11.2 |
| 2 |
2 |
0 |
1.1 |
4.14 |
10.6 |
| 3 |
2 |
182 |
0.8 |
3.6 |
11 |
| 4 |
2 |
365 |
1 |
3.55 |
11.6 |
| 5 |
2 |
768 |
1.9 |
3.92 |
10.6 |
| 6 |
2 |
1790 |
2.6 |
3.32 |
11.3 |
| 7 |
2 |
2151 |
3.6 |
2.92 |
11.5 |
| 1942 |
312 |
390 |
7.4 |
3.56 |
11.7 |
| 1943 |
312 |
775 |
16.3 |
3.34 |
13 |
| 1944 |
312 |
1075 |
23.4 |
3.42 |
13.4 |
Each subject appears once per visit. day is the time of that visit (in days from enrollment), and the other columns hold the lab values measured that day. futime and status record the final outcome and are the same for every row belonging to the same subject.
Building the counting-process dataset
split_episodes() takes a per-subject baseline table and a repeated-measurements table and merges them into the counting-process long format. The id column links the two, and visit_time names the column in visits that carries measurement times. covariates= chooses which columns of visits become time-varying covariates.
A visit exactly at the end of follow-up does not create a zero-length interval. The preceding interval still carries the subject’s event indicator and uses the covariate value measured before that endpoint. If the only visit is at the endpoint, carry_forward=True keeps (0, time] with missing covariates and the original event.
pbcseq already repeats each subject’s follow-up time and status on every visit row, so the same table can serve as both inputs. split_episodes() then uses the first row of each subject as its baseline, after checking that futime and status really are constant within the subject.
import pandas as pd
# One long table: one row per visit, with follow-up and status repeated for each subject
pbcseq_pd = gw.load_dataset("pbcseq", backend="pandas")
# Split into counting-process intervals, carrying bilirubin, albumin, and protime forward
long = gw.split_episodes(
baseline=pbcseq_pd,
visits=pbcseq_pd,
id="id",
time="futime",
event="status",
visit_time="day",
covariates=["bili", "albumin", "protime"],
format="pandas",
)
# Preview the counting-process dataset
long
PandasRows1,945Columns7 |
|
|
|
|
|
|
|
|
| 0 |
1 |
0 |
192 |
0 |
14.5 |
2.6 |
12.2 |
| 1 |
1 |
192 |
400 |
2 |
21.3 |
2.94 |
11.2 |
| 2 |
2 |
0 |
182 |
0 |
1.1 |
4.14 |
10.6 |
| 3 |
2 |
182 |
365 |
0 |
0.8 |
3.6 |
11 |
| 4 |
2 |
365 |
768 |
0 |
1 |
3.55 |
11.6 |
| 5 |
2 |
768 |
1790 |
0 |
1.9 |
3.92 |
10.6 |
| 6 |
2 |
1790 |
2151 |
0 |
2.6 |
3.32 |
11.3 |
| 7 |
2 |
2151 |
2515 |
0 |
3.6 |
2.92 |
11.5 |
| 1942 |
312 |
390 |
775 |
0 |
7.4 |
3.56 |
11.7 |
| 1943 |
312 |
775 |
1075 |
0 |
16.3 |
3.34 |
13 |
| 1944 |
312 |
1075 |
1457 |
0 |
23.4 |
3.42 |
13.4 |
Each row is now one interval for one subject. Subject 1, for example, contributes two rows: the first covers days 0–192 (bilirubin 14.5), and the second covers days 192–400 (bilirubin 21.3). The event column is 0 for the first row and carries the original status code (2 = death) for the last.
# Subject 1's intervals
long[long["id"] == 1][["id", "tstart", "tstop", "status", "bili", "albumin", "protime"]]
PandasRows2Columns7 |
|
|
|
|
|
|
|
|
| 0 |
1 |
0 |
192 |
0 |
14.5 |
2.6 |
12.2 |
| 1 |
1 |
192 |
400 |
2 |
21.3 |
2.94 |
11.2 |
A few intervals have missing lab values, and status uses the pbc coding (2 is death). Neither needs a manual cleaning step. The comparison status == 2 in the formula marks deaths, and when the model is fit, every row with a missing value in a column it uses is dropped. The printout reports how many.
Fitting the TVC Cox model
The counting-process response names tstart and tstop before the event, as R’s Surv(start, stop, event) does. Pass the formula to CoxPH with the frame, exactly as with a right-censored model.
# A status of 2 marks a death at the end of an interval, and rows with missing lab values are
# dropped at fit time
cox_tvc = gw.CoxPH().fit(
"Surv(tstart, tstop, status == 2) ~ bili + albumin + protime", data=long
)
cox_tvc
CoxPH (Cox proportional hazards model, ties='efron')
coef exp(coef) se(coef) z p
bili 0.1378 1.148 0.01044 13.203 8.467e-40
albumin -2.016 0.1331 0.1649 -12.232 2.107e-34
protime 0.1962 1.217 0.03393 5.783 7.356e-09
n = 1945, events = 140
Likelihood ratio test = 440.4 on 3 df, p = 3.852e-95
For comparison, fit the same model using only the baseline lab values from pbc.
# Fit the same covariates using only baseline values from the cross-sectional dataset
pbc = gw.load_dataset("pbc", backend="pandas")
# The same death endpoint, now with one row per subject
cox_static = gw.CoxPH().fit("Surv(time, status == 2) ~ bili + albumin + protime", data=pbc)
cox_static
CoxPH (Cox proportional hazards model, ties='efron')
coef exp(coef) se(coef) z p
bili 0.1202 1.128 0.0126 9.547 1.332e-21
albumin -1.226 0.2933 0.198 -6.193 5.902e-10
protime 0.1864 1.205 0.04962 3.756 0.0001729
n = 416, events = 160 (2 observations deleted due to missingness)
Likelihood ratio test = 146.5 on 3 df, p = 1.481e-31
The TVC model gives notably larger hazard ratios for bilirubin (exp(0.14) ≈ 1.15 vs exp(0.09) ≈ 1.09 in the static model). Updating the covariate at each visit lets the model use the most recent measurement, which better reflects the current disease severity and recovers more of the true effect.
Predicting survival along a covariate trajectory
For a subject whose covariates change over time you cannot describe them with a single static vector. Instead, pass a trajectory DataFrame to predict(). Each row of the trajectory defines a constant-covariate interval (tstart, tstop]. The covariate values in that interval are used to accumulate hazard up to any requested time point.
Subject 1 of the PBC study had two visits. Their trajectory therefore has two rows.
# Build subject 1's two-visit covariate trajectory
tvc_path = pd.DataFrame({
"tstart": [0.0, 192.0],
"tstop": [192.0, 400.0],
"bili": [14.5, 21.3],
"albumin": [2.60, 2.94],
"protime": [12.2, 11.2],
})
tvc_path
PandasRows2Columns5 |
|
|
|
|
|
|
| 0 |
0 |
192 |
14.5 |
2.6 |
12.2 |
| 1 |
192 |
400 |
21.3 |
2.94 |
11.2 |
# Predict survival at selected time points along the trajectory
surv_tvc = cox_tvc.predict(
trajectory=tvc_path,
times=[0, 50, 100, 192, 300, 400],
format="pandas",
)
surv_tvc
PandasRows6Columns2 |
|
|
|
| 0 |
0 |
1 |
| 1 |
50 |
0.966310651788 |
| 2 |
100 |
0.851568181027 |
| 3 |
192 |
0.667410151188 |
| 4 |
300 |
0.562411125081 |
| 5 |
400 |
0.395677691279 |
Notice the kink near day 192: bilirubin rises from 14.5 to 21.3 at the second visit, which steepens the hazard and causes survival to drop more quickly after that point. A static model using only the baseline value of 14.5 would miss this deterioration entirely.
To see how different the trajectories are, compare the TVC prediction to a static prediction that fixes bilirubin at 14.5 for the entire follow-up:
# Static prediction: fix covariates at their baseline (day-0) values
static_path = pd.DataFrame({
"bili": [14.5], "albumin": [2.60], "protime": [12.2]
})
# Predict survival with covariates fixed at baseline values
surv_static = cox_tvc.predict(
static_path,
type="survival",
times=[0, 50, 100, 192, 300, 400],
format="pandas",
)
# Side-by-side comparison of TVC vs. static predictions
import polars as pl
comparison = pl.from_pandas(surv_tvc).rename({"subject_1": "tvc"}).with_columns(
pl.from_pandas(surv_static)["subject_1"].alias("static (baseline bili only)")
)
comparison
PolarsRows6Columns3 |
|
|
|
static (baseline bili only) f64 |
| 0 |
0 |
1 |
1 |
| 1 |
50 |
0.966310651788 |
0.966310651788 |
| 2 |
100 |
0.851568181027 |
0.851568181027 |
| 3 |
192 |
0.667410151188 |
0.667410151188 |
| 4 |
300 |
0.562411125081 |
0.567617113732 |
| 5 |
400 |
0.395677691279 |
0.406971009427 |
At day 400 (the subject’s actual exit time) the TVC prediction gives a survival probability roughly 10 percentage points lower than the static prediction, reflecting the bilirubin rise captured during follow-up.
Visualizing the survival trajectory
Build predictions at a finer time resolution and melt them into a tidy long frame, then pass it to Altair to draw overlapping survival curves.
import altair as alt
import polars as pl
# Predict at fine resolution for smooth plotting
times = list(range(0, 421, 5))
surv_tvc_fine = cox_tvc.predict(
trajectory=tvc_path, times=times, format="polars"
)
surv_static_fine = cox_tvc.predict(
static_path, type="survival", times=times, format="polars"
)
# Combine into a tidy long frame for Altair
plot_data = (
surv_tvc_fine.rename({"subject_1": "TVC (updating bilirubin)"})
.join(
surv_static_fine.rename({"subject_1": "Static (baseline bilirubin)"}),
on="time",
)
.unpivot(index="time", variable_name="model", value_name="survival")
)
# Draw the two survival curves
alt.Chart(plot_data).mark_line().encode(
x=alt.X("time:Q", title="Time (days)"),
y=alt.Y("survival:Q", title="Survival probability", scale=alt.Scale(domain=[0, 1])),
color=alt.Color("model:N", title="Model"),
)
The TVC curve (updating bilirubin at day 192) separates noticeably from the static curve after the second visit. Both start from the same point because the baseline covariate values are identical. The divergence begins exactly where the trajectory changes.
Summary
| Load repeated-measurement data |
gw.load_dataset("pbcseq") |
| Build counting-process dataset |
gw.split_episodes(baseline=, visits=, id=, time=, event=, visit_time=, covariates=[...]) |
| Fit the TVC Cox model |
gw.CoxPH().fit("Surv(tstart, tstop, status == 2) ~ bili + ...", data=long) |
| Predict along a trajectory |
cox.predict(trajectory=path, times=...) |
Time-varying covariates require more data preparation than a static model, but the tradeoff is hazard ratios that reflect the covariate at the time of each event rather than a single snapshot.