Bayesian analysis with Stan

Bayesian fits use the same Stan likelihood as MLE but sample from the posterior with rstan::sampling.

Before you model

Brief exploratory context for the larynx data (see Getting started for full EDA):

library(spsurv)
library(generics)
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) +
  labs(x = "Time (years)", y = "Survival probability", color = "Stage") +
  theme_bw() +
  theme(legend.position = "bottom")

Kaplan-Meier by stage (unadjusted).

What to decide

Stage separation motivates the same formula used in the getting-started vignette.
Small sample → posteriors may be wide; priors matter more.

Fit and summarize

fit <- bpph(
  Surv(time, delta) ~ age + stage,
  degree = 5,
  data = larynx,
  approach = "bayes",
  iter = 400,
  warmup = 200,
  chains = 1,
  cores = 1,
  init = 0,
  priors = list(
    beta = c("normal(0,4)"),
    gamma = c("lognormal(0,4)")
  )
)
#> 
#> SAMPLING FOR MODEL 'spbp' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 3.5e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.35 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 400 [  0%]  (Warmup)
#> Chain 1: Iteration:  40 / 400 [ 10%]  (Warmup)
#> Chain 1: Iteration:  80 / 400 [ 20%]  (Warmup)
#> Chain 1: Iteration: 120 / 400 [ 30%]  (Warmup)
#> Chain 1: Iteration: 160 / 400 [ 40%]  (Warmup)
#> Chain 1: Iteration: 200 / 400 [ 50%]  (Warmup)
#> Chain 1: Iteration: 201 / 400 [ 50%]  (Sampling)
#> Chain 1: Iteration: 240 / 400 [ 60%]  (Sampling)
#> Chain 1: Iteration: 280 / 400 [ 70%]  (Sampling)
#> Chain 1: Iteration: 320 / 400 [ 80%]  (Sampling)
#> Chain 1: Iteration: 360 / 400 [ 90%]  (Sampling)
#> Chain 1: Iteration: 400 / 400 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.2 seconds (Warm-up)
#> Chain 1:                0.161 seconds (Sampling)
#> Chain 1:                0.361 seconds (Total)
#> Chain 1:
summary(fit)
#> Call:
#> bpph(formula = Surv(time, delta) ~ age + stage, degree = 5, data = larynx, 
#>     approach = "bayes", iter = 400, warmup = 200, chains = 1, 
#>     cores = 1, init = 0, priors = list(beta = c("normal(0,4)"), 
#>         gamma = c("lognormal(0,4)")), model = "ph")
#> 
#> Bayesian Bernstein PH model: 
#> Regression coefficients:
#>        Estimate    2.5%   97.5% Std. Error
#> age      0.0213 -0.0061  0.0500        0.0
#> stage2   0.1945 -0.5993  1.1772        0.5
#> stage3   0.7536  0.1206  1.5273        0.4
#> stage4   1.7399  0.9831  2.6690        0.4
#> 
#> Exponentiated coefficients:
#>        Estimate 2.5% 97.5%
#> age        1.02 0.99   1.1
#> stage2     1.36 0.27   2.7
#> stage3     2.26 0.82   3.9
#> stage4     6.27 1.97  12.6
#> 
#> --- 
#> DIC = 296   WAIC = -149

fit_mle <- bpph(
  Surv(time, delta) ~ age + stage,
  degree = 5,
  data = larynx,
  approach = "mle",
  init = 0
)

Visual check — posterior forest plot

td_b <- tidy(fit, conf.int = TRUE, exponentiate = TRUE)
hpd <- credint(fit, prob = 0.95, type = "HPD")
td_b$hpd.low <- exp(hpd[, 1])
td_b$hpd.high <- exp(hpd[, 2])
td_b$term <- factor(td_b$term, levels = rev(td_b$term))

ggplot(td_b, aes(x = estimate, y = term)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
  geom_pointrange(aes(xmin = hpd.low, xmax = hpd.high), linewidth = 0.4, color = "steelblue") +
  labs(x = "Hazard ratio (posterior median)", y = NULL, title = "95% HPD intervals") +
  theme_bw()

Posterior medians with 95% HPD intervals (exponentiated).

Visual check — MLE vs Bayes

td_m <- tidy(fit_mle, conf.int = TRUE, exponentiate = TRUE)
cmp <- merge(
  td_m[, c("term", "estimate", "conf.low", "conf.high")],
  td_b[, c("term", "estimate", "hpd.low", "hpd.high")],
  by = "term",
  suffixes = c("_mle", "_bayes")
)
cmp_long <- rbind(
  data.frame(term = cmp$term, method = "MLE", low = cmp$conf.low, high = cmp$conf.high,
             est = cmp$estimate_mle),
  data.frame(term = cmp$term, method = "Bayes", low = cmp$hpd.low, high = cmp$hpd.high,
             est = cmp$estimate_bayes)
)
cmp_long$term <- factor(cmp_long$term, levels = rev(unique(cmp$term)))

ggplot(cmp_long, aes(x = est, y = term, color = method)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey50") +
  geom_pointrange(aes(xmin = low, xmax = high), position = position_dodge(width = 0.5), linewidth = 0.4) +
  labs(x = "Hazard ratio", y = NULL, color = NULL) +
  theme_bw()

MLE 95% CI vs Bayesian 95% HPD (exponentiated).

Why this is insightful

Posterior intervals show whether data overwhelm priors. Agreement between MLE and Bayes supports data-driven conclusions; wide HPDs flag weak information.

What to decide

Posterior ≈ MLE → conclusion robust to prior choice.
Wide HPD + low ESS → do not report precise HRs; tune Stan (adapt_delta, more iter).
WAIC/DIC for comparing models within the same family at the same degree.

Credible intervals and criteria

credint(fit, prob = 0.95, type = "HPD")
#>               lower      upper
#> age    -0.006129856 0.05002859
#> stage2 -0.599284367 1.17716377
#> stage3  0.120596698 1.52734682
#> stage4  0.983148869 2.66902842
#> attr(,"Probability")
#> [1] 0.95

sm <- summary(fit)
c(DIC = sm$dic, WAIC = sm$waic, LPML = sm$lpml)
#>       DIC      WAIC      LPML 
#>  295.8220 -148.6405 -148.7211

Convergence (optional deep dive)

Install bayesplot for trace and pairs plots:

library(bayesplot)
mcmc_pairs(fit$posterior$beta, off_diag_fun = "hex")

Always inspect divergences, split R-hat, and ESS before interpreting results.

tidybayes posterior draws

Install posterior and tidybayes for long-format draws and interval plots. After library(spsurv), spread_draws() and tidy_draws() dispatch on Bayes fits when tidybayes is loaded.

library(tidybayes)
dr <- as_draws_df.spbp(fit)
spread_draws(fit, `beta[age]`)
#> # A tibble: 200 × 4
#>    .chain .iteration .draw `beta[age]`
#>     <int>      <int> <int>       <dbl>
#>  1      1          1     1     0.0756 
#>  2      1          2     2     0.0251 
#>  3      1          3     3     0.00827
#>  4      1          4     4     0.0417 
#>  5      1          5     5     0.0394 
#>  6      1          6     6     0.0251 
#>  7      1          7     7     0.0196 
#>  8      1          8     8     0.0184 
#>  9      1          9     9     0.0237 
#> 10      1         10    10     0.0356 
#> # ℹ 190 more rows

Draw-level survival curves:

head(spread_surv_draws.spbp(fit, times = c(1, 2, 3), newdata = larynx[1, ]))
#>     stage time age diagyr delta .chain .iteration .draw time      surv id
#> 1       1  0.6  77     76     1      1          1     1    1 0.8920691  1
#> 1.1     1  0.6  77     76     1      1          1     1    2 0.8188651  1
#> 1.2     1  0.6  77     76     1      1          1     1    3 0.7591480  1
#> 1.3     1  0.6  77     76     1      1          2     2    1 0.9089544  1
#> 1.4     1  0.6  77     76     1      1          2     2    2 0.8008039  1
#> 1.5     1  0.6  77     76     1      1          2     2    3 0.7007909  1

See vignette("tidymodels", package = "spsurv") for workflows and parsnip engines.

Survival prediction

Posterior mean survival is a Bernstein polynomial in time, so the curve is smooth (geom_line). Kaplan–Meier in the EDA plot above remains a step function.

pr <- predict(fit, times = seq(0, max(larynx$time), length.out = 121))
head(pr)
#>   id       time      surv     lower     upper     cumhaz     std.err
#> 1  1 0.00000000 1.0000000 1.0000000 1.0000000 0.00000000 0.000000000
#> 2  1 0.08916667 0.9887786 0.9815992 0.9946361 0.01129081 0.003436125
#> 3  1 0.17833333 0.9778741 0.9639561 0.9889651 0.02239698 0.006591042
#> 4  1 0.26750000 0.9672653 0.9470323 0.9830119 0.03333048 0.009488356
#> 5  1 0.35666667 0.9569322 0.9329433 0.9785885 0.04410319 0.012150139
#> 6  1 0.44583333 0.9468562 0.9180277 0.9720775 0.05472681 0.014597013
ggplot(pr, aes(x = time, y = surv)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.6) +
  labs(x = "Time (years)", y = "Survival probability") +
  theme_bw()

Smooth posterior survival at the mean covariate profile with 95% HPD ribbon.

See the Survival prediction and ggplot vignette for covariate-specific ribbons.