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 2.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 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.156 seconds (Warm-up)
#> Chain 1:                0.103 seconds (Sampling)
#> Chain 1:                0.259 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.019 -0.006  0.048        0.0
#> stage2    0.161 -0.697  1.023        0.5
#> stage3    0.660 -0.014  1.418        0.4
#> stage4    1.753  0.799  2.570        0.5
#> 
#> Exponentiated coefficients:
#>        Estimate 2.5% 97.5%
#> age        1.02 0.99   1.0
#> stage2     1.31 0.42   2.7
#> stage3     2.08 0.75   3.6
#> stage4     6.38 1.85  12.3
#> 
#> --- 
#> 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.005975181 0.04755877
#> stage2 -0.696719165 1.02336736
#> stage3 -0.013954095 1.41836202
#> stage4  0.798983503 2.57018978
#> attr(,"Probability")
#> [1] 0.95

sm <- summary(fit)
c(DIC = sm$dic, WAIC = sm$waic, LPML = sm$lpml)
#>       DIC      WAIC      LPML 
#>  295.6199 -148.4887 -148.5117

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.0209  
#>  2      1          2     2    0.000523
#>  3      1          3     3    0.00608 
#>  4      1          4     4    0.0316  
#>  5      1          5     5    0.0192  
#>  6      1          6     6    0.0128  
#>  7      1          7     7    0.0406  
#>  8      1          8     8    0.00790 
#>  9      1          9     9    0.00378 
#> 10      1         10    10    0.0296  
#> # ℹ 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.9376593  1
#> 1.1     1  0.6  77     76     1      1          1     1    2 0.8897280  1
#> 1.2     1  0.6  77     76     1      1          1     1    3 0.8471786  1
#> 1.3     1  0.6  77     76     1      1          2     2    1 0.9410663  1
#> 1.4     1  0.6  77     76     1      1          2     2    2 0.8943769  1
#> 1.5     1  0.6  77     76     1      1          2     2    3 0.8527992  1

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

Survival prediction

pr <- predict(fit, times = seq(0, max(larynx$time), length.out = 50))
head(pr)
#>   id      time      surv     lower     upper     cumhaz     std.err
#> 1  1 0.0000000 1.0000000 1.0000000 1.0000000 0.00000000 0.000000000
#> 2  1 0.2183673 0.9730279 0.9560485 0.9930562 0.02738597 0.009080395
#> 3  1 0.4367347 0.9477253 0.9130656 0.9801310 0.05383983 0.016379217
#> 4  1 0.6551020 0.9238187 0.8753110 0.9665140 0.07952982 0.022235645
#> 5  1 0.8734694 0.9010734 0.8390094 0.9477490 0.10461775 0.026938269
#> 6  1 1.0918367 0.8792876 0.8204992 0.9457512 0.12925881 0.030730311

See the Survival prediction and ggplot vignette for ribbon workflows.