---
title: "Customization recipes"
author: "survminer"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Customization recipes}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  comment = "#>",
  collapse = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 5
)
```

`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.

```{r setup}
library(survminer)
library(survival)
library(dplyr)
library(ggplot2)
```

# Dashed curves once few patients remain at risk

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](https://pubmed.ncbi.nlm.nih.gov/18355567/) recommend a dashed
line once the number at risk falls below 10% of the initial cohort (survminer
issue [#559](https://github.com/kassambara/survminer/issues/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
`linetype`s. Repeating the boundary point in the solid part makes the two step
segments meet with no gap.

```{r dashed-helper}
# 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()
}
```

```{r dashed-plot}
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:

- The threshold is each group's own initial cohort (`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.
- This rebuilds the line layer to control `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))`).

# Step-shaped legend keys

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](https://github.com/kassambara/survminer/issues/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.

```{r step-key-helper}
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
}
```

```{r step-key-plot}
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).

# Mark the median survival and its confidence interval

`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](https://github.com/kassambara/survminer/issues/345)), use
`surv_median()`, which returns `median`, `lower` and `upper` per group, and add
a horizontal band at survival = 0.5.

```{r median-ci}
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)
}
```

```{r median-ci-plot}
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)
p
```

The 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.

# Plot math (expressions) in the legend

`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](https://github.com/kassambara/survminer/issues/350)).

```{r plotmath-legend}
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:

```{r plotmath-legend-risktable}
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))
p
```

# One shared legend across a grid of plots

Several separate survival analyses that use the same grouping are often shown in
a grid with a single common legend (survminer issue
[#340](https://github.com/kassambara/survminer/issues/340)). Pull each
`ggsurvplot`'s `$plot` and arrange them with `ggpubr::ggarrange()`, which draws
one shared legend and keeps every panel the same size.

```{r shared-legend, fig.height = 6}
colon2  <- colon[!is.na(colon$sex), ]
subsets <- split(colon2, colon2$extent)
subsets <- subsets[sapply(subsets, nrow) > 30]
subsets <- subsets[seq_len(min(4, length(subsets)))]

plots <- lapply(names(subsets), function(nm) {
  d   <- subsets[[nm]]
  fit <- survfit(Surv(time, status) ~ sex, data = d)
  ggsurvplot(fit, data = d, legend.title = "Sex",
             legend.labs = c("Male", "Female"), title = paste("extent", nm))$plot
})

ggpubr::ggarrange(plotlist = plots, ncol = 2, nrow = 2,
                  common.legend = TRUE, legend = "bottom")
```

# Risk table under Cox-adjusted survival curves

Cox-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](https://github.com/kassambara/survminer/issues/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:

```{r cox-risktable, fig.height = 6}
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.

# Risk table under `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](https://github.com/kassambara/survminer/issues/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.

```{r ggadjustedcurves-risktable, fig.height = 6}
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.

# Single event and a number-at-risk table for competing risks

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](https://github.com/kassambara/survminer/issues/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.

```{r cuminc-risktable, fig.height = 6, eval = requireNamespace("cmprsk", quietly = TRUE)}
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.

# Forest plot of several univariate Cox models

`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](https://github.com/kassambara/survminer/issues/459),
[#310](https://github.com/kassambara/survminer/issues/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.

```{r univariate-forest, fig.height = 4.5}
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.

# Smooth parametric (survreg / Weibull) survival curves

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](https://github.com/kassambara/survminer/issues/276)). A `survreg` curve is
only defined once the covariates are fixed, so choose the profile(s) to draw:

```{r survreg-smooth, fig.height = 4}
# 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:

```{r survreg-overlay, fig.height = 4}
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.
