SemiCompRisks (version 2.0)

BayesSurv: The function to perform Bayesian parametric and semi-parametric regression analysis for univariate time-to-event data.

Description

Univariate regression analysis of time-to-event data is performed using Cox proportional hazards model in a Bayesian framework. The prior of baseline hazard function can be specified by either parametric Weibull or non-parametric mixture of piecewise exponential models (PEM).

Usage

BayesSurv(survData, hyperParams, startValues, mcmcParams, numReps, thin, path, 
burninPerc = 0.5, hz.type = "Weibull", nChain = 1)

Arguments

survData

The data frame containing univariate survival data with covariate matrix from n subjects. It is of dimension $n\times(2+p)$, where $p$ is the number of covariates. The first two columns correspond to $y$, $\delta$.

hyperParams

a vector containing hyperparameter values. For PEM model (hz.type="PEM"), it is of length 4; c($a$, $b$, $\alpha$, $c_{\lambda}$). For Weibull model (hz.type="Weibull"), it is of length 4; c($a$, $b$, $c$, $d$).

startValues

a list containing vectors of starting values of parameters. The length of list must be equal to the number of chains (nChain). For PEM model (hz.type="PEM"), each of list components is of length $p+2(J+1)+3$; c($\beta$, $J$, $\mu

mcmcParams

a vector containing variables required for Metropolis-Hastings-Green (MHG) algorithm. For PEM model (hz.type="PEM"), it is of length $(6+num\_s\_propBI+nTime\_lambda$); c($C$, delPert, num_s_propBI, $J_{max}$, $s_{max}$, nTime_lambda, s_prop

numReps

total number of scans

thin

extent of thinning

path

the name of directory where the results are saved

burninPerc

the proportion of burn-in

hz.type

prior specification for baseline hazard functions : either "Weibull" or "PEM"

nChain

the number of chains

Value

• BayesSurv returns an object of class BayesSurv. names(object$chain1) shows the list of posterior samples of model parameters, the number of acceptance in MH algorithm, etc.

Details

The function BayesSurv implements Bayesian semi-parametric (a mixture of piecewise constant) and parametric (Weibull) models. Let $t_{i}$ denote time to event of interest from subject $i=1,...,n$. The covariates $x_{i}$ are incorporated via Cox proportional hazards model: $$h(t_{i} | x_{i}) = h_{0}(t_{i})\exp(\beta x_{i}), t_{i}>0,$$ where $h_0$ is an unspecified baseline hazard function and $\beta$ is a vector of $p$ log-hazard ratio regression parameters. For "semi-parametric analysis", let $s_{\max}$ denote the largest observed event time. Then, consider the finite partition of the relevant time axis into $J + 1$ disjoint intervals: $0

where $I(\cdot)$ is the indicator function and $s_0 \equiv 0$. In our proposed Bayesian framework, our prior choices are:

$$\pi(\beta) \propto 1,$$ $$\lambda | J, \mu_{\lambda}, \sigma_{\lambda}^2 \sim \mathcal{N}_{J+1}(\mu_{\lambda}1, \sigma_{\lambda}^2\Sigma_{\lambda}),$$ $$J \sim \mathcal{P}(\alpha),$$ $$\pi(s | J) \propto \frac{(2J+1)! \prod_{j=1}^{J+1}(s_j-s_{j-1})}{(s_{J+1})^{(2J+1)}},$$ $$\pi(\mu_{\lambda}) \propto 1,$$ $$\sigma_{\lambda}^{-2} \sim \mathcal{G}(a, b).$$ Note that $J$ and $s$ are treated as random and the priors for $J$ and $s$ jointly form a time-homogeneous Poisson process prior for the partition. The number of time splits and their positions are therefore updated within our computational scheme using reversible jump MCMC. For "parametric analysis", we assume the Weibull model for $h_{0}(t) = \alpha \kappa t^{\alpha-1}$. In our Bayesian framework, our prior choices are: $$\pi(\beta) \propto 1,$$ $$\pi(\alpha) \sim \mathcal{G}(a, b),$$ $$\pi(\kappa) \sim \mathcal{G}(c, d).$$

References

Lee, K. H., Haneuse, S., Schrag, D., and Dominici, F. (2014), Bayesian semiparametric analysis of semicompeting risks data: investigating hospital readmission after a pancreatic cancer diagnosis, Journal of the Royal Statistical Society: Series C, in press.

See Also

print.BayesSurv, summary.BayesSurv, plot.BayesSurv

Examples

# loading a data set	
data(survData)

n = dim(survData)[1]
p = dim(survData)[2]-2

###############################
# setting hyperparameter values

a	<- 0.5	# prior parameters for 1/sigma^2
b	<- 0.01
alpha	<- 5	# prior parameter for J
c_lam 	<- 1	# prior parameter for MVN-ICAR specification

hyperParams <- c(a, b, alpha, c_lam)

#########################
# setting starting values

s_max	<- max(survData$time[survData$event == 1])

beta	<- rep(1, p)
s	<- c(seq(8, max(survData$time[survData$event == 1]), 8), s_max);
J	<- length(s) - 1
lambda	<- runif(J+1, -3, 0)
sigSq_lam<- var(lambda)
mu_lam	<- mean(lambda)

# chain 1

startValues <- list()
startValues[[1]] <- as.vector(c(beta, J, mu_lam, sigSq_lam, lambda, s))

# chain 2

beta	<- rep(0.2, p)
lambda	<- runif(J+1, -3, -1)

startValues[[2]] <- as.vector(c(beta, J, mu_lam, sigSq_lam, lambda, s))

##################################################
# setting variable values needed for MHG algorithm

C		<- 0.70
delPert		<- 0.5
s_propBI	<- floor(sort(unique(survData$time[survData$event == 1])))
s_propBI	<- s_propBI[s_propBI < s_max]
num_s_propBI	<- length(s_propBI)
J_max 		<- 30 					
time_lambda	<- seq(1:39)
nTime_lambda	<- length(time_lambda)

mcmcParams <- c(C, delPert, num_s_propBI, J_max, s_max, nTime_lambda, s_propBI, time_lambda)

##################################################
# number of chains

numReps	= 2e6
thin	= 20
burninPerc	= 0.5
path1	= "outcome/"
hz.type = "PEM"
nChain	= 2

# fitting Bayesian semi-parametric regression model to univariate survival data	

fitSurv <- BayesSurv(survData, hyperParams, startValues, mcmcParams, numReps, thin, path1, 
burninPerc, hz.type, nChain)
		
print(fitSurv)
summary(fitSurv)		
		
## plot for estimates of baseline hazard function
plot(fitSurv)