ggsurvplot() returns
a compound object – a list holding the survival curve
($plot), the risk table ($table), the data
actually plotted ($data.survplot), and more – and
surv_summary() returns a plain data frame with one row per
event/censoring time (time, surv,
n.risk, n.event, n.censor,
strata, …). Because both are ordinary ggplot2 / data-frame
objects, a lot of “can survminer do X?” questions are answered by
composing on what these already give you rather than by a
dedicated argument.
This page collects such recipes. Each is self-contained and uses the
lung data set from the survival package.
Some reporting guidelines ask that a survival curve be drawn with a dashed line beyond the point where few patients remain at risk – for example the guidelines for reporting mortality and morbidity after cardiac valve interventions recommend a dashed line once the number at risk falls below 10% of the initial cohort (survminer issue #559).
surv_summary() carries n.risk per stratum,
so you can find, for each curve, the first time it drops below the
threshold, split the curve into a solid “dense” part and a dashed
“sparse” tail, and draw the two with different linetypes.
Repeating the boundary point in the solid part makes the two step
segments meet with no gap.
# KM curves that switch solid -> dashed once the number at risk drops below
# `frac` of each group's initial cohort (default 10%).
ggsurv_dashed_below <- function(fit, data, frac = 0.10, ...) {
ss <- surv_summary(fit, data = data) %>%
group_by(strata) %>%
mutate(onset = { s <- n.risk < frac * max(n.risk)
if (any(s)) min(time[s]) else Inf },
phase = ifelse(time >= onset, "sparse", "dense")) %>%
ungroup()
# repeat the first sparse point in the solid phase so the two steps join
bridge <- ss %>% filter(phase == "sparse") %>%
group_by(strata) %>% slice_min(time, n = 1) %>% ungroup() %>%
mutate(phase = "dense")
ggplot(bind_rows(ss, bridge), aes(time, surv, colour = strata)) +
geom_step(aes(linetype = phase), linewidth = 0.9) +
geom_point(data = filter(ss, n.censor > 0), shape = 3, size = 2,
show.legend = FALSE) +
scale_linetype_manual(values = c(dense = "solid", sparse = "dashed"),
breaks = "sparse",
labels = sprintf("< %g%% at risk", 100 * frac),
name = NULL) +
coord_cartesian(ylim = c(0, 1)) +
labs(x = "Time", y = "Survival probability", ...) +
theme_survminer()
}fit <- survfit(Surv(time, status) ~ sex, data = lung)
ggsurv_dashed_below(fit, lung, frac = 0.10) +
scale_colour_discrete(name = "Sex", labels = c("Male", "Female"))Each curve turns dashed at its own crossing point. A few adjustments:
max(n.risk) per stratum). To key it to the pooled study
size instead, replace max(n.risk) with the overall starting
n and every curve dashes at the same count.frac is an argument, so the same helper produces the
guideline’s 10% cutoff or any other fraction.linetype, so it
does not go through ggsurvplot() directly. If you also want
the number-at-risk table underneath, build it separately with
ggrisktable() on the same fit and stack the two (for
example with
ggpubr::ggarrange(..., ncol = 1, heights = c(3, 1))).By default the legend key is a straight horizontal line. Some
publications draw the key as a small step that echoes the Kaplan-Meier
shape (survminer issue #537).
ggsurvplot() returns the curve as $plot, whose
first layer draws the survival lines, so you can swap that layer’s
legend key for a custom step glyph – overriding only the key, leaving
the curve drawing untouched.
library(grid)
# Draw the legend key of a ggsurvplot's curves as a small step (stair) that
# echoes the Kaplan-Meier shape, instead of the default straight line.
step_legend_key <- function(p) {
`%||%` <- function(x, y) if (is.null(x)) y else x # base R has this since 4.4
draw_key_step <- function(data, params, size) {
grid::linesGrob(
x = c(0, 0.45, 0.45, 1), y = c(0.25, 0.25, 0.7, 0.7),
gp = grid::gpar(col = data$colour %||% "black",
lwd = (data$linewidth %||% data$size %||% 0.75) * .pt,
lty = data$linetype %||% 1, lineend = "butt"))
}
# the survival curves are the first layer; subclass its geom to override ONLY
# the legend key, so the curve drawing itself is unchanged
g <- p$plot$layers[[1]]$geom
p$plot$layers[[1]]$geom <-
ggplot2::ggproto(paste0(class(g)[1], "Key"), g, draw_key = draw_key_step)
p
}fit2 <- survfit(Surv(time, status) ~ sex, data = lung)
p <- ggsurvplot(fit2, data = lung, palette = c("#E7302A", "#2E43F3"),
legend.title = "Group", legend.labs = c("High", "Low"))
step_legend_key(p)Adjust the x/y coordinates inside
draw_key_step() to change the step shape (for example a
taller riser, or a two-step stair).
surv.median.line = "hv" already drops a line at each
group’s median survival. To also show the confidence interval
of the median (survminer issue #345), use
surv_median(), which returns median,
lower and upper per group, and add a
horizontal band at survival = 0.5.
add_median_ci <- function(p, fit, y = 0.5) {
med <- surv_median(fit)
# one colour per group, read from the drawn curve layer (no internal helpers)
b <- ggplot2::ggplot_build(p$plot)$data[[1]]
cols <- b$colour[match(sort(unique(b$group)), b$group)]
# drop groups whose median is not reached (surv_median() returns NA there)
keep <- !is.na(med$median)
med <- med[keep, ]
cols <- cols[keep]
p$plot +
geom_segment(data = med, aes(x = lower, xend = upper, y = y, yend = y),
colour = cols, linewidth = 0.9, inherit.aes = FALSE) +
geom_point(data = med, aes(x = median, y = y),
colour = cols, size = 2.4, inherit.aes = FALSE)
}fit3 <- survfit(Surv(time, status) ~ sex, data = lung)
p <- ggsurvplot(fit3, data = lung, palette = c("#E7302A", "#2E43F3"),
surv.median.line = "hv",
legend.labs = c("Male", "Female"), legend.title = "Sex")
p$plot <- add_median_ci(p, fit3)
pThe horizontal bar spans the 95% CI of the median survival time
(surv_median()’s lower/upper);
the point marks the median itself. A group whose survival never drops to
0.5 has no median (surv_median() returns NA),
so it is simply left without a bar.
legend.labs = c(expression(...)) alone does not render
as math, and + does not apply to the compound
ggsurvplot object. Add the scale to $plot (a
plain ggplot) instead, supplying the same palette so the colours are
unchanged (survminer issue #350).
pal <- c("#E7302A", "#2E43F3")
p <- ggsurvplot(fit3, data = lung, palette = pal, legend.title = "Group")
p$plot <- p$plot +
scale_colour_manual(values = pal,
labels = c(expression(x^2 ~ y^2), expression(alpha[i] >= beta)),
name = "Group")
p(The “Scale for colour is already present” message is expected – we are replacing survminer’s default colour scale with one that carries the math labels.)
With a risk table (risk.table = TRUE), the table’s
strata labels are a separate discrete y-axis, so they keep the raw
sex=1 / sex=2 names unless you relabel them
too. Carry the same expressions to $table with
scale_y_discrete(). The table lists the strata
top-to-bottom but the discrete y-axis orders its breaks bottom-to-top,
so pass the labels reversed to keep each label on its own row:
math.labs <- c(expression(x^2 ~ y^2), expression(alpha[i] >= beta))
p <- ggsurvplot(fit3, data = lung, palette = pal, risk.table = TRUE,
legend.title = "Group")
p$plot <- p$plot +
scale_colour_manual(values = pal, labels = math.labs, name = "Group")
p$table <- p$table + scale_y_discrete(labels = rev(math.labs))
pCox-model adjusted curves are drawn with
survfit(coxph_fit, newdata = ...), one curve per covariate
profile. Asking for risk.table = TRUE there gives a single
pooled row, because those predicted curves have no per-group number at
risk – the profiles are hypothetical covariate combinations, not subject
groups (survminer issue #231).
When the profiles correspond to a real grouping
variable (here sex), the meaningful at-risk counts
are the Kaplan-Meier ones for that variable. Build them with
ggrisktable() from a KM fit and drop them into the compound
object’s $table slot, matching legend.labs and
sharing an xlim so the two panels line up:
km.fit <- survfit(Surv(time, status) ~ sex, data = lung)
res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
sex_df <- with(lung, data.frame(sex = c(1, 2),
age = mean(age, na.rm = TRUE), ph.ecog = c(1, 1)))
fit <- survfit(res.cox, newdata = sex_df)
# note: a survfit() from a coxph fit carries no data, so pass `data =` explicitly
cox.surv <- ggsurvplot(fit, data = sex_df, conf.int = TRUE,
legend.labs = c("sex=1", "sex=2"))
cox.surv$table <- ggrisktable(km.fit, data = lung, color = "strata",
legend.labs = c("sex=1", "sex=2"),
xlim = c(0, max(fit$time)))
print(cox.surv, risk.table.height = 0.28)Be careful with the meaning: the table is the Kaplan-Meier
number at risk for the grouping variable, not an at-risk count
for the model-based adjusted curves (adjusted curves are expectations
for hypothetical profiles and have no literal number at risk). It is
honest only when each newdata row is a real subgroup of the
data (e.g. sex = 1 / sex = 2); for arbitrary
covariate combinations there is no per-profile risk table and this
overlay would mislead.
ggadjustedcurves()ggadjustedcurves() draws Cox-adjusted curves by a
grouping variable and returns a plain ggplot, so
risk.table = TRUE has no effect there (survminer issue #286). The
same idea as the previous recipe applies, the other way round: build a
normal Kaplan-Meier ggsurvplot() (which comes with a
correctly aligned number-at-risk table for the grouping variable) and
swap its $plot for the adjusted-curves plot, matching
palette, xlim and x-breaks so the panels line up.
lung2 <- transform(lung, sex = factor(sex, labels = c("Male", "Female")))
pal <- c("#E7302A", "#2E43F3")
xmax <- max(lung2$time, na.rm = TRUE)
# a KM ggsurvplot -> aligned curve + number-at-risk table by group
km.surv <- ggsurvplot(survfit(Surv(time, status) ~ sex, data = lung2), data = lung2,
risk.table = TRUE, palette = pal,
legend.labs = c("Male", "Female"), legend.title = "Sex",
break.time.by = 250, xlim = c(0, xmax))
# Cox-adjusted curves (a plain ggplot); match palette / xlim / breaks
res.cox <- coxph(Surv(time, status) ~ sex + age + ph.ecog, data = lung2)
adj <- ggadjustedcurves(res.cox, variable = "sex", data = lung2, palette = pal) +
coord_cartesian(xlim = c(0, xmax)) +
scale_x_continuous(breaks = seq(0, xmax, 250)) +
labs(x = "Time", y = "Adjusted survival probability") +
guides(colour = guide_legend("Sex"), fill = guide_legend("Sex")) +
theme_survminer()
# swap the KM curve for the adjusted curve; keep the aligned KM risk table
km.surv$plot <- adj
print(km.surv, risk.table.height = 0.26)The same honesty caveat applies: the table is the KM number at risk for the grouping variable, not for the adjusted curves, and is meaningful only when the adjustment variable is a real group.
For a cmprsk::cuminc() competing-risks fit, two common
requests are to plot only the event of interest
(dropping the competing events) and to add a number-at-risk
table underneath (survminer issue #223).
Both compose from the cuminc object and a companion
survfit(): the cumulative-incidence estimates (and their
variance) for the event of interest come straight out of
cuminc, and the number at risk – subjects still event-free
and under observation, the usual convention for a competing-risks figure
– is the n.risk of a survfit() that treats
any event as the endpoint.
library(cmprsk)
# competing-risks data: event 1 = PCM (of interest), 2 = death (competing)
mg <- transform(mgus2,
etime = ifelse(pstat == 1, ptime, futime),
event = ifelse(pstat == 1, 1, 2 * death))
pal <- c("#E7302A", "#2E43F3"); xmax <- max(mg$etime); eoi <- "1"
ci <- cuminc(mg$etime, mg$event, group = mg$sex)
# cumulative incidence of the event of interest, per group, with 95% CI bands
comps <- names(ci)[grepl(paste0(" ", eoi, "$"), names(ci))]
cif <- do.call(rbind, lapply(comps, function(nm) {
d <- as.data.frame(ci[[nm]][c("time", "est", "var")])
d$group <- sub(paste0(" ", eoi, "$"), "", nm); d
}))
cif$lower <- pmax(0, cif$est - 1.96 * sqrt(cif$var))
cif$upper <- pmin(1, cif$est + 1.96 * sqrt(cif$var))
p_cif <- ggplot(cif, aes(time, est, colour = group, fill = group)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.15, colour = NA) +
geom_step(linewidth = 0.9) +
scale_colour_manual(values = pal) + scale_fill_manual(values = pal) +
scale_x_continuous(breaks = seq(0, xmax, 60)) +
coord_cartesian(xlim = c(0, xmax), ylim = c(0, 1)) +
labs(x = "Time (months)", y = "Cumulative incidence of PCM") +
guides(colour = guide_legend("Sex"), fill = guide_legend("Sex")) +
theme_survminer()
# number-at-risk table: survfit on "any event" gives the event-free at-risk count
km.surv <- ggsurvplot(survfit(Surv(etime, event > 0) ~ sex, data = mg), data = mg,
risk.table = TRUE, palette = pal, legend.labs = c("F", "M"),
legend.title = "Sex", break.time.by = 60, xlim = c(0, xmax))
km.surv$plot <- p_cif # keep the aligned number-at-risk table
print(km.surv, risk.table.height = 0.26)To pick a different event of interest, change eoi (the
numeric fstatus code as a string). The confidence bands use
cuminc’s var; note the at-risk counts are the
standard event-free-under-observation convention, not a competing-risks
weighted count.
ggforest() draws one fitted model. For the common
screening step of showing the hazard ratio from a separate
univariate Cox model for each candidate covariate on one forest
plot (survminer issues #459, #310), fit
one model per covariate, collect the HR / 95% CI / p from each, and draw
the forest with geom_pointrange() on a log scale.
d <- lung
d$sex <- factor(d$sex, labels = c("Male", "Female"))
d$ph.ecog <- factor(d$ph.ecog)
vars <- c("sex", "age", "ph.ecog", "wt.loss") # candidate covariates
# one univariate Cox model per covariate -> HR, 95% CI, p (one row per coefficient)
uni <- do.call(rbind, lapply(vars, function(v) {
s <- summary(coxph(reformulate(v, response = "Surv(time, status)"), data = d))
data.frame(variable = v, term = rownames(s$conf.int),
HR = s$conf.int[, "exp(coef)"],
lower = s$conf.int[, "lower .95"],
upper = s$conf.int[, "upper .95"],
p = s$coefficients[, "Pr(>|z|)"], row.names = NULL)
}))
# row label: variable name, plus the factor level for multi-level terms
level <- mapply(function(v, t) sub(paste0("^", v), "", t), uni$variable, uni$term)
uni$label <- ifelse(level == "", uni$variable, paste0(uni$variable, ": ", level))
uni$label <- factor(uni$label, levels = rev(uni$label))
ggplot(uni, aes(HR, label)) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
geom_pointrange(aes(xmin = lower, xmax = upper)) +
geom_text(aes(label = sprintf("%.2f (%.2f-%.2f)", HR, lower, upper)),
vjust = -0.8, size = 3) +
scale_x_log10() +
labs(x = "Hazard ratio (95% CI, log scale)", y = NULL,
title = "Univariate Cox regression") +
theme_survminer()Each row is its own univariate model (not a single multivariable
fit); a multi-level factor contributes one row per level. Swap in a
multivariable
ggforest(coxph(Surv(time, status) ~ ., data = ...)) when
you want one adjusted model instead.
A survreg() accelerated-failure-time fit predicts a
smooth survival curve; feeding its predictions to
ggsurvplot() as a data frame draws them as a step function
by default. Use surv.geom = geom_line to draw the curve as
a smooth line instead (survminer issue #276). A
survreg curve is only defined once the covariates are
fixed, so choose the profile(s) to draw:
# a parametric (Weibull) AFT model with a covariate
wb <- survreg(Surv(time, status) ~ sex, data = lung)
# predict(type = "quantile", p = s) returns the time at which survival = 1 - s,
# i.e. the inverse of the fitted survival function; pick the profiles to show.
s <- seq(0.01, 0.99, by = 0.01)
prof <- data.frame(sex = c(1, 2))
pred <- predict(wb, newdata = prof, type = "quantile", p = s) # rows = profiles
wbdf <- data.frame(
time = as.vector(t(pred)),
surv = rep(1 - s, times = nrow(prof)),
strata = factor(rep(c("sex=1", "sex=2"), each = length(s))),
upper = NA, lower = NA, std.err = NA # a model curve carries no CI band
)
pal <- c("#E7302A", "#2E43F3")
# surv.geom = geom_line -> smooth line, not a step; conf.int = FALSE is required.
ggsurvplot(wbdf, conf.int = FALSE, surv.geom = geom_line, palette = pal,
legend.title = "Weibull", legend.labs = c("sex=1", "sex=2"))To judge how well the parametric model fits, overlay the curves (dashed) on the empirical Kaplan-Meier:
km <- ggsurvplot(survfit(Surv(time, status) ~ sex, data = lung), data = lung,
palette = pal, legend.labs = c("sex=1", "sex=2"))
km$plot + geom_line(data = wbdf, aes(time, surv, colour = strata),
linetype = "dashed", linewidth = 1)The curve is the model’s prediction at the covariate profile you
choose – for a multi-covariate model, fix the other covariates in
newdata and decide which one indexes the curves. It is a
deterministic parametric fit, not an empirical estimate, so it has no
confidence band (upper/lower are
NA and conf.int = FALSE is required). For a
Weibull fit over the Kaplan-Meier with built-in confidence bands,
ggflexsurvplot() (on a flexsurv::flexsurvreg()
model) is an alternative.
ggadjustedcurves()