Time-Varying Covariates in a Cox Model in R (and Immortal-Time Bias)

Include a covariate that changes during follow-up — build counting-process (start, stop, event) data with tmerge, fit the model, and avoid the immortal-time bias that wrecks the naive analysis

Model covariates whose value changes during follow-up with a Cox model in R. Build the counting-process (start, stop, event) data layout with tmerge, fit coxph(Surv(start, stop, event) ~ …), interpret the time-varying hazard ratio, and see how the naive baseline-only analysis creates immortal-time bias on the Stanford heart-transplant data — plus when to prefer a simpler landmark analysis.

TipKey takeaways

• A time-varying covariate is a predictor whose value changes during follow-up — a treatment received later, a biomarker that drifts, a status that flips. A baseline Cox model can’t capture it.
• The fatal mistake is immortal-time bias: treating a future status (e.g. “ever transplanted”) as if it were known at baseline. A patient must survive long enough to reach that status, which manufactures a fake survival advantage. On the heart-transplant data the naive model reports an 82% risk reduction (HR ≈ 0.18) that vanishes (HR ≈ 1.0) once you do it right.
• Do it right with the counting-process layout: split each patient into (start, stop, event) intervals so the covariate takes its correct value during each interval. Build it with tmerge().
• Fit it with coxph(Surv(start, stop, event) ~ ...) — the formula and interpretation are otherwise identical to an ordinary Cox model.
• When the counting-process build is awkward, landmark analysis is the simpler robust alternative: pick a landmark time, classify everyone by their status at that time, and analyse only the survivors.
• A time-varying covariate (a changing value) is not the same as a time-varying coefficient (a fixed value whose effect changes — a PH violation).

Introduction

Does receiving a heart transplant improve survival? The obvious analysis — split patients into “transplanted” and “not transplanted”, then compare — gives a spectacular answer: transplant slashes the risk of death. It is also wrong. A patient can only be counted as “transplanted” if they lived long enough to receive the organ; anyone who died while waiting is forced into the “not transplanted” group. The transplant group is handed a stretch of guaranteed survival — immortal time — that has nothing to do with the operation. The comparison is rigged before it starts.

This is the central problem of time-varying covariates: predictors whose value changes during follow-up. Transplant status (0 then 1), a biomarker that rises, a treatment switched on partway through — none of these is a fixed baseline number, and pretending otherwise produces immortal-time bias, one of the most common and damaging errors in survival analysis.

The fix is to let the covariate change inside the model. You restructure each patient into short time intervals, give the covariate its correct value in each one, and fit a Cox model on that. This lesson shows the whole recipe in R with survival (the tmerge() builder and coxph()) and survminer (the forest plot), using the classic Stanford heart-transplant data so you can watch the biased result collapse into the honest one. You will learn to:

• recognise a time-varying covariate and the immortal-time bias that comes with mishandling it;
• build the counting-process (start, stop, event) data layout — by hand and with tmerge();
• fit and interpret coxph(Surv(start, stop, event) ~ ...);
• quantify the bias by comparing the naive and correct models;
• and choose landmark analysis when it is the better tool.

NoteThe two things people call “time-varying”

Don’t confuse them — they need different fixes:

• A time-varying covariate (this lesson): the value changes over time (transplant 0 → 1). You fix it with the counting-process layout.
• A time-varying coefficient: the value is fixed, but its effect (the hazard ratio) changes over time — a violation of the proportional-hazards assumption. You fix it with tt(), a time interaction, or stratification.

The trap: immortal-time bias, made visible

Start with the wrong analysis, because seeing it fail is the fastest way to understand the fix. The jasa dataset (Stanford heart-transplant study) ships with the survival package with one row per patient: a follow-up time (futime), a death indicator (fustat), and a transplant flag marking whether the patient was ever transplanted.

Fit a Cox model treating “ever transplanted” as a fixed baseline covariate — exactly the analysis instinct tells you to run:

library(survival)

# WRONG: `transplant` here = "ever transplanted", a FUTURE status used as a baseline value
bad.cox <- coxph(Surv(futime, fustat) ~ transplant + age + surgery, data = jasa)
summary(bad.cox)$conf.int

           exp(coef) exp(-coef) lower .95 upper .95
transplant 0.1795839   5.568429 0.1040357 0.3099933
age        1.0606546   0.942814 1.0298264 1.0924057
surgery    0.6576938   1.520464 0.3177452 1.3613459

Read the transplant row: HR ≈ 0.18, a 95% CI of roughly 0.10–0.31, p < 0.0001. Taken at face value, transplant cuts the hazard of death by 82% — a miracle drug. It is an artefact. Every patient in the transplanted group survived until their operation; that waiting time is immortal, and the model credits it to the transplant.

WarningWhy this is biased, in one sentence

To be labelled “transplanted”, a patient had to stay alive until the transplant — so the label itself selects for survival, and the hazard ratio measures that selection, not the treatment.

The fix: the counting-process layout

The honest model lets transplant be 0 before the operation and 1 after it — for the same patient. You do that by splitting each patient into time intervals and giving the covariate its correct value in each. This is the counting-process (or (start, stop, event)) data layout.

The heart dataset is the same Stanford study already in this layout. Look at three patients:

library(survival)

# Patients 2 (never transplanted), 3 and 4 (transplanted partway through)
heart[heart$id %in% c(2, 3, 4), c("id", "start", "stop", "event", "transplant", "age", "surgery")]

  id start stop event transplant       age surgery
2  2     0    6     1          0  3.835729       0
3  3     0    1     0          0  6.297057       0
4  3     1   16     1          1  6.297057       0
5  4     0   36     0          0 -7.737166       0
6  4    36   39     1          1 -7.737166       0

Read the rows:

• Patient 2 has one interval (0, 6], transplant = 0 throughout — never transplanted, died at day 6.
• Patient 4 has two intervals: (0, 36] with transplant = 0 (still waiting), then (36, 39] with transplant = 1 (after the operation), with the event on the second. The covariate changes value at day 36 — exactly what a single baseline row cannot express.

Each interval is a separate row, but event is 1 only on the interval where the event actually happened (every earlier interval is event = 0, a censored fragment that still contributes its person-time). That is the whole trick: the patient is “at risk with transplant = 0” for the first stretch and “at risk with transplant = 1” for the second.

Build it yourself with tmerge()

You rarely receive data pre-split. The survival package’s tmerge() builds the counting-process layout from a baseline table plus the times at which things change. Start from a per-patient baseline (jasa), compute each patient’s follow-up time and transplant time, then merge:

library(survival)

# One row per patient: id, follow-up time, transplant time (the day status flips), death
jasa$subject <- seq_len(nrow(jasa))
tdata <- with(jasa, data.frame(
  subject = subject,
  futime  = pmax(0.5, fu.date - accept.dt),                                 # follow-up days
  txtime  = ifelse(is.na(tx.date), (fu.date - accept.dt) - 1,               # never transplanted
                                   (tx.date - accept.dt)),                  # day of transplant
  fustat  = fustat
))

# tmerge: lay down the death event, then add `trt` as a time-dependent covariate (tdc)
sdata <- tmerge(jasa, tdata, id = subject,
                death = event(futime, fustat),
                trt   = tdc(txtime))

# the resulting counting-process rows
head(sdata[, c("subject", "tstart", "tstop", "death", "trt")], 8)

  subject tstart tstop death trt
1       1      0    48     0   0
2       1     48    49     1   1
3       2      0     4     0   0
4       2      4     5     1   1
5       3      0    15     1   1
6       4      0    35     0   0
7       4     35    38     1   1
8       5      0    16     0   0

Two helpers do the work: event(time, status) marks where the event lands, and tdc(time) — time-dependent covariate — flips trt from 0 to 1 at txtime. tmerge() returns one row per interval, with tstart/tstop boundaries and the covariate carrying its correct value in each. A patient transplanted on day 36 gets a trt = 0 interval up to day 36 and a trt = 1 interval after — the same structure you saw in heart, built from scratch.

TipReading the tmerge result

tmerge() prints a tcount attribute summarising how each addition landed (how many changes fell within follow-up, before it, etc.). Check it: changes that land after a patient’s follow-up, or stack on the same day, are usually a sign the input times are off.

Fit the time-varying Cox model

With the data in (start, stop, event) form, fitting is the same coxph() call as an ordinary Cox model — only the Surv() response gains a start time. Use heart (it carries age and surgery too):

library(survival)

# RIGHT: transplant is 0 before the operation, 1 after — within the same patient
tv.cox <- coxph(Surv(start, stop, event) ~ transplant + age + surgery, data = heart)
summary(tv.cox)

Call:
coxph(formula = Surv(start, stop, event) ~ transplant + age + 
    surgery, data = heart)

  n= 172, number of events= 75 

                coef exp(coef) se(coef)      z Pr(>|z|)  
transplant1  0.01610   1.01623  0.30859  0.052   0.9584  
age          0.03054   1.03101  0.01389  2.198   0.0279 *
surgery     -0.77333   0.46147  0.35967 -2.150   0.0315 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

            exp(coef) exp(-coef) lower .95 upper .95
transplant1    1.0162     0.9840     0.555    1.8606
age            1.0310     0.9699     1.003    1.0595
surgery        0.4615     2.1670     0.228    0.9339

Concordance= 0.6  (se = 0.036 )
Likelihood ratio test= 10.72  on 3 df,   p=0.01
Wald test            = 9.68  on 3 df,   p=0.02
Score (logrank) test = 10  on 3 df,   p=0.02

Interpret the output:

• transplant: HR ≈ 1.02, 95% CI 0.55–1.86, p = 0.96. Once transplant status is handled correctly, there is no detectable survival benefit — the HR sits on 1 and the CI is wide on both sides. The miraculous 82% reduction was entirely immortal-time bias.
• age: HR ≈ 1.03 per year, p = 0.03. Older patients have a modestly higher hazard (age here is centred at 48, so the HR is per extra year).
• surgery: HR ≈ 0.46, p = 0.03. Prior bypass surgery is associated with a lower hazard.
• n = 172 rows for 103 patients, 75 events — the rows are intervals, not people; coxph() correctly treats the multiple intervals from one patient as the same subject’s person-time.

The bias, quantified

Put the two transplant hazard ratios side by side — this contrast is the lesson:

library(survival)
library(survminer)

bad.cox <- coxph(Surv(futime, fustat) ~ transplant + age + surgery, data = jasa)
tv.cox  <- coxph(Surv(start, stop, event) ~ transplant + age + surgery, data = heart)

# Pull the transplant HR + 95% CI from each model into a small data frame
hr_row <- function(model, term, label) {
  ci <- summary(model)$conf.int[term, c("exp(coef)", "lower .95", "upper .95")]
  data.frame(model = label, hr = ci[1], lo = ci[2], hi = ci[3])
}
df <- rbind(
  hr_row(bad.cox, "transplant",  "Naive baseline-only\n(immortal-time bias)"),
  hr_row(tv.cox,  "transplant1", "Correct time-varying")
)

ggplot(df, aes(x = hr, y = model)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
  geom_errorbar(aes(xmin = lo, xmax = hi), width = 0.15, orientation = "y", colour = "#3a86d4") +
  geom_point(size = 3, colour = "#3a86d4") +
  scale_x_log10(breaks = c(0.1, 0.2, 0.5, 1, 2)) +
  labs(x = "Hazard ratio for transplant (log scale)", y = NULL) +
  theme_minimal(base_size = 12)

The naive point sits far left of 1 (HR ≈ 0.18, a strong “protective” effect); the correct time-varying point sits on the line (HR ≈ 1.0, no effect). Same data, same patients — the only difference is whether transplant was allowed to change value at the right moment. That gap is the bias.

Report

Transplant status was modelled as a time-varying covariate using the counting-process formulation (coxph(Surv(start, stop, event) ~ ...)) on the Stanford heart-transplant data (103 patients, 75 deaths). After correctly accounting for the time of transplantation, transplant was not associated with survival (HR = 1.02, 95% CI 0.55–1.86, p = 0.96). A naive analysis treating transplantation as a fixed baseline characteristic spuriously suggested a large benefit (HR = 0.18, 95% CI 0.10–0.31), an artefact of immortal-time bias.

Landmark analysis — the simpler robust alternative

Building counting-process data is fiddly, and when the question is really “given a patient’s status at a fixed point, what happens next?”, landmark analysis is cleaner and harder to get wrong. The recipe:

1. Pick a landmark time (e.g. day 30) — chosen for clinical/biological reasons, not from the data.
2. Drop everyone who died or was censored before it — the analysis starts the clock at the landmark.
3. Classify the survivors by their covariate value as of the landmark.
4. Fit an ordinary Cox model from the landmark onward.

Because every subject was alive at the landmark and is classified by their then-current status, there is no immortal time to bias the comparison. Here it is on heart, landmarking at day 30:

library(survival)

# Build a per-patient table: status at day 30, and survival measured FROM day 30
landmark <- 30

per_patient <- do.call(rbind, lapply(split(heart, heart$id), function(d) {
  d <- d[order(d$stop), ]
  total_stop  <- max(d$stop)                                   # last follow-up time
  died        <- d$event[which.max(d$stop)]                    # 1 if died at the end
  tx_by_30    <- as.integer(any(d$transplant == 1 & d$start < landmark))
  data.frame(id = d$id[1], stop = total_stop, event = died,
             tx30 = tx_by_30, age = d$age[1], surgery = d$surgery[1])
}))

# Keep only patients alive past the landmark; survival time is measured FROM the landmark
lm_data <- per_patient[per_patient$stop > landmark, ]
lm_data$time_from_lm <- lm_data$stop - landmark

lm.cox <- coxph(Surv(time_from_lm, event) ~ tx30 + age + surgery, data = lm_data)
summary(lm.cox)$conf.int

        exp(coef) exp(-coef) lower .95 upper .95
tx30    0.9587380  1.0430378 0.5497392 1.6720267
age     1.0378955  0.9634881 1.0036068 1.0733557
surgery 0.4405816  2.2697269 0.1961135 0.9897949

The landmark model asks a well-posed question — among patients alive at day 30, does having been transplanted by then predict later survival? — and, like the time-varying model, finds no clear benefit for tx30. The trade-off: landmark analysis discards early events and the information before the landmark, and the answer can depend on the landmark you pick (report a sensitivity analysis at a couple of times). Use it when there is a natural decision point; use the counting-process model when you want every event and a covariate that updates continuously.

TipWhich one should I use?

• Counting-process / tmerge → the covariate updates repeatedly (a biomarker measured at many visits), you want maximum efficiency (every event counts), and you can build the intervals reliably.
• Landmark analysis → there is a natural decision point, you want a simple, hard-to-bias answer, and you can afford to drop pre-landmark events.

A note on time-varying coefficients

A time-varying covariate is a changing value. A different problem — a fixed covariate whose effect changes over time — is a time-varying coefficient, i.e. a violation of the proportional-hazards assumption. You diagnose it with cox.zph() and fix it with tt(), a time × covariate interaction, or stratification. That is a separate workflow covered in Testing the proportional hazards assumption. Keep the two straight: value changes → counting process (this lesson); effect changes → tt() / non-PH remedies.

Try it live

Build the time-varying model and compare it with the biased one. Try a different landmark time, or add year (time since the study began) to the model. The sandbox boots on first Run.

🟢 With an AI agent

Ask Prova “my covariate changes during follow-up — help me build counting-process (start, stop, event) data with tmerge and fit a time-varying Cox model, and check I’m not creating immortal-time bias” — it answers with survival code you can run on your own data. The runtime is the judge. Ask Prova →

Common issues

You forgot Surv(start, stop, event) and used Surv(time, event). With counting-process data the response needs three arguments — a start, a stop, and the event. Passing only two collapses every interval to “time from zero”, silently double-counts person-time, and gives wrong estimates. The giveaway: your n equals the number of intervals, not the number you intended, and standard errors look too small.

The event is coded on the wrong interval. Only the interval where the event actually happened may carry event = 1; every earlier interval for that patient is event = 0. If you copy the patient’s final status onto all their rows, you record several “events” for one death. Let tmerge’s event() place it for you, or set it explicitly on the last interval only.

Overlapping or zero-length intervals after tmerge. Intervals must be non-overlapping and have stop > start. Two changes on the same day, or a transplant time equal to the death time, produce a zero-length interval that coxph() drops or errors on. Check the tcount attribute, nudge tied times apart by a fraction of a day, and make sure each tstop exceeds its tstart.

Frequently asked questions

NoteWhat is immortal time bias in survival analysis?

Immortal-time bias arises when a stretch of follow-up during which the event cannot occur is misassigned to a treatment group. The classic case: classifying patients by a status they can only reach by surviving (e.g. “received a transplant”, “responded to treatment at 12 weeks”). Those patients had to stay alive to be labelled, so the label selects for survival and manufactures a fake benefit. Fix it by modelling the status as a time-varying covariate (counting-process layout) or with landmark analysis — never by treating a future status as a baseline value.

NoteHow do I include a time-varying covariate in a Cox model in R?

Restructure the data into the counting-process layout — one row per (start, stop, event) interval, with the covariate taking its correct value in each interval — then fit coxph(Surv(start, stop, event) ~ covariates, data). Build the intervals with tmerge() from the survival package: event() places the outcome, and tdc() adds a time-dependent covariate that flips value at a given time. The model formula and the hazard-ratio interpretation are otherwise identical to an ordinary Cox model.

NoteWhat is the counting-process data format?

It is a long layout where each subject is split into one or more time intervals, each a row with a start time, a stop time, and an event indicator (1 only on the interval containing the event). Covariates can take a different value on each interval, so a status that changes mid-follow-up (transplant 0 → 1) is represented exactly. You fit it with Surv(start, stop, event). The survival package’s heart dataset is a worked example; tmerge() builds the layout from baseline data plus change times.

NoteWhat is the difference between a time-varying covariate and a time-varying coefficient?

A time-varying covariate has a value that changes over time (a biomarker, a treatment switched on later) — you handle it with the counting-process layout and tmerge(). A time-varying coefficient is a fixed covariate whose effect (hazard ratio) changes over time — a violation of the proportional-hazards assumption — which you handle with tt(), a time interaction, or stratification (see testing the PH assumption). Different problems, different fixes; don’t conflate them.

NoteWhen should I use landmark analysis instead of a time-varying Cox model?

Use landmark analysis when there is a natural decision point and you want a simple, hard-to-bias answer: pick a landmark time, keep only patients alive at it, classify them by their status as of that time, and fit an ordinary Cox model from the landmark onward. It avoids immortal time by construction but discards pre-landmark events and depends on the landmark chosen. Prefer the counting-process / time-varying model when the covariate updates repeatedly or you want every event to contribute.

Test your understanding

ImportantPractice

1. Run it. In the live cell below, fit the time-varying Cox model on heart, but leave out surgery (keep transplant and age). Fill the blank. Does the transplant hazard ratio stay near 1?
2. Interpret. A colleague reports that “patients who switched to the new drug during the trial lived much longer (HR = 0.4)” — and they classified the switch from study entry. Name the bias and the fix.

NoteHint

Fill the blank with transplant. The response is always Surv(start, stop, event) for counting-process data. Look at the transplant1 row’s exp(coef) and its p-value — a HR near 1 with a non-significant p-value means no detectable effect.

NoteSolution

tv.cox <- coxph(Surv(start, stop, event) ~ transplant + age, data = heart)
summary(tv.cox)
#> transplant1 HR is still ~1 (p large) — no transplant benefit once status is time-varying.
#> Dropping `surgery` barely moves it: the transplant null is robust, not an artefact of adjustment.

For question 2: this is immortal-time bias — patients had to survive long enough to switch drugs, so the “switcher” group is selected for survival, faking the HR = 0.4. The fix is to model the switch as a time-varying covariate (counting-process layout with tmerge()), so each patient contributes drug = old person-time before the switch and drug = new after — or to use landmark analysis at a fixed time.

TipQuick check

A study classifies patients as “responders” vs “non-responders” based on response measured at 8 weeks, then compares survival from study entry. Why is the responders’ survival advantage likely overstated?

NoteShow answer

Immortal-time bias. To be a “responder”, a patient had to survive to the 8-week assessment — anyone who died earlier is automatically a non-responder. The responder group is therefore guaranteed at least 8 event-free weeks that have nothing to do with responding, inflating its apparent survival. Correct it by treating response as a time-varying covariate (it becomes 1 only at week 8) or with a landmark analysis at week 8 using only the survivors.

Conclusion

You can now handle a covariate that changes during follow-up without falling into the trap that catches so many survival analyses. The key insight is that a future status — transplant, response, a drug switch — cannot be a baseline value: pretending it is creates immortal-time bias and a hazard ratio that is pure artefact (the heart-transplant data’s “82% benefit” that evaporated to nothing). The honest analysis puts the covariate in the counting-process (start, stop, event) layout — built by hand or, in practice, with tmerge() — and fits the same coxph() you already know. When that build is awkward or the question has a natural decision point, landmark analysis gives a simpler, bias-resistant answer. Keep the distinction from a time-varying coefficient (a PH violation) clear, and you have the full toolkit for survival data that refuses to stand still.

Related lessons

• Cox proportional hazards model — the base model this extends; fit and interpret hazard ratios with 95% CIs and the ggforest plot. · Testing the proportional hazards assumption — the other “time-varying” problem (a coefficient whose effect changes), with cox.zph and the tt() remedy. · Stratified Cox models — give a nuisance covariate its own baseline hazard, an alternative when PH fails. · Cox model building and variable selection — choose which covariates (time-varying or not) belong in the model.
• Go further: Survival Analysis — the full series. · Biostatistics — the pillar.

Was this page helpful?

Prove you can do it. Master the whole Survival Analysis in R series — track your path, build projects, and earn a certificate.

Start free →

Go Pro — unlimited Prova on your own data and a verifiable certificate that proves the skill.

from $15/mo billed yearly

Go Pro →

✓ You're Pro — keep going. The runtime is the judge.

Ready to level up?

Try this on your own data with Prova — the runtime is the judge.Open Prova →

Get new R & Python lessons by email

Practical, reproducible, no spam. Unsubscribe anytime.

Double opt-in. We never share your email.

Share this pageXLinkedInRedditHN

Note

Every result on this page was produced by the code shown, run at build time against a pinned R environment — edit any block and Run to reproduce it yourself.