Survival prediction and ggplot workflows

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Unlike stepwise Cox survfit curves, spsurv evaluates survival on a dense time grid using the fitted Bernstein baseline. Kaplan–Meier stays a step function; Bernstein-polynomial curves are always smooth lines (geom_line, or plot(survfit(fit))). This vignette connects observed KM to model-based predictions.

Before you model

library(spsurv)
library(KMsurv)
library(survival)
library(ggplot2)
data(larynx)
larynx$stage <- factor(larynx$stage)

km_stage <- survfit(Surv(time, delta) ~ stage, data = larynx)
km_long <- data.frame(
  time = km_stage$time,
  surv = km_stage$surv,
  stage = rep(levels(larynx$stage), km_stage$strata)
)
ggplot(km_long, aes(x = time, y = surv, color = stage)) +
  geom_step(linewidth = 0.6) +
  geom_point(
    data = subset(larynx, delta == 0),
    aes(x = time, y = 0, color = stage),
    shape = 3,
    size = 1.5,
    alpha = 0.5,
    inherit.aes = FALSE
  ) +
  labs(x = "Time (years)", y = "Survival probability", color = "Stage") +
  theme_bw() +
  theme(legend.position = "bottom")

Kaplan-Meier by stage with censoring marks (+).

We compare predictions at age = 70 across all stages — a clinically interpretable profile that holds age fixed while varying disease stage.

Why this is insightful

KM curves show unadjusted stage separation. Censoring marks (+) remind you that late follow-up estimates rest on fewer subjects.

What to decide

Clear stage ordering in KM → including stage in the model is justified.
Crossing KM curves → proportional hazards may be questionable; see the Choosing PH, PO, and AFT models vignette.

Fit and summarize

fit <- bpph(
  Surv(time, delta) ~ age + stage,
  degree = 5,
  data = larynx,
  approach = "mle",
  init = 0
)
summary(fit)
#> Call:
#> bpph(formula = Surv(time, delta) ~ age + stage, degree = 5, data = larynx, 
#>     approach = "mle", init = 0, model = "ph")
#> 
#> Bernstein PH model: 
#> Regression coefficients:
#>        Estimate    2.5%   97.5% Std. Error z value Pr(>|z|)    
#> age      0.0197 -0.0085  0.0478     0.0143     1.4     0.17    
#> stage2   0.1730 -0.7324  1.0783     0.4619     0.4     0.71    
#> stage3   0.6521 -0.0450  1.3492     0.3557     1.8     0.07 .  
#> stage4   1.7778  0.9471  2.6086     0.4239     4.2    3e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Exponentiated coefficients:
#>        Estimate 2.5% 97.5%
#> age        1.02 0.99   1.0
#> stage2     1.19 0.48   2.9
#> stage3     1.92 0.96   3.9
#> stage4     5.92 2.58  13.6
#> 
#> --- 
#> loglik = -141   AIC = 299

Prediction grid

plot_times <- seq(0, max(larynx$time), length.out = 121)
newdata <- data.frame(age = 70, stage = levels(larynx$stage))
newdata
#>   age stage
#> 1  70     1
#> 2  70     2
#> 3  70     3
#> 4  70     4

`predict()` — tidy data frame

pr <- predict(fit, newdata = newdata, times = plot_times)
head(pr)
#>     age stage id       time      surv     lower upper      cumhaz    std.err
#> 1    70     1  1 0.00000000 1.0000000 1.0000000     1 0.000000000 0.00000000
#> 1.1  70     1  1 0.08916667 0.9922984 0.9666007     1 0.007731445 0.01338715
#> 1.2  70     1  1 0.17833333 0.9847620 0.9390056     1 0.015355286 0.02427520
#> 1.3  70     1  1 0.26750000 0.9773821 0.9160002     1 0.022877617 0.03309298
#> 1.4  70     1  1 0.35666667 0.9701500 0.8964400     1 0.030304583 0.04031670
#> 1.5  70     1  1 0.44583333 0.9630573 0.8792141     1 0.037642381 0.04647250

Visual check — KM vs model overlay

Panel A: observed KM (steps) vs model prediction (smooth dashed) per stage.

pr$stage <- newdata$stage[match(as.character(pr$id), as.character(seq_len(nrow(newdata))))]

ggplot() +
  geom_step(
    data = km_long,
    aes(x = time, y = surv, color = stage),
    linewidth = 0.5,
    alpha = 0.8
  ) +
  geom_line(
    data = pr,
    aes(x = time, y = surv, color = stage),
    linetype = "dashed",
    linewidth = 0.6
  ) +
  labs(x = "Time (years)", y = "Survival probability", color = "Stage") +
  theme_bw() +
  theme(legend.position = "bottom")

Observed KM (steps) vs smooth Bernstein prediction at age 70 (dashed).

Panel B: model-based ribbons (adjusted prediction with uncertainty), still drawn with geom_line rather than geom_step.

ggplot(pr, aes(x = time, y = surv, color = stage, fill = stage)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.5) +
  labs(x = "Time (years)", y = "Survival probability", color = NULL, fill = NULL) +
  theme_bw() +
  theme(legend.position = "bottom")

Smooth stage-specific Bernstein survival with 95% intervals.

Why this is insightful

Panel A connects the smooth BP fit to observed risk sets. Panel B shows covariate-adjusted uncertainty — ribbons often widen late in follow-up when data are sparse.

What to decide

Model curves systematically below KM → check covariate adjustment or misspecification.
Wide bands late in follow-up → avoid over-interpreting tail survival.
For publication → export predict() output directly to your ggplot pipeline.

`survfit()` with `tidy = TRUE`

sf <- survfit(fit, newdata = newdata, times = plot_times, tidy = TRUE)
head(sf)
#>     age stage id       time      surv     lower upper      cumhaz    std.err
#> 1    70     1  1 0.00000000 1.0000000 1.0000000     1 0.000000000 0.00000000
#> 1.1  70     1  1 0.08916667 0.9922984 0.9666007     1 0.007731445 0.01338715
#> 1.2  70     1  1 0.17833333 0.9847620 0.9390056     1 0.015355286 0.02427520
#> 1.3  70     1  1 0.26750000 0.9773821 0.9160002     1 0.022877617 0.03309298
#> 1.4  70     1  1 0.35666667 0.9701500 0.8964400     1 0.030304583 0.04031670
#> 1.5  70     1  1 0.44583333 0.9630573 0.8792141     1 0.037642381 0.04647250

Bayesian fits

For approach = "bayes", credible bands use posterior draws with interval.type = "hpd" and monotone = TRUE (default for Bayes).

fit_bayes <- bpph(
  Surv(time, delta) ~ age + stage,
  degree = 5,
  data = larynx,
  approach = "bayes",
  iter = 400,
  chains = 1,
  cores = 1
)
predict(fit_bayes, newdata = newdata, times = plot_times, interval.type = "hpd")

Interval types (MLE)

Log–log intervals are undefined at S = 1 (time = 0), so we evaluate at time > 0 only:

loglog_times <- plot_times[plot_times > 0]
head(predict(fit, newdata = newdata, times = loglog_times, type = "log-log"))
#>     age stage id       time      surv     lower     upper      cumhaz
#> 1    70     1  1 0.08916667 0.9922984 0.7943630 0.9997404 0.007731445
#> 1.1  70     1  1 0.17833333 0.9847620 0.7115221 0.9993075 0.015355286
#> 1.2  70     1  1 0.26750000 0.9773821 0.6772829 0.9986577 0.022877617
#> 1.3  70     1  1 0.35666667 0.9701500 0.6629301 0.9977685 0.030304583
#> 1.4  70     1  1 0.44583333 0.9630573 0.6549430 0.9966575 0.037642381
#> 1.5  70     1  1 0.53500000 0.9560957 0.6459815 0.9953977 0.044897252
#>        std.err
#> 1   0.01338715
#> 1.1 0.02427520
#> 1.2 0.03309298
#> 1.3 0.04031670
#> 1.4 0.04647250
#> 1.5 0.05212577

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.

Survival prediction and ggplot workflows

Before you model

Why this is insightful

What to decide

Fit and summarize

Prediction grid

predict() — tidy data frame

Visual check — KM vs model overlay

Why this is insightful

What to decide

survfit() with tidy = TRUE

Bayesian fits

Interval types (MLE)

`predict()` — tidy data frame

`survfit()` with `tidy = TRUE`