Time-Varying Covariates in Survival Models

The Cox model becomes:

If \(e^\beta=0.70\), then at any time \(t\), patients currently receiving treatment have an estimated hazard 30% lower than comparable patients currently not receiving treatment, according to this model.

That interpretation is an association conditional on the model, not automatically a causal treatment effect.

Why "Ever Treated" Is Usually Wrong

Consider this incorrect coding:

treated_ever <- ifelse(patient_ever_received_treatment, 1, 0)

Suppose a patient receives treatment on day 100.

If the patient is classified as treated from day 0, the analysis incorrectly attributes the first 100 days of follow-up to the treated group.

The patient could not have benefited from a treatment that had not yet been administered.

This can produce substantial bias.

Immortal Time Bias

The problem becomes particularly important when exposure status is determined by an event that occurs after baseline.

Suppose researchers classify patients according to whether they eventually receive a treatment.

A patient who eventually receives treatment must survive long enough to receive it.

That period before treatment initiation can therefore be an immortal period for the eventually-treated group: the patient had to remain alive and event-free long enough to become classified as treated.

The Correct Time-Varying Solution

Instead of:

from baseline, define treatment according to the actual initiation time.

and:

This preserves the temporal ordering of exposure and outcome.

A Complete Worked Example

Consider a hypothetical Phase III oncology study evaluating time to progression.

Patients are followed from randomization until progression, death, or censoring.

Some patients receive a dose modification during follow-up.

We want to investigate whether the current treatment dose is associated with progression hazard.

Suppose Patient 25 has the following history:

The patient begins at dose level 1, increases to dose level 2 on day 28, and returns to dose level 1 on day 56.

Progression occurs at day 83.

The patient therefore contributes three intervals.

Model Specification

Let \(D_i(t)\) represent the current dose level.

A simple model is:

The hazard ratio for a one-unit increase in dose is:

For example, if:

then:

Thus, a one-unit higher current dose level is associated with approximately a 28% higher instantaneous hazard, assuming the model is otherwise correctly specified.

Time-Varying Biomarkers

Another common application is repeated laboratory measurement.

Suppose C-reactive protein is measured at baseline and every four weeks.

A model might use:

If:

then each one-unit increase in the current CRP value is associated with a 15% higher hazard, conditional on the other variables in the model.

The scale of the biomarker matters enormously.

For a biomarker measured over a wide range, it may be more interpretable to use a clinically meaningful increment.

For example:

would represent the hazard ratio associated with a 10-unit increase.

Categorical Time-Varying Covariates

A time-varying covariate does not have to be continuous.

For example, disease status could be categorized as:

The model could then estimate:

with:

The interpretation compares the hazard while a patient is in the specified disease state with the reference state.

Multiple Time-Varying Covariates

A model can contain several changing predictors simultaneously.

For example:

• Current treatment exposure
• Current dose
• Current biomarker
• Current disease status

can all be represented in the same model.

Start-Stop Intervals

A common data representation is:

For example:

Each patient can have multiple rows.

The event indicator is typically 1 only on the final interval containing the event.

R: Using the survival Package

In R, time-varying Cox models can be fit using the coxph() function from the survival package.

library(survival)

fit <- coxph(
  Surv(start, stop, event) ~ treatment,
  data = tv_data
)

summary(fit)

The key difference from a standard Cox model is the survival response:

Surv(start, stop, event)

rather than:

Surv(time, event)

The first form tells R that the data are represented using start-stop intervals.

R Worked Example

Consider:

tv_data <- data.frame(
  id = c(1,1,1,2,2,3,3),
  start = c(0,30,60,0,40,0,50),
  stop = c(30,60,90,40,100,50,120),
  treatment = c(0,1,1,0,1,0,1),
  event = c(0,0,1,0,1,0,1)
)

The Cox model is:

fit <- coxph(
  Surv(start, stop, event) ~ treatment,
  data = tv_data
)

summary(fit)

The treatment coefficient can be exponentiated to obtain the hazard ratio:

exp(coef(fit))

and a confidence interval can be obtained with:

exp(confint(fit))

Adding Baseline Covariates

Time-varying covariates can be combined with ordinary baseline covariates.

fit <- coxph(
  Surv(start, stop, event) ~
    treatment +
    age +
    sex +
    stage,
  data = tv_data
)

Here, treatment varies over time while age, sex, and baseline stage remain fixed.

The model therefore combines:

Using the tt() Function

R's coxph() function also supports transformations in which the effect of a baseline covariate changes as a function of time.

For example:

fit <- coxph(
  Surv(time, event) ~ age + tt(age),
  data = data,
  tt = function(x, t, ...) x * log(t)
)

This is a different use of "time-dependent" modeling.

The patient's age itself does not change during the study.

Instead, the coefficient associated with age is allowed to have a time-dependent component.

SAS: PROC PHREG

SAS PROC PHREG can analyze counting-process data using the start and stop times.

proc phreg data=tv_data;
    model (start, stop)*event(0) =
        treatment
        age
        sex
        stage;
run;

The notation:

(start, stop)*event(0)

defines the interval-censored risk representation used for the counting-process formulation.

SAS Example with a Time-Varying Treatment

proc phreg data=tv_data;
    model (start, stop)*event(0) = treatment age sex stage;
run;

The treatment variable can take different values on different rows for the same patient.

For example:

The model therefore uses treatment status appropriate to each interval.

Creating Start-Stop Data

A major practical challenge is constructing the analysis dataset correctly.

Suppose a patient has:

• Enrollment at day 0
• Treatment initiation at day 35
• Biomarker measurement at day 60
• Treatment discontinuation at day 90
• Event at day 120

The analysis dataset may need intervals such as:

Every time a covariate changes, a new interval begins.

The General Data-Construction Algorithm

Data Quality Checks

Time-varying survival analyses are particularly sensitive to data-structure errors.

At minimum, verify:

• Intervals are ordered chronologically.
• There are no unintended overlapping intervals.
• There are no unexplained gaps in the patient's risk period.
• The final interval ends at the event or censoring time.
• Only the event-containing interval has an event indicator of 1.
• Covariate changes occur at the correct times.
• No covariate value uses information observed after the interval.
• Patients are not duplicated accidentally.

A Useful Validation Check in R

tv_data[order(tv_data$id, tv_data$start), ]

For each patient, inspect whether:

when contiguous intervals are expected.

You can also check for multiple events:

aggregate(
  event ~ id,
  data = tv_data,
  sum
)

For ordinary time-to-first-event analysis, the total event count per patient should generally be at most one.

Time-Varying Covariates and Causality

This is one of the most important advanced issues.

A time-varying covariate may itself be affected by prior treatment or evolving disease status.

For example:

The time-varying biomarker is therefore both affected by treatment and predictive of future treatment and outcome.

A conventional time-varying Cox regression may estimate an association, but it does not automatically recover a causal treatment effect in the presence of such time-dependent confounding.

Time-Dependent Confounding

Suppose disease severity \(L(t)\) changes over time.

Treatment \(A(t)\) is influenced by severity:

But treatment also influences future severity:

and severity influences future event risk:

This creates feedback between treatment and the time-varying confounder.

A standard regression adjustment can be problematic because the confounder may simultaneously be:

• A predictor of future treatment
• A predictor of the outcome
• A consequence of previous treatment

Landmark Analysis

A landmark analysis is another approach for handling certain time-dependent exposure questions.

Choose a landmark time \(t_L\).

Only patients alive and event-free at \(t_L\) are included.

Exposure status at the landmark is then used to define groups for subsequent survival.

The analysis then asks a question such as: Among patients alive and event-free at day 90, is survival after day 90 associated with treatment status at day 90?

Time-Varying Covariate vs. Landmark Analysis

Time-Varying Covariates and the Proportional Hazards Assumption

It is important to understand that adding a time-varying covariate does not automatically solve proportional hazards violations.

Consider:

The coefficient \(\beta\) can still be assumed constant.

If instead the treatment effect itself changes with time, the model may need:

or a model containing an interaction between treatment and a function of time.

These are separate modeling decisions.

A Combined Model

In some situations both the covariate and its effect can vary over time:

This is considerably more complex and should be motivated by the scientific question and supported by sufficient data.

Time-Dependent Covariates Based on Events

A covariate can change when an intermediate event occurs.

For example, define:

Before hospitalization:

After hospitalization:

This allows the analysis to compare the hazard before and after hospitalization.

But again, the interpretation is associational.

Patients who become hospitalized may already be becoming sicker, and hospitalization may therefore be a marker of disease severity rather than a randomized exposure.

Internal vs. External Time-Varying Covariates

A useful conceptual distinction is between external and internal time-varying covariates.

External Covariate

An external covariate evolves independently of the patient's event process, at least conceptually.

Examples might include:

• Calendar-time environmental exposure
• Ambient temperature
• External treatment availability

Internal Covariate

An internal covariate is generated by the patient's own biological or clinical process.

Examples include:

• Tumor burden
• Biomarker concentration
• Current disease status
• Physiological measurements

Internal covariates require particular care because their evolution can be closely connected to the outcome process.

Last Observation Carried Forward: A Warning

Suppose a biomarker is measured every four weeks.

If the measurement at week 8 is 5.2, the analyst might use 5.2 until the week 12 measurement becomes available.

That can be a reasonable working definition of the current observed biomarker.

However, it should not be confused with a general solution to missing data.

Measurement Timing Matters

Suppose an event occurs on day 45 and the biomarker is measured on day 60.

The day-60 biomarker cannot be used as the covariate value for the day-45 risk set.

Doing so introduces future information into the model.

A safe temporal ordering is:

Lagged Time-Varying Covariates

In some applications, the biologically relevant question concerns the prior measurement rather than the instantaneous current measurement.

For example:

where \(L\) is a prespecified lag.

For example, a 7-day lag might be used when the scientific hypothesis is that the biomarker measured one week earlier predicts subsequent risk.

Lagging can also help ensure that the covariate clearly precedes the outcome.

Why the Choice of Lag Should Be Prespecified

Trying many lag periods and selecting the one producing the strongest association can inflate false-positive findings and produce unstable results.

A clinically justified lag should ideally be selected before examining the final outcome associations.

Time-Varying Treatment and Clinical Trials

In randomized trials, treatment assignment at baseline is ordinarily fixed.

Therefore, the primary intention-to-treat survival analysis usually uses treatment assignment as a baseline covariate.

A time-varying treatment variable answers a different question.

For example:

• What is the association between current treatment exposure and risk?
• What happens after treatment discontinuation?
• Is cumulative exposure associated with outcome?
• Is dose intensity associated with risk?

These analyses can be useful as secondary or exploratory analyses but should not automatically replace the prespecified primary analysis.

Cumulative Exposure

Sometimes the relevant variable is cumulative exposure rather than current exposure.

Define:

where \(A_i(u)\) represents treatment exposure.

If treatment is binary, \(C_i(t)\) can represent cumulative time exposed.

For example, a patient treated continuously for 60 days has:

while a patient treated for only 20 days has:

The survival model can then include cumulative exposure as a time-varying covariate.

Current Exposure vs. Cumulative Exposure

These are different estimands and should not be treated as interchangeable.

Time Since Exposure Initiation

Another useful time-varying variable is time since treatment initiation.

Define:

library(survival)
library(splines)

fit <- coxph(
  Surv(start, stop, event) ~
    ns(biomarker, df = 3),
  data = tv_data
)

library(survival)

# Time-varying Cox model
fit <- coxph(
  Surv(start, stop, event) ~
    treatment +
    biomarker +
    age +
    sex +
    stage,
  data = tv_data
)

# Model summary
summary(fit)

# Hazard ratios
exp(coef(fit))

# 95% confidence intervals
exp(confint(fit))

proc phreg data=tv_data;

    model (start, stop)*event(0) =
        treatment
        biomarker
        age
        sex
        stage;

run;

tv_data[
  order(tv_data$id, tv_data$start),
  c("id","start","stop","treatment","biomarker","event")
]