ExpertIndependently reviewed, not yet spot-checked by a human

Time-dependent covariates

When exposure changes during follow-up, Cox needs a (start, stop] long table rather than one row per person. How tmerge() builds it, why an Andersen-Gill model for recurrent events must use a robust variance, whether a repeated lab value should enter as its baseline or its current reading, when a landmark analysis is the better choice, and why a time-dependent hazard ratio cannot be used for prediction.

Two things that get conflated

They share a name and they are entirely different problems:

	Time-dependent covariate	Time-varying coefficient
What moves	$x$ moves: this person’s exposure status today differs from three months from now	$\beta$ moves: exposure is fixed, but its effect grows or decays over time
The model	$h(t \mid x) = h_0(t) \exp(\beta \, x(t))$	$h(t \mid x) = h_0(t) \exp(\beta(t) \, x)$
Typical examples	Receiving a transplant, starting dialysis, today’s serum bilirubin, current medication	Early surgical hazard, delayed immunotherapy effect, waning vaccine protection
What it solves	“The value at baseline does not represent the whole of follow-up”	“Proportional hazards has been violated”
Where it lives	This page	B3-04

The (start, stop] long table: one person, many rows

Ordinary survival data gives each person one row: a time, an event indicator, and some fixed covariates. Once exposure moves, that shape is not enough, because “what is this person’s $x$ ” has different answers at different moments.

The fix is to cut each person into segments inside which $x$ is constant. One row per segment, with columns for the half-open interval (tstart, tstop], whether an event occurred at the end of that segment, and the covariate values that applied during it.

Two timeline panels side by side. On the left, six patients from the Stanford heart transplant study, one horizontal line each; a patient who receives a transplant changes colour at that moment, splitting the line into two segments, and a red dot at the end of a line marks death. On the right, six patients with primary biliary cholangitis, each line cut into many short segments separated by thin white gaps, with the shade of each segment showing which tertile of the distribution that period's serum bilirubin fell into. — Left: an exposure that switches at most once, so a person has at most two rows. Right: a lab value measured at scheduled visits, so a new segment starts at every visit — one person can easily have a dozen rows. Both feed the model exactly the same data structure.Plotting script `figures/scripts/B3-06-time-varying.R`

What matters is that this structure does not turn one person into several people. Cox estimation happens entirely within the risk set at each event time (see B3-01): at time $t$ the model only asks “of the people still under follow-up right now, what is each one’s $x$ at this instant”. Every row of the long table is just a lookup entry answering that question, and a person has exactly one row covering time $t$ , so they are counted once.

`tmerge()`: turning a wide table into a long one

Cutting intervals by hand goes wrong easily — overlaps, missing segments, an event recorded on the wrong row — and tmerge() in the survival package exists to do exactly this. Its mental model: first build a skeleton holding only the follow-up interval and the outcome, then layer events and time-dependent variables onto it one at a time.

Three main verbs:

Function	Meaning	Used for
`event(time, status)`	Record an event at `time` and cut follow-up there	The primary outcome
`tdc(time, value)`	From `time` onward the covariate becomes `value` (time-dependent covariate)	Lab values, medications, surgery
`cumevent(time)`	A running count, incremented from `time` onward	“How many have happened so far”

Run it on the Stanford heart transplant data (103 patients on the waiting list, 69 of whom received a transplant) and you get 170 rows: one each for the 34 who were never transplanted, two each for the 67 who were, and then 2 patients transplanted on day 0 — tdc() cannot cut an interval at time zero, so they get a single row too, and that row already carries transplant equal to 1.

34 + 2 = 36 patients with one row, plus 67 with two rows, is exactly the 103 patients in the study; and 67 + 2 = 69 received a transplant. Do both sums every time you build a long table. The common misuses table below says that a badly built long table throws no error and simply returns the wrong answer quietly — those two day-0 patients are that sentence made concrete.

One critical detail about tdc(): by default the covariate takes effect only after the stated moment. A patient transplanted today still counts as untransplanted at that instant. The convention is deliberate. If a change and an event were allowed to take effect simultaneously, and the event were itself the reason for the change (“deteriorating, so operated on emergently, died the same day”), you would be smuggling outcome information into the exposure.

The same logic applies to lab values: blood drawn today cannot be used to predict an event that happens today. One of the most common contaminations in real data is the deranged lab value from the day before death being recorded as “the covariate for that period”, after which the model appears to predict superbly when all it has read is the fact of dying. The remedy is to lag the lab value — by a day, by a week — and run it as a sensitivity analysis.

Repeated measurements: baseline value or current value

pbcseq is the follow-up data from a primary biliary cholangitis cohort: 312 patients with scheduled blood tests, a median of 5 visits and a maximum of 16, followed for up to 12.5 years, during which 125 died. Reshaped into a long table it becomes 1807 rows.

Bilirubin is the single most important prognostic marker in this disease, and it rises as the disease advances. Using the baseline measurement and using “the most recent value as of right now” give appreciably different answers:

Covariate	Baseline value: HR	95% CI	Current value: HR	95% CI
log(bilirubin)	2.920	2.387–3.572	4.004	3.169–5.058
log(albumin)	0.032	0.008–0.122	0.008	0.003–0.021
Age at entry	1.040	1.023–1.057	1.047	1.028–1.067
Concordance (C statistic)	0.833		0.910

Two panels side by side. On the left, serum bilirubin trajectories on a log scale for 24 randomly selected patients with primary biliary cholangitis: red lines are patients who died during follow-up and mostly trend upward, climbing steeply near the end, while pale blue lines are patients still alive and stay roughly flat. On the right, two horizontal confidence intervals: the upper one is the hazard ratio of 2.920 using baseline bilirubin only, the lower one is 4.004 using the current value, and the intervals overlap slightly between 3.169 and 3.572 while the point estimate for the current value falls outside the baseline interval. — Left: bilirubin does not stay where it started, and the trajectory shapes differ between those who died and those who did not. Right: same patients, same outcome, one change in how the covariate is defined, and the hazard ratio for log bilirubin moves from 2.920 to 4.004.Plotting script `figures/scripts/B3-06-time-varying.R`

The hazard ratio from baseline bilirubin is the smaller one, and there is no mystery about why: the baseline value is a proxy for the current value that gradually stops working. The longer follow-up runs, the further that one early blood draw sits from the patient’s condition right now, and the association is diluted — regression dilution. Using the current value removes that layer of dilution.

Recurrent events: Andersen-Gill and the robust variance

Some outcomes happen again and again: infections, asthma attacks, heart failure readmissions, seizures. One row per person keeps only the first and discards the rest.

cgd is a randomised trial in chronic granulomatous disease comparing interferon gamma with placebo for the prevention of serious infection. 128 children, 76 infections in total during follow-up, 17 of the children having more than one and the worst-affected child having 7.

Keeping only the first event throws away 32 events, which is 42.1% of the total.

The Andersen-Gill (AG) model takes an unglamorous approach: cut each person’s follow-up into intervals running from one event to the next, throw them all into a single Cox model, and treat them as separate observations. The data structure is exactly the long table from the previous section.

Analysis	Events used	HR (interferon vs placebo)	95% CI	SE of the coefficient
First event only	44	0.335	0.174–0.645	0.335
Andersen-Gill, no robust variance	76	0.334 (same as the row below)	0.201–0.558 (too narrow)	0.261 (understated)
Andersen-Gill with a robust (cluster) variance	76	0.334	0.181–0.616	0.312

The default variance of an AG model is wrong, and must be replaced by a robust one. The standard Cox variance assumes every row is an independent observation, but repeated infections in the same child are obviously correlated — an infection-prone child keeps getting infected. Ignoring that amounts to claiming those 203 rows are 203 independent observations, when they come from 128 children and the 76 infections occurred in only 44 of them.

The cost is visible: the standard error goes from 0.261 to 0.312, inflating by 19.5%. In R this is one added cluster(id) (as + cluster(id) or cluster = id); in SAS it is covs(aggregate). Forget it and your confidence intervals are systematically too narrow.

Note that even after paying that cost, the AG standard error (0.312) is still slightly smaller than the first-event-only analysis (0.335) — the extra 32 events did buy some precision. But the gain is far smaller than “nearly twice as many events” would suggest, because information from repeated events is discounted to begin with.

AG is not the only recurrent event model, and the choice turns on who belongs in the risk set for the second event.

Model	Assumption	Suits
Andersen-Gill	The baseline hazard is the same for every event and event order carries no effect	Events that are independent, where having one does not change the propensity for the next (accidental injury, say)
PWP (Prentice-Williams-Peterson)	Stratified by event number; the $k$ -th event has its own baseline hazard	Conditions where having had one makes another likelier (most chronic disease progression)
WLW (Wei-Lin-Weissfeld)	Every event is timed from time zero	Events of different kinds, each with its own effect estimate
Frailty models	Each person carries a random effect (their frailty)	When you want to estimate the variation in individual susceptibility explicitly

This is a substantive modelling choice, driven by clinical mechanism and written down in advance. A paper that says only “a Cox model for recurrent events was used” leaves the reader with nothing to judge.

Landmark analysis: the blunter option that is easier to explain

A time-dependent model is not the only way out. Landmark analysis picks a timepoint $L$ , keeps only the people who survived to $L$ , groups them by their status at $L$ , and restarts the clock there.

Why that avoids immortal time bias: everything needed to assign groups is already known at $L$ , and time before $L$ is excluded from both groups alike.

Using the Stanford heart transplant data with a landmark at day 30:

Two panels side by side. On the left, two Kaplan-Meier curves after the day 30 landmark: a blue curve for patients not yet transplanted and a green curve for those already transplanted, crossing back and forth throughout follow-up without clearly separating. On the right, hazard ratios with confidence intervals from three landmark timepoints (days 14, 30 and 60); the point estimate moves from above 1 to below 1 and all three intervals cross 1. — Left: the 79 patients who survived to day 30, split by whether they had been transplanted by then. Right: change the landmark and the point estimate moves with it — not merely statistical noise, since the three analyses define three different populations and three different groupings.Plotting script `figures/scripts/B3-06-time-varying.R`

Landmark	Kept	Dropped	Transplanted by then	HR	95% CI	p
Day 14	88	15	21	1.441	0.812–2.558	0.212
Day 30	79	24	35	0.915	0.527–1.589	0.752
Day 60	64	39	45	0.808	0.407–1.604	0.542

None of the three landmarks reaches statistical significance — every confidence interval crosses 1 — which agrees with the time-dependent model in B3-01. What agrees is the absence of a detected difference, not the direction of the point estimate: B3-01’s hazard ratio is above 1, whereas two of these three are below it. And the point estimate travels from 1.441 at day 14 to 0.808 at day 60, while the number of patients retained falls from 88 to 64.

	Time-dependent covariate	Landmark analysis
Data used	All of it	Only those alive at $L$ , and only follow-up after $L$
Group definition	Updated at every instant	Frozen at $L$ (anyone who switches later is misclassified)
Statistical efficiency	Higher	Lower, and the later $L$ is, the more is discarded
Interpretability for readers	Poorer — the hazard ratio matches no baseline grouping	Better — it is just two groups’ KM curves
Main risk	Data-wrangling errors (overlapping intervals, covariates peeking at the future)	The choice of $L$ ; picking $L$ after the fact is data dredging

An interpretive limit: a time-dependent hazard ratio cannot be used for prediction

This is the most important section on the page, and the one most often skipped.

The hazard ratio for a time-dependent covariate is a contemporaneous association: it says that at time $t$ , people whose $x(t)$ is high have a higher hazard at that instant. It is not a predictive statement.

The reason: to use this model to compute someone’s probability of surviving the next five years, you would need to know their $x(t)$ at every instant of those five years — which is exactly what you do not have. How bilirubin will move is itself unknown; if you knew it, you would probably already know the outcome.

Three concrete consequences:

Do not use a time-dependent model to produce survival curves or $n$ -year survival probabilities. survfit() applied to one will return something, but that curve describes a hypothetical person whose covariates never move, which is usually no real patient.
Do not use a time-dependent model’s C statistic to claim predictive performance. The rise from 0.833 to 0.910 above comes from the model having seen information generated after the moment of prediction. To compare predictive performance, both models may use only what was known at the moment they predict.
To build a prediction model around a marker that moves, use a different tool: landmark prediction (fitting a separate model at each landmark using only information available then), or a joint model, which models the longitudinal biomarker trajectory and survival time together (the JM and JMbayes2 packages).

There is also a subtler problem: if the time-dependent covariate is itself affected by the exposure and is also a cause of the outcome — antihypertensive drug → blood pressure → stroke — then putting it in a Cox model both blocks part of the real effect and opens a new biasing path. This “time-dependent mediator that is also a confounder” cannot be handled by standard Cox regression; the correct tool is a marginal structural model with IPTW, covered in B6-08.

Run it yourself

library(survival)

# -- 1. What a long table looks like --
data(cgd, package = "survival")
head(cgd[cgd$id %in% c(1, 2), c("id", "tstart", "tstop", "status", "treat", "enum")])

# -- 2. tmerge(): building a long table from a wide one --
# Skeleton: one row per person, the follow-up interval and the outcome
base <- subset(pbc, id <= 312, select = c(id:sex, stage))
d <- tmerge(base, base, id = id, death = event(time, status == 2))
# Layer on the time-dependent labs; tdc() takes effect AFTER the stated moment
d <- tmerge(d, pbcseq, id = id,
            bili = tdc(day, bili), albumin = tdc(day, albumin))
head(d[d$id == 2, c("id", "tstart", "tstop", "death", "bili", "albumin")])

# -- 3. Baseline value vs current value --
b <- merge(base, pbcseq[!duplicated(pbcseq$id), c("id", "bili", "albumin")], by = "id")
coxph(Surv(time, status == 2) ~ log(bili) + log(albumin) + age, data = b)
coxph(Surv(tstart, tstop, death) ~ log(bili) + log(albumin) + age, data = d)

# -- 4. Recurrent events: Andersen-Gill --
coxph(Surv(tstart, tstop, status) ~ treat, data = cgd)                  # SE is wrong
coxph(Surv(tstart, tstop, status) ~ treat + cluster(id), data = cgd)    # robust variance
# PWP: stratify by event number, so the k-th event has its own baseline hazard
coxph(Surv(tstart, tstop, status) ~ treat + strata(enum) + cluster(id), data = cgd)

# -- 5. Landmark analysis --
data(heart, package = "survival")
L <- 30
lm <- subset(jasa, futime >= L)
lm$g <- factor(ifelse(!is.na(lm$wait.time) & lm$wait.time <= L, "tx", "no"))
coxph(Surv(futime - L, fustat) ~ g, data = lm)
plot(survfit(Surv(futime - L, fustat) ~ g, data = lm))

# -- 6. Check the reshaping for errors (never skip this) --
# Intervals must not overlap, must not have tstop <= tstart, and must join up
with(cgd, tapply(seq_along(id), id, function(i)
  all(tstart[i][-1] == tstop[i][-length(i)]))) |> table()

Verified with R 4.6.0 and survival 3.8.6

import numpy as np
import pandas as pd
from lifelines import CoxTimeVaryingFitter

# Where the long table comes from: the repeated labs in pbcseq, joined to
# the follow-up endpoint in pbc.
RD  = "https://vincentarelbundock.github.io/Rdatasets/csv/"
pbc = pd.read_csv(RD + "survival/pbc.csv")
seq = pd.read_csv(RD + "survival/pbcseq.csv")
base = (pbc.loc[pbc["id"] <= 312, ["id", "time", "status", "age"]]
           .rename(columns={"time": "futime"}))
measurements = (seq[["id", "day", "bili"]].merge(base, on="id")
                   .assign(log_bili=lambda d: np.log(d["bili"])))

# Long table columns: id / start / stop / event / covariates; many rows per person.
# lifelines has no tmerge(); cut the intervals yourself (pandas merge_asof helps).
long = (measurements
        .sort_values(["id", "day"])
        .assign(start=lambda d: d["day"],
                stop=lambda d: d.groupby("id")["day"].shift(-1)))
long["stop"] = long["stop"].fillna(long["futime"])
long["event"] = ((long["stop"] == long["futime"]) & (long["status"] == 2)).astype(int)
long = long[long["stop"] > long["start"]]          # empty intervals raise an error

ctv = CoxTimeVaryingFitter()
ctv.fit(long[["id", "start", "stop", "event", "log_bili", "age"]],
        id_col="id", event_col="event", start_col="start", stop_col="stop")
ctv.print_summary()

# Recurrent events: same API, but a robust variance is required
# (otherwise the standard errors are understated)
# Recurrent events: on the R side this is coxph(..., cluster = id).
# lifelines cannot do it -- CoxTimeVaryingFitter raises
# NotImplementedError for robust=True -- and skipping the robust variance
# understates the standard errors. Go back to R, or write the sandwich.

lifelines' CoxTimeVaryingFitter takes a long table with exactly the same column roles as R, but there is no equivalent of tmerge() — you assemble the intervals yourself. For recurrent events the robust variance is out of reach here: robust=True raises NotImplementedError, so that step goes back to R's coxph(..., cluster = id).

How to read this in a paper

When was the exposure determined. Whenever the grouping variable can only be known after baseline — received surgery, responded, completed the course, started a drug — go to the Methods and find out whether they used a time-dependent variable, a landmark, or nothing at all.
When the time-dependent covariate was measured, and how it entered the model. “The most recent value”, “the mean over the interval” and “the value lagged by a week” are three different definitions and give different results. None stated means not reproducible.
Whether recurrent events used a robust variance. Look for “robust”, “sandwich”, “cluster”, “covs(aggregate)”. Absent, widen the confidence intervals by roughly a fifth in your head.
Which recurrent event model. AG, PWP, WLW and frailty models rest on different assumptions; “a Cox model for recurrent events” is not enough.
How $L$ was chosen for a landmark analysis, and how many people it discarded. Was it in the protocol; and are the characteristics of those excluded — the people who died before $L$ — described.
Whether a time-dependent model was used for prediction. A sentence like “this model’s C statistic was X, better than the baseline model” attached to a model containing time-dependent covariates is an unfair comparison.

Common misuses

Misuse	Why it is wrong
Grouping at baseline by a status only known later	Immortal time bias; see B3-01
Confusing a time-dependent covariate with a time-varying coefficient	The first is a moving exposure, the second a violated proportional hazards assumption; the fixes are entirely different
Building the long table and never checking the intervals	Overlaps, negative lengths and events on the wrong row raise no error — they just quietly return the wrong answer
Letting a covariate take effect at the instant of the event	That writes outcome information into the exposure; a dying state gets read as a strong predictor
Andersen-Gill for recurrent events without a robust variance	Correlated repeated observations are counted as independent information, and standard errors are clearly understated
Keeping only the first event without saying so	The discarded events are usually the majority, and mostly come from the sickest patients
Choosing between AG, PWP and WLW without justification	Their assumptions and risk set definitions differ, and results can differ substantially
Choosing the landmark after seeing the data	The point estimate moves with L, so a post hoc choice is selective reporting
A landmark analysis that never says how many people were excluded	Readers cannot judge whether the remaining population still addresses the original question
Using a time-dependent model for survival curves or n-year survival	Those require future covariate trajectories, which is exactly what is unknown
Claiming better prediction from a time-dependent model’s C statistic	It uses information generated after the moment of prediction, so the comparison is unfair
Adjusting for an exposure-affected time-dependent variable as an ordinary covariate	It is both mediator and confounder; standard Cox cannot handle it, a marginal structural model can

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B3-06-time-varying.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

The two variance estimators for an Andersen-Gill model give the same hazard ratio but different standard errors. Which should be reported, and why?

Show the answer and why

Correct answer: 0.312 - it does not assume that one patient's repeated infections are independent

0.312 is the robust standard error. The naive 0.261 assumes a patient's repeated infections are independent, and they are not - some children simply relapse often, so their second event carries less new information than independence implies. Ignoring that makes intervals too narrow and p-values too small while the point estimate is unchanged, so nothing in the hazard ratio column looks wrong. 0.335 comes from the first-event-only model, which does avoid the correlation, at the cost of discarding more than half the events - a different question, not a fix for this one.

If the same CGD data are analysed using only each patient's first infection, what is given up?

Show the answer and why

Correct answer: 32 events are discarded, all of them somebody's second or later infection

76 events in total, 44 once only first events count, and the difference of 32 is entirely second-and-later infections - close to half the events, and not a random half: they fall precisely on the children who relapse most. So first-event-only analysis does not merely shrink the sample, it removes the tendency to relapse from the analysis altogether. Papers taking this route rarely say in the methods what was dropped; the reader sees only a smaller event count.

Landmark analysis excludes anyone whose event occurred before the landmark day. What happens when the landmark moves from day 14 to day 30?

Show the answer and why

Correct answer: 24 are excluded, and they are exactly the patients whose events came earliest

A landmark at day 30 excludes 24 patients; the earlier landmark excludes only 15 - later landmark, more exclusions. What matters is not the count but who: those excluded are the ones whose events came first, so the post-landmark population is healthier than the one enrolled. That is not a bias but the design working as intended, and what it buys is the removal of immortal time bias. The price is that conclusions apply to patients who survived to the landmark day and cannot be carried back to the original cohort. 79 is how many the day-30 version kept.

Chapters that use this method

Watch next

繁中Ming-Chieh Shih· 18 minIn Mandarin, continuing from part one into the Cox model and data structure. It has no section on time-dependent covariates as such, but it makes clear what shape of data the model consumes — a warm-up for the second section here.

繁中MeDA（臺大公衛洪弘教授）· 19 minA hands-on R session in Mandarin showing how survival data is prepared and modelled. Watch it first if you want the code on this page to actually run.

ENVery Normal· 15 minIn English, rebuilding survival analysis from the definition of the hazard function. The claim on this page that "the risk set is reassembled at every event time" lands much better after watching it.

Sources and licences

This page is original writing

Two things that get conflated

The (start, stop] long table: one person, many rows

tmerge(): turning a wide table into a long one

Repeated measurements: baseline value or current value

Recurrent events: Andersen-Gill and the robust variance

Landmark analysis: the blunter option that is easier to explain

An interpretive limit: a time-dependent hazard ratio cannot be used for prediction

Run it yourself

How to read this in a paper

Common misuses

Reproducing every number on this page

Read the figure

Chapters that use this method

Watch next

Sources and licences

`tmerge()`: turning a wide table into a long one