Writing

Survival analysis in R for molecular subgroups

November 8, 2022 · R, survival analysis, survminer, biostatistics

In brain tumor genomics, the question that follows every new classification is the same: do these groups do differently? IDH-mutant versus IDH-wildtype gliomas, 1p/19q codeleted versus intact, MGMT promoter methylated versus unmethylated, and now the DNA methylation classes. Defining the subgroup is usually the hard, interesting part. The survival analysis afterwards looks routine, and that is exactly where a lot of good molecular work picks up avoidable mistakes.

This post walks through the survival analysis part in R with the survival and survminer packages. I can’t share patient data, so everything runs on a dataset that ships with R. That dataset has no molecular measurements, and I’m not going to pretend otherwise. Instead I use a real binary tumor characteristic as a stand-in for a “subgroup”, and later a clearly simulated marker to demonstrate one particular pitfall. Swap in IDH status or a methylation class and the code is the same.

Kaplan-Meier curves for glioma-specific (top) and progression-free (bottom) survival of NF1-associated gliomas, split by molecular low-grade versus high-grade group, with numbers at risk underneath. Figure 3 from Lucas CHG, Sloan EA, Gupta R, et al. Multiplatform molecular analyses refine classification of gliomas arising in patients with neurofibromatosis type 1. Acta Neuropathologica 144, 747-765 (2022). Open access, CC BY.

Setup

install.packages(c("survival", "survminer"))
library(survival)
library(survminer)

survival is Terry Therneau’s package and comes with every R installation as a recommended package. survminer adds ggplot2-based plotting on top of it.

The data

survival::colon comes from an adjuvant chemotherapy trial in resected colon cancer, with three arms: observation, levamisole, and levamisole plus 5-FU. Run ?colon for the variable descriptions and the original trial references.

colon <- survival::colon
head(colon)
table(colon$etype)

The dataset has two rows per patient: etype == 1 is recurrence and etype == 2 is death. Mixing them up is an easy way to get nonsense, so the first step is to keep one endpoint. I use overall survival (death from any cause).

d <- subset(colon, etype == 2)

d$years    <- d$time / 365.25
d$age10    <- d$age / 10
d$grade    <- factor(d$differ, levels = 1:3,
                     labels = c("well", "moderate", "poor"))
d$subgroup <- factor(d$node4, levels = 0:1,
                     labels = c("<=4 nodes", ">4 nodes"))

str(d[, c("years", "status", "subgroup", "age10", "grade", "rx")])

node4 indicates more than four positive lymph nodes. It is not molecular, but it is structurally the same kind of variable as IDH status in a glioma cohort: a binary property of the tumor, measured at baseline, that we suspect is prognostic. I call it subgroup in the code so the analogy stays visible. The covariates are age (in decades, so the hazard ratio is per 10 years), tumor differentiation as a stand-in for grade, and treatment arm rx.

The survival object

Survival data has two parts: a time, and whether the event happened at that time. A patient who is still alive at their last follow-up is censored: we know they survived at least that long, but not how much longer. Throwing censored patients away, or treating their last follow-up as a death, both bias the answer. Survival methods use the partial information correctly.

head(Surv(d$years, d$status), 10)

Censored times print with a +. The event indicator is normally 0/1 (or FALSE/TRUE), but Surv() also accepts 1/2 coding, where 1 means censored and 2 means the event. survival::lung uses that coding, which surprises people the first time they see it. Check your own event coding with a table() before anything else.

All of this assumes censoring is non-informative: patients who drop out of follow-up are not systematically sicker or healthier than those who stay. No statistical test can fully check that for you, so it’s worth thinking through how follow-up actually ended in your cohort.

One number I report for every cohort is median follow-up, estimated with the reverse Kaplan-Meier method (flip the event indicator so censoring becomes the “event”):

survfit(Surv(years, 1 - status) ~ 1, data = d)

Kaplan-Meier curves

The Kaplan-Meier estimator gives the probability of surviving past each time point, accounting for censoring. With a grouping variable, you get one curve per group.

km <- survfit(Surv(years, status) ~ subgroup, data = d)
km
summary(km, times = c(1, 3, 5))

Printing km gives the number of patients, the number of events and the median survival for each group (with a confidence interval). summary(..., times =) gives survival probabilities at specific landmarks, which are often more useful than medians when one group never reaches 50%.

For the plot, ggsurvplot() does almost everything you’d want in one call:

ggsurvplot(km, data = d,
           conf.int = TRUE,
           pval = TRUE,
           risk.table = TRUE,
           surv.median.line = "hv",
           break.time.by = 1,
           xlab = "Years since study entry",
           legend.title = "",
           legend.labs = levels(d$subgroup))

Always show the risk table. The right-hand end of a Kaplan-Meier curve is often based on a handful of patients, and the risk table tells the reader how much to trust it. A curve that drops sharply after year six because two of the last five patients died is not a finding.

The log-rank test

The pval = TRUE above is a log-rank test. To run it directly:

survdiff(Surv(years, status) ~ subgroup, data = d)

The output compares the observed number of events in each group with the number expected if all groups shared the same survival curve, and turns that into a chi-square statistic. The log-rank test weights all time points equally and is most powerful when one group’s hazard is a constant multiple of the other’s. It tells you whether the curves differ, not by how much. survdiff(..., rho = 1) gives the Peto-Peto variant, which puts more weight on early differences.

With more than two groups, the default test is a global one: are any of the curves different? Pairwise comparisons need a multiple testing correction, which survminer handles:

survdiff(Surv(years, status) ~ rx, data = d)
pairwise_survdiff(Surv(years, status) ~ rx, data = d,
                  p.adjust.method = "BH")

The bigger limitation of a Kaplan-Meier comparison is that it is unadjusted. Molecular subgroups are rarely independent of clinical factors. IDH-mutant glioma patients are, on average, considerably younger than IDH-wildtype patients, and age is strongly prognostic on its own. A beautiful separation of Kaplan-Meier curves can be partly or mostly an age effect. To ask whether a subgroup carries information beyond the known factors, you need a regression model.

Cox proportional hazards models

The Cox model (Cox 1972) relates covariates to the hazard, the instantaneous risk of the event at time t among those still at risk:

h(t | x) = h0(t) * exp(b1*x1 + b2*x2 + ...)

The baseline hazard h0(t) is left completely unspecified, which is why the model is called semi-parametric. Each exp(b) is a hazard ratio: the multiplicative change in hazard for a one-unit change in that covariate, holding the others fixed. A hazard ratio of 2 means twice the instantaneous risk at every time point, not half the survival time.

Before fitting, I make a complete-case data frame. coxph() silently drops rows with missing values, and if models are fit on different subsets of patients, comparing them is meaningless.

dc <- na.omit(d[, c("years", "status", "subgroup", "age10", "grade", "rx")])
nrow(d); nrow(dc)

fit0 <- coxph(Surv(years, status) ~ age10 + grade + rx, data = dc)
fit1 <- coxph(Surv(years, status) ~ subgroup + age10 + grade + rx, data = dc)
summary(fit1)

summary() prints the coefficients, hazard ratios (exp(coef)), standard errors, Wald tests and 95% confidence intervals for each term, plus three global tests (likelihood ratio, Wald and score) for the model as a whole and a concordance statistic. The number of events is printed near the top. Look at it, because it is the effective sample size, not the number of patients.

To test whether the subgroup adds information beyond the clinical model, compare the nested models with a likelihood ratio test:

anova(fit0, fit1)

This is the honest version of “subgroup is an independent prognostic factor”. It only works because both models were fit to the same rows.

Hazard ratios and confidence intervals

For a table in a paper, I pull the hazard ratios and Wald confidence intervals directly:

hr <- exp(cbind(HR = coef(fit1), confint(fit1)))
round(hr, 2)

For factors, each hazard ratio is relative to the reference level, so grade gets one row for moderate versus well and one for poor versus well. Set the reference level deliberately with relevel() (in a glioma cohort I’d usually make the group with the best expected prognosis the reference, so the hazard ratios read as “how much worse”). A confidence interval that runs from 0.6 to 4.0 is telling you that the data are compatible with a modest protective effect and a large harmful one. Report the interval, not just whether it crosses 1.

Forest plots

ggforest(fit1, data = dc)

ggforest() draws every term with its hazard ratio, confidence interval, p-value and the number of patients per level, and prints the global statistics underneath. Pass the same data frame you used to fit the model, so the counts match.

Checking proportional hazards

The Cox model assumes each hazard ratio is constant over time. That is often roughly true and sometimes clearly false: a treatment effect that fades after the first couple of years, or a marker that matters early and not late. The standard check is based on scaled Schoenfeld residuals (Grambsch and Therneau 1994):

zph <- cox.zph(fit1)
zph
plot(zph, var = 1)   # first term in the model
ggcoxzph(zph)

The table gives a test for each term and a global test. The plots show the estimated coefficient as a function of time, with a smooth; a flat line is what proportional hazards looks like. Don’t rely on the p-values alone. In a large cohort, a trivial and harmless deviation can be “significant”, and in a small cohort a real problem can easily miss the cutoff. The plot shows you the size and shape of the deviation, which is what you need to decide whether it matters.

If a term clearly violates the assumption, there are a few standard options:

• If the variable is a nuisance you don’t need a hazard ratio for, stratify on it (next section).
• Split follow-up into periods and estimate a separate hazard ratio in each.
• Model the coefficient as an explicit function of time with tt() in coxph().

The period split is easy to explain to clinicians. survSplit() cuts each patient’s follow-up into pieces at the times you choose:

ds <- survSplit(Surv(years, status) ~ ., data = dc,
                cut = 2, episode = "period")
fit_pw <- coxph(Surv(tstart, years, status) ~ subgroup:strata(period) +
                  age10 + grade + rx, data = ds)
summary(fit_pw)

This gives one subgroup hazard ratio for the first two years and another for the time after, using the counting process form Surv(start, stop, event). Pick the cut points from clinical reasoning or from the residual plot, and say so in the methods.

Stratification

A stratified Cox model lets each stratum have its own baseline hazard, while the other coefficients are shared across strata:

fit_strat <- coxph(Surv(years, status) ~ subgroup + age10 + rx +
                     strata(grade), data = dc)
summary(fit_strat)

You no longer get a hazard ratio for the stratification variable, which is the price. In exchange, you make no proportional hazards assumption about it. I use stratification all the time for things like cohort or study site. If you combine a public dataset with an institutional one, their baseline survival can differ for reasons that have nothing to do with biology (different eras, different referral patterns), and strata(cohort) handles that without forcing a hazard ratio onto it. Note that ggforest() doesn’t draw strata terms, since there’s nothing to draw.

Stratification in the plotting sense, Kaplan-Meier curves within levels of another variable, is also worth doing:

ggsurvplot_facet(km, data = d, facet.by = "rx", pval = TRUE)

If the subgroup effect looks very different across treatment arms, that’s a question about interaction, which you can test directly:

fit_int <- coxph(Surv(years, status) ~ subgroup * rx + age10 + grade,
                 data = dc)
anova(fit1, fit_int)

This is the distinction between a prognostic marker (associated with outcome regardless of treatment) and a predictive one (associated with how much a treatment helps). MGMT promoter methylation and temozolomide in glioblastoma is the classic example of the second kind. Interaction tests need a lot of events to have any power, so a non-significant interaction in a small cohort is weak evidence of no interaction.

Pitfalls

These are the problems I see most often in molecular subgroup papers, roughly in order of how much damage they do.

Immortal time bias

Immortal time is follow-up during which, by construction, a patient could not have had the event. The classic case is comparing “patients who received treatment X” with “patients who didn’t” when X starts months after diagnosis: to be in the treated group, you had to survive long enough to get treated, so that group gets a free head start.

Genomics studies have their own version. Suppose tumors were profiled from tissue collected at a second surgery, or patients joined a sequencing study some time after diagnosis. If you then measure survival from diagnosis, everyone in the profiled cohort was guaranteed to survive until profiling, and the cohort looks better than it is. If profiling rates differ between subgroups, the comparison is biased too. The fixes are to measure time from a common origin at which everyone’s group is already known (a landmark analysis), to handle delayed entry with left truncation, or to treat a group that changes over time as a time-dependent covariate. Schematically, with your own column names:

# Schematic: substitute your own data and column names
# delayed entry: each patient is at risk from entry until exit
coxph(Surv(entry_time, exit_time, event) ~ subgroup + age, data = cohort)

tmerge() in the survival package is the tool for building time-dependent covariates, and the “timedep” vignette (vignette("timedep", package = "survival")) walks through it.

Small subgroups

Precision in survival analysis comes from events, not patients. A subgroup with forty patients and six deaths produces a hazard ratio with a confidence interval wide enough to be compatible with almost anything. Check the events per group before modeling:

with(dc, table(subgroup, status))
with(dc, table(grade, status))

A common rule of thumb is around ten events per estimated coefficient, and every factor level counts as a coefficient. Rare methylation classes hit this limit quickly. Collapsing rare levels with a clear biological rationale, or reporting them descriptively without a hazard ratio, is better than a forest plot full of intervals that span two orders of magnitude.

Dichotomizing at an optimized cutoff

A continuous marker (an expression score, a methylation beta value, a mutation count) gets turned into “high” versus “low” at whatever cutoff gives the smallest log-rank p-value. The resulting p-value is not valid, because it is the minimum over many tests that were run on the same data, and the effect size is inflated for the same reason (Altman et al. 1994).

This is easy to demonstrate with a marker that is pure noise. The marker below is simulated, has no relationship to survival by construction, and is attached to the real colon survival times:

min_logrank_p <- function(marker, time, status,
                          probs = seq(0.1, 0.9, by = 0.05)) {
  cuts <- quantile(marker, probs)
  p <- sapply(cuts, function(k) {
    high <- marker > k
    test <- survdiff(Surv(time, status) ~ high)
    pchisq(test$chisq, df = 1, lower.tail = FALSE)
  })
  min(p)
}

set.seed(2019)
p_min <- replicate(200, {
  noise <- rnorm(nrow(d))   # simulated: unrelated to outcome
  min_logrank_p(noise, d$years, d$status)
})
mean(p_min < 0.05)

If the cutoff search were harmless, about 5% of these noise markers would come out “significant”. You’ll get a fraction well above that. survminer::surv_cutpoint() is often used for this kind of search; it is based on maximally selected rank statistics, which can be corrected for the search, but the p-value people usually report afterwards is the naive log-rank p-value on the chosen split. Better options are to keep the marker continuous (with a spline via pspline() if the relationship isn’t linear), to use a cutoff fixed in advance from prior work, or to choose the cutoff in one cohort and test it in an independent one.

Multiple comparisons

Testing ten methylation classes against each other gives forty-five pairwise comparisons. Testing every gene for association with survival gives twenty thousand. Either way, some p-values will be small by chance. Use a global test first, adjust pairwise comparisons (pairwise_survdiff() with p.adjust.method = "BH", or p.adjust() on your own vector of p-values), and in methods sections say how many comparisons you made, including the ones that didn’t make it into the figures.

Competing risks

If the endpoint is something other than death from any cause, such as recurrence, then death without recurrence is a competing event: once it happens, recurrence can no longer be observed. Treating those deaths as censored and plotting one minus Kaplan-Meier overestimates the cumulative incidence of recurrence. The survival package handles this with a multi-state Surv() object, where the event is a factor whose first level means censored, and survfit() then returns Aalen-Johansen estimates:

rec <- subset(colon, etype == 1)
dth <- subset(colon, etype == 2)
cr <- merge(rec[, c("id", "time", "status", "node4")],
            dth[, c("id", "time", "status")],
            by = "id", suffixes = c(".rec", ".dth"))

cr$event <- with(cr, ifelse(status.rec == 1, 1, ifelse(status.dth == 1, 2, 0)))
cr$etime <- with(cr, ifelse(status.rec == 1, time.rec, time.dth)) / 365.25
cr$event <- factor(cr$event, levels = 0:2,
                   labels = c("censored", "recurrence", "death"))

cif <- survfit(Surv(etime, event) ~ node4, data = cr)
cif
plot(cif, xlab = "Years", ylab = "Cumulative incidence")

For regression, finegray() in survival and cmprsk::crr() fit Fine-Gray subdistribution hazard models, and cause-specific Cox models are another option. Which one answers your question depends on whether you care about etiology or prediction, and that is a longer post than this one.

Checklist before you send the figure

1. Pick one endpoint and one time origin, and make sure every patient’s subgroup was known at that origin.
2. Check event coding with table(); Surv() accepts 0/1 and 1/2, and they mean different things.
3. Report median follow-up (reverse Kaplan-Meier) and the number of events per group.
4. Show Kaplan-Meier curves with a risk table and confidence bands.
5. Use the log-rank test for the unadjusted comparison, with corrected pairwise tests if there are more than two groups.
6. Fit a Cox model with the clinical factors your audience would ask about (in gliomas: at least age, grade, extent of resection and treatment), on a complete-case data frame.
7. Test the subgroup’s added value with a likelihood ratio test between nested models.
8. Check proportional hazards with cox.zph() and look at the plots, not just the p-values.
9. Stratify on nuisance variables like cohort or site.
10. Don’t optimize cutoffs on the data you’re testing, and count every comparison you made.

For more depth, Therneau and Grambsch’s book is still the reference I reach for, and the vignettes that come with the survival package (browseVignettes("survival")) cover time-dependent covariates, multi-state models and competing risks with worked examples.

References

• Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association 53, 457-481 (1958).
• Cox DR. Regression models and life-tables. Journal of the Royal Statistical Society, Series B 34, 187-220 (1972).
• Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika 81, 515-526 (1994).
• Altman DG, Lausen B, Sauerbrei W, Schumacher M. Dangers of using “optimal” cutpoints in the evaluation of prognostic factors. Journal of the National Cancer Institute 86, 829-835 (1994).
• Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer, New York (2000).

← All writing