Bayes factors for Cox regression

bayesian_survival.R

# approach from
# https://www.medrxiv.org/content/10.1101/2022.11.02.22281762v1.full.pdf

library(tidyverse)
library(survival)
library(ggsurvfit)
#devtools::install_github("maxlinde/baymedr")
library(baymedr)
# library(BayesSurvival)

d <- tribble(~ condition, ~time, ~event,
             "reliable", .1, 0,
             "reliable", .4, 0,
             "reliable", .4, 0,
             "reliable", .6, 0,
             "reliable", .6, 0,
             "reliable", 1, 0,
             "reliable", 1, 0,
             "reliable", 1.5, 0,
             "reliable", 1.5, 0,
             "reliable", 1.5, 0,
             "reliable", 2, 0,
             "reliable", 3, 1,
             "reliable", 3, 1,
             "unreliable", .4, 0,
             "unreliable", .5, 0, 
             "unreliable", .6, 0,
             "unreliable", 1, 0,
             "unreliable", 2, 0,
             "unreliable", 3, 1,
             "unreliable", 3, 1,
             "unreliable", 3, 1,
             "unreliable", 3, 1,
             "unreliable", 3, 1,
             "unreliable", 3, 1,
             "unreliable", 3, 1,
             "unreliable", 3, 1)

d$event <- 1-d$event

Surv(d$time, d$event)

model <- survfit(Surv(time, event)~ condition, data = d)

summary(model)

ggsurvfit(model) +
  add_confidence_interval() + 
  add_risktable()

survdiff(Surv(time, event) ~ condition, data = d)

# bayesian version
df_cox <- d[,c(2,3,1)] |>
  mutate(condition = 1 - (as.numeric(as.factor(condition))-1))

coxph_bf(data = df_cox, 
         null_value = 0, 
         alternative = "one.sided", 
         direction = "high", 
         prior_mean = 0, 
         prior_sd = 1)