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 Surv.counting() wraps it into the response object.
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")
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 the per-subject baseline table and the repeated-measurements table and merges them into the counting-process long format. The id column links the two. visit_time names the column in visits that carries measurement times.
import pandas as pd
# Build a per-subject baseline: one row per subject with follow-up time and outcome.
# pbcseq repeats futime/status for every visit, so drop_duplicates gives one row per id.
pbcseq_pd = gw.load_dataset("pbcseq", backend="pandas")
base = (
pbcseq_pd.drop_duplicates("id")[["id", "futime", "status"]]
.rename(columns={"futime": "time"})
)
# Split into counting-process intervals, carrying bilirubin, albumin, and protime forward.
long = gw.split_episodes(
base,
pbcseq_pd[["id", "day", "bili", "albumin", "protime"]],
id="id",
time="time",
event="status",
visit_time="day",
format="pandas",
)
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 |
Before fitting we drop the small number of rows with missing covariate values and create a binary event indicator.
long_clean = long.dropna(subset=["bili", "albumin", "protime"]).copy()
long_clean["event"] = (long_clean["status"] == 2).astype(int)
print(f"{len(long_clean)} intervals, {long_clean['event'].sum()} events, "
f"{long_clean['id'].nunique()} subjects")
1945 intervals, 140 events, 312 subjects
Fitting the TVC Cox model
The counting-process Surv object wraps tstart and tstop alongside the event indicator. Pass the resulting response and covariate columns to CoxPH exactly as with a right-censored model.
y = gw.Surv.counting(long_clean["tstart"], long_clean["tstop"], long_clean["event"])
cox_tvc = gw.CoxPH().fit(y, long_clean[["bili", "albumin", "protime"]])
cox_tvc
/var/folders/s8/bj_jsx3d7jqd2btw7bwm6yx80000gp/T/ipykernel_11628/3376967461.py:2: UserWarning: Subjects in counting-process data have different start times, some much larger than 0. This may indicate that start/stop times are calendar time rather than subject-relative time. Each subject's timeline should begin at 0. If you have calendar dates, subtract each subject's entry date from their start/stop times before fitting.
cox_tvc = gw.CoxPH().fit(y, long_clean[["bili", "albumin", "protime"]])
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.
pbc_pd = gw.load_dataset("pbc", backend="pandas")
pbc_clean = pbc_pd.dropna(subset=["bili", "albumin", "protime"]).copy()
pbc_clean["event"] = (pbc_clean["status"] == 2).astype(int)
y_static = gw.Surv.right(pbc_clean["time"], event=pbc_clean["event"])
cox_static = gw.CoxPH().fit(y_static, pbc_clean[["bili", "albumin", "protime"]])
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
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.
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.56241112508 |
| 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]
})
surv_static = cox_tvc.predict(
static_path,
type="survival",
times=[0, 50, 100, 192, 300, 400],
format="pandas",
)
# Side-by-side comparison
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.56241112508 |
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
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")
)
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(base, visits, id=, time=, event=, visit_time=) |
| Wrap the response |
gw.Surv.counting(tstart, tstop, event) |
| Fit the TVC Cox model |
gw.CoxPH().fit(y, covariates) |
| 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.