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.4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.34 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.19 seconds (Warm-up)
#> Chain 1:                0.101 seconds (Sampling)
#> Chain 1:                0.291 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.02009 -0.00093  0.05001        0.0
#> stage2  0.11440 -1.02630  1.02433        0.5
#> stage3  0.63381  0.03522  1.23235        0.3
#> stage4  1.70501  0.91802  2.65480        0.4
#> 
#> Exponentiated coefficients:
#>        Estimate 2.5% 97.5%
#> age        1.02 1.00   1.1
#> stage2     1.27 0.36   2.8
#> stage3     2.00 1.04   3.4
#> stage4     6.05 1.30  11.4
#> 
#> --- 
#> DIC = 296   WAIC = -148

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.0009266491 0.05001317
#> stage2 -1.0262966911 1.02432978
#> stage3  0.0352201064 1.23235423
#> stage4  0.9180206255 2.65479802
#> attr(,"Probability")
#> [1] 0.95

sm <- summary(fit)
c(DIC = sm$dic, WAIC = sm$waic, LPML = sm$lpml)
#>       DIC      WAIC      LPML 
#>  296.2896 -148.3970 -148.3833

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.0229
#>  2      1          2     2      0.0299
#>  3      1          3     3      0.0288
#>  4      1          4     4      0.0143
#>  5      1          5     5      0.0273
#>  6      1          6     6      0.0315
#>  7      1          7     7      0.0123
#>  8      1          8     8      0.0237
#>  9      1          9     9      0.0187
#> 10      1         10    10     -0.0184
#> # ℹ 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.9466978  1
#> 1.1     1  0.6  77     76     1      1          1     1    2 0.8945650  1
#> 1.2     1  0.6  77     76     1      1          1     1    3 0.8343935  1
#> 1.3     1  0.6  77     76     1      1          2     2    1 0.9207739  1
#> 1.4     1  0.6  77     76     1      1          2     2    2 0.8540969  1
#> 1.5     1  0.6  77     76     1      1          2     2    3 0.7844553  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.9890002 0.9819762 0.9958358 0.01106839 0.003863031
#> 3  1 0.17833333 0.9782980 0.9646813 0.9913873 0.02196982 0.007429415
#> 4  1 0.26750000 0.9678711 0.9480810 0.9865386 0.03271780 0.010722792
#> 5  1 0.35666667 0.9576989 0.9321426 0.9818371 0.04332550 0.013765001
#> 6  1 0.44583333 0.9477619 0.9168340 0.9773529 0.05380576 0.016576190
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.