Bayesian Cox Models with graph-structure priors

This is a R/Rcpp package BayesSurvive for Bayesian survival models with graph-structured selection priors for sparse identification of high-dimensional features predictive of survival (Madjar et al., 2021) and its extensions with the use of a fixed graph via a Markov Random Field (MRF) prior for capturing known structure of high-dimensional features, e.g. disease-specific pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

Examples

Simulate data

library("BayesSurvive")
# Load the example dataset
data("simData", package = "BayesSurvive")
dataset = list("X" = simData[[1]]$X, 
               "t" = simData[[1]]$time,
               "di" = simData[[1]]$status)

Run a Bayesian Cox model

## Initial value: null model without covariates
initial = list("gamma.ini" = rep(0, ncol(dataset$X)))
# Prior parameters
hyperparPooled = list(
  "c0"     = 2,                      # prior of baseline hazard
  "tau"    = 0.0375,                 # sd (spike) for coefficient prior
  "cb"     = 20,                     # sd (slab) for coefficient prior
  "pi.ga"  = 0.02,                   # prior variable selection probability for standard Cox models
  "a"      = -4,                     # hyperparameter in MRF prior
  "b"      = 0.1,                    # hyperparameter in MRF prior
  "G"      = simData$G               # hyperparameter in MRF prior
)   

## run Bayesian Cox with graph-structured priors
fit <- BayesSurvive(survObj = dataset, model.type = "Pooled", MRF.G = TRUE, 
                    hyperpar = hyperparPooled, initial = initial, nIter = 100)

## show posterior mean of coefficients and 95% credible intervals
library("GGally")
plot(fit) + 
  coord_flip() + 
  theme(axis.text.x = element_text(angle = 90, size = 7))

#plot(fit$output$beta.p[,1], type="l")
#fit$output$beta.margin
#fit$output$gamma.margin
#simData[[1]]$trueB

Plot time-dependent Brier scores

The function BayesSurvive::plotBrier() can show the time-dependent Brier scores based on posterior mean of coefficients or Bayesian model averaging.

plotBrier(fit, , survObj.new = dataset)

The integrated Brier score (IBS) can be obtained by the function BayesSurvive::predict().

predict(fit, survObj.new = dataset)

##                     IBS
## Null model          0.09147208
## Bayesian Cox model  0.03433363

Predict survival probabilities and cumulative hazards

The function BayesSurvive::predict() can estimate the survival probabilities and cumulative hazards.

predict(fit, survObj.new = dataset, type = c("cumhazard", "survival"))

##        observation times cumhazard survival
##              <int> <num>     <num>    <num>
##     1:           1   3.3  2.11e-04 1.00e+00
##     2:           2   3.3  3.29e-01 7.20e-01
##     3:           3   3.3  2.06e-06 1.00e+00
##     4:           4   3.3  1.19e-02 9.88e-01
##     5:           5   3.3  5.36e-04 9.99e-01
##   ---                                     
##  9996:          96   9.5  2.67e+01 2.57e-12
##  9997:          97   9.5  1.08e+03 0.00e+00
##  9998:          98   9.5  2.23e+00 1.08e-01
##  9999:          99   9.5  3.72e+00 2.42e-02
## 10000:         100   9.5  3.37e+01 2.38e-15

Run a ‘Pooled’ Bayesian Cox model with graphical learning

hyperparPooled <- append(hyperparPooled, list("lambda" = 3, "nu0" = 0.05, "nu1" = 5))
fit2 <- BayesSurvive(survObj = list(dataset), model.type = "Pooled", MRF.G = FALSE,
                     hyperpar = hyperparPooled, initial = initial, nIter = 10)

Run a Bayesian Cox model with subgroups using fixed graph

# specify a fixed joint graph between two subgroups
hyperparPooled$G <- Matrix::bdiag(simData$G, simData$G)
dataset2 <- simData[1:2]
dataset2 <- lapply(dataset2, setNames, c("X", "t", "di", "X.unsc", "trueB"))
fit3 <- BayesSurvive(survObj = dataset2, 
                     hyperpar = hyperparPooled, initial = initial, 
                     model.type="CoxBVSSL", MRF.G = TRUE, 
                     nIter = 10, burnin = 5)

Run a Bayesian Cox model with subgroups using graphical learning

fit4 <- BayesSurvive(survObj = dataset2, 
                     hyperpar = hyperparPooled, initial = initial, 
                     model.type="CoxBVSSL", MRF.G = FALSE, 
                     nIter = 3, burnin = 0)