| Title: | Penalized Semiparametric Bayesian Cox Models |
|---|---|
| Description: | Algorithms to speed up the Bayesian Lasso Cox model (Lee et al., Int J Biostat, 2011 <doi:10.2202/1557-4679.1301>) and the Bayesian Lasso Cox with mandatory variables (Zucknick et al. Biometrical J, 2015 <doi:10.1002/bimj.201400160>). |
| Authors: | Zhi Zhao [aut, cre], Manuela Zucknick [aut], Maral Saadati [aut], Axel Benner [aut] |
| Maintainer: | Zhi Zhao <[email protected]> |
| License: | GPL-3 |
| Version: | 2.0.7 |
| Built: | 2024-10-30 06:49:27 UTC |
| Source: | CRAN |
psbcSpeedUp
Extract the point estimates of the regression coefficients
## S3 method for class 'psbcSpeedUp' coef(object, type = "mean", ...)## S3 method for class 'psbcSpeedUp' coef(object, type = "mean", ...)
object |
an object of class |
type |
type of point estimates of regressions. One of
|
... |
not used |
Estimated coefficients are from an object of class psbcSpeedUp.
If the psbcSpeedUp specified data standardization, the fitted values
are base based on standardized data.
# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1 ) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) coef(fitBayesCox)# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1 ) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) coef(fitBayesCox)
Simulated data set for a quick test. The data set is a list with six
components: survival times "t", event status "di", covariates
"x", number of genomics variables "p", number of clinical
variables "1" and true effects of covariates "beta_true".
The R code for generating the simulated data is given in the Examples
paragraph.
exampleDataexampleData
An object of class list of length 6.
# Load the example dataset data("exampleData", package = "psbcSpeedUp") str(exampleData) # =============== # The code below is to show how to generate the dataset "exampleData.rda" # =============== requireNamespace("MASS", quietly = TRUE) ########################### Predefined Functions Expo <- function(times, surv) { z1 <- -log(surv[1]) t1 <- times[1] lambda <- z1 / (t1) list(rate = lambda) } Weibull <- function(times, surv) { z1 <- -log(surv[1]) z2 <- -log(surv[2]) t1 <- times[1] t2 <- times[2] gamma <- log(z2 / z1) / log(t2 / t1) lambda <- z1 / (t1^gamma) list(scale = lambda, shape = gamma) } ########################### Problem Dimensions n <- 200 p <- 30 q <- 5 s <- 10 ############################ Simulate a set of n x p covariates # effects bg <- c(0.75, -0.75, 0.5, -0.5, 0.25, -0.25, rep(0, p - 6)) bc <- c(-1.0, 1.0, 0.3, 0, -0.3) bX <- c(bg, bc) # covariates # genomic means <- rep(0, p) Sigma <- diag(1, p) Xg <- MASS::mvrnorm(n, means, Sigma) # clinical x1 <- rbinom(n = n, size = 1, prob = 0.7) x2 <- rbinom(n = n, size = 1, prob = 0.3) x3 <- rnorm(n = n, mean = 0, sd = 1) x4 <- rnorm(n = n, mean = 0, sd = 1) x5 <- rnorm(n = n, mean = 0, sd = 1) Xc <- cbind(x1, x2, x3, x4, x5) # all X <- data.frame(Xg, Xc) names(X) <- c(paste("G", 1:p, sep = ""), paste("C", 1:q, sep = "")) X <- scale(X) # censoring function # - follow-up time 36 to 72 months # - administrative censoring: uniform data entry (cens1) # - loss to follow-up: exponential, 20% loss in 72 months (cens2) ACT <- 36 FUT <- 72 cens.start <- FUT cens.end <- ACT + FUT cens1 <- runif(n, cens.start, cens.end) loss <- Expo(times = 72, surv = 0.8) cens2 <- rexp(n, rate = loss$rate) cens <- pmin(cens1, cens2) # survival distribution (Weibull, survival probs 0.5 and 0.9 at 12 and 36 months) h0 <- round(log(2) / 36, 2) surv <- Weibull(times = c(12, 36), surv = c(0.9, 0.5)) dt <- (-log(runif(n)) * (1 / surv$scale) * exp(-as.matrix(X) %*% bX))^(1 / surv$shape) # survival object status <- ifelse(dt <= cens, 1, 0) os <- pmin(dt, cens) exampleData <- list("t" = os, "di" = status, "x" = X, "beta_true" = bX)# Load the example dataset data("exampleData", package = "psbcSpeedUp") str(exampleData) # =============== # The code below is to show how to generate the dataset "exampleData.rda" # =============== requireNamespace("MASS", quietly = TRUE) ########################### Predefined Functions Expo <- function(times, surv) { z1 <- -log(surv[1]) t1 <- times[1] lambda <- z1 / (t1) list(rate = lambda) } Weibull <- function(times, surv) { z1 <- -log(surv[1]) z2 <- -log(surv[2]) t1 <- times[1] t2 <- times[2] gamma <- log(z2 / z1) / log(t2 / t1) lambda <- z1 / (t1^gamma) list(scale = lambda, shape = gamma) } ########################### Problem Dimensions n <- 200 p <- 30 q <- 5 s <- 10 ############################ Simulate a set of n x p covariates # effects bg <- c(0.75, -0.75, 0.5, -0.5, 0.25, -0.25, rep(0, p - 6)) bc <- c(-1.0, 1.0, 0.3, 0, -0.3) bX <- c(bg, bc) # covariates # genomic means <- rep(0, p) Sigma <- diag(1, p) Xg <- MASS::mvrnorm(n, means, Sigma) # clinical x1 <- rbinom(n = n, size = 1, prob = 0.7) x2 <- rbinom(n = n, size = 1, prob = 0.3) x3 <- rnorm(n = n, mean = 0, sd = 1) x4 <- rnorm(n = n, mean = 0, sd = 1) x5 <- rnorm(n = n, mean = 0, sd = 1) Xc <- cbind(x1, x2, x3, x4, x5) # all X <- data.frame(Xg, Xc) names(X) <- c(paste("G", 1:p, sep = ""), paste("C", 1:q, sep = "")) X <- scale(X) # censoring function # - follow-up time 36 to 72 months # - administrative censoring: uniform data entry (cens1) # - loss to follow-up: exponential, 20% loss in 72 months (cens2) ACT <- 36 FUT <- 72 cens.start <- FUT cens.end <- ACT + FUT cens1 <- runif(n, cens.start, cens.end) loss <- Expo(times = 72, surv = 0.8) cens2 <- rexp(n, rate = loss$rate) cens <- pmin(cens1, cens2) # survival distribution (Weibull, survival probs 0.5 and 0.9 at 12 and 36 months) h0 <- round(log(2) / 36, 2) surv <- Weibull(times = c(12, 36), surv = c(0.9, 0.5)) dt <- (-log(runif(n)) * (1 / surv$scale) * exp(-as.matrix(X) %*% bX))^(1 / surv$shape) # survival object status <- ifelse(dt <= cens, 1, 0) os <- pmin(dt, cens) exampleData <- list("t" = os, "di" = status, "x" = X, "beta_true" = bX)
Plot point estimates of regression coefficients and 95% credible intervals
## S3 method for class 'psbcSpeedUp' plot(x, type = "mean", interval = TRUE, ...)## S3 method for class 'psbcSpeedUp' plot(x, type = "mean", interval = TRUE, ...)
x |
an object of class |
type |
type of point estimates of regression coefficients. One of
|
interval |
logical argument to show 95% credible intervals. Default
is |
... |
additional arguments sent to |
ggplot object
# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1 ) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) plot(fitBayesCox, color = "blue")# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1 ) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) plot(fitBayesCox, color = "blue")
Predict time-dependent Brier scores based on Cox regression models
plotBrier(object, survObj.new = NULL, method = "mean", times = NULL, ...)plotBrier(object, survObj.new = NULL, method = "mean", times = NULL, ...)
object |
fitted object obtained with |
survObj.new |
a list containing observed data from new subjects with
components |
method |
option to use the posterior mean ( |
times |
maximum time point to evaluate the prediction |
... |
not used |
psbcSpeedUp
# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) # predict survival probabilities of the train data plotBrier(fitBayesCox, times = 80)# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) # predict survival probabilities of the train data plotBrier(fitBayesCox, times = 80)
Predict survival probability, (cumulative) hazard or (integrated) Brier scores based on Cox regression models
## S3 method for class 'psbcSpeedUp' predict( object, survObj.new = NULL, type = "brier", method = "mean", times = NULL, ... )## S3 method for class 'psbcSpeedUp' predict( object, survObj.new = NULL, type = "brier", method = "mean", times = NULL, ... )
object |
fitted object obtained with |
survObj.new |
a list containing observed data from new subjects with
components |
type |
option to chose for predicting survival probabilities (one of
|
method |
option to use the posterior mean ( |
times |
time points at which to evaluate the risks. If |
... |
not used |
psbcSpeedUp
# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1) # run Bayesian Lasso Cox library("psbcSpeedUp") library("survival") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) # predict survival probabilities of the train data predict(fitBayesCox)# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1) # run Bayesian Lasso Cox library("psbcSpeedUp") library("survival") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) # predict survival probabilities of the train data predict(fitBayesCox)
This a speed-up and extended version of the function psbcGL() in the R package psbcGrouup
psbcSpeedUp( survObj = NULL, p = 0, q = 0, hyperpar = list(), nIter = 1, burnin = 0, thin = 1, rw = FALSE, outFilePath, tmpFolder = "tmp/" )psbcSpeedUp( survObj = NULL, p = 0, q = 0, hyperpar = list(), nIter = 1, burnin = 0, thin = 1, rw = FALSE, outFilePath, tmpFolder = "tmp/" )
survObj |
a list containing observed data from |
p |
number of covariates for variable selection |
q |
number of mandatory covariates |
hyperpar |
a list containing prior parameter values; among
|
nIter |
the number of iterations of the chain |
burnin |
number of iterations to discard at the start of the chain. Default is 0 |
thin |
thinning MCMC intermediate results to be stored |
rw |
when setting to "TRUE", the conventional random walk Metropolis Hastings algorithm is used. Otherwise, the mean and the variance of the proposal density is updated using the jumping rule described in Lee et al. (2011) |
outFilePath |
path to where the output files are to be written |
tmpFolder |
the path to a temporary folder where intermediate data
files are stored (will be erased at the end of the chain). It is specified
relative to |
psbcSpeedUp
t |
a vector of n times to the event |
di |
a vector of n censoring indicators for the event time (1=event occurred, 0=censored) |
x |
covariate matrix, n observations by p variables |
groupInd |
a vector of p group indicator for each variable |
beta.ini |
the starting values for coefficients |
eta0 |
scale parameter of gamma process prior for the cumulative baseline hazard, |
kappa0 |
shape parameter of gamma process prior for the cumulative baseline hazard, |
c0 |
the confidence parameter of gamma process prior for the cumulative baseline hazard, |
r |
the shape parameter of the gamma prior for |
delta |
the rate parameter of the gamma prior for |
lambdaSq |
the starting value for |
sigmaSq |
the starting value for |
tauSq |
the starting values for |
s |
the set of time partitions for specification of the cumulative baseline hazard function |
h |
the starting values for |
beta.prop.var |
the variance of the proposal density for |
beta.clin.var |
the starting value for the variance of |
An object of class psbcSpeedUp is saved as
obj_psbcSpeedUp.rda in the output file, including the following components:
input - a list of all input parameters by the user
output - a list of the all output estimates:
"beta.p" - a matrix with MCMC intermediate estimates of the regression coefficients.
"h.p" - a matrix with MCMC intermediate estimates of the increments in the cumulative baseline hazard in each interval.
"tauSq.p" - a vector MCMC intermediate estimates of the hyperparameter "tauSq".
"sigmaSq.p" - a vector MCMC intermediate estimates of the hyperparameter "sigmaSq".
"lambdaSq.p" - a vector MCMC intermediate estimates of the hyperparameter "lambdaSq".
"accept.rate" - a vector acceptance rates of individual regression coefficients.
call - the matched call.
Lee KH, Chakraborty S, and Sun J (2011). Bayesian Variable Selection in Semiparametric Proportional Hazards Model for High Dimensional Survival Data. The International Journal of Biostatistics, 7(1):1-32.
Zucknick M, Saadati M, and Benner A (2015). Nonidentical twins: Comparison of frequentist and Bayesian lasso for Cox models. Biometrical Journal, 57(6):959–81.
# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1 ) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) plot(fitBayesCox, color = "blue")# Load the example dataset data("exampleData", package = "psbcSpeedUp") p <- exampleData$p q <- exampleData$q survObj <- exampleData[1:3] # Set hyperparameters mypriorPara <- list( "groupInd" = 1:p, "eta0" = 0.02, "kappa0" = 1, "c0" = 2, "r" = 10 / 9, "delta" = 1e-05, "lambdaSq" = 1, "sigmaSq" = runif(1, 0.1, 10), "beta.prop.var" = 1, "beta.clin.var" = 1 ) # run Bayesian Lasso Cox library("psbcSpeedUp") set.seed(123) fitBayesCox <- psbcSpeedUp(survObj, p = p, q = q, hyperpar = mypriorPara, nIter = 10, burnin = 0, outFilePath = tempdir() ) plot(fitBayesCox, color = "blue")