Learn R Programming

psbcGroup (version 1.7)

psbcEN: Function to Fit the Penalized Semiparametric Bayesian Cox Model with Elastic Net Prior

Description

Penalized semiparametric Bayesian Cox (PSBC) model with elastic net prior is implemented to analyze survival data with high-dimensional covariates.

Usage

psbcEN(survObj, priorPara, initial, rw=FALSE, mcmcPara, num.reps, 
		thin, chain = 1, save = 1000)

Value

psbcEN returns an object of class psbcEN


beta.p

posterior samples for \(\beta\)

h.p

posterior samples for \(h\)

tauSq.p

posterior samples for \(\tau^2\)

mcmcOutcome

The list containing posterior samples for the remaining model parameters

Arguments

survObj

The list containing observed data from n subjects; t, di, x

priorPara

The list containing prior parameter values; eta0, kappa0, c0, r1, r2, delta1, delta2, s

initial

The list containing the starting values of the parameters; beta.ini, lambda1Sq, lambda2, sigmaSq, tauSq, h

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).

mcmcPara

The list containing the values of options for Metropolis-Hastings step for \(\beta\); numBeta, beta.prop.var

num.reps

the number of iterations of the chain

thin

thinning

chain

the numeric name of chain in the case when running multiple chains.

save

frequency of storing the results in .Rdata file. For example, by setting "save = 1000", the algorithm saves the results every 1000 iterations.

Author

Kyu Ha Lee, Sounak Chakraborty, (Tony) Jianguo Sun

Details

ta vector of n times to the event
dia vector of n censoring indicators for the event time (1=event occurred, 0=censored)
xcovariate matrix, n observations by p variables
eta0scale parameter of gamma process prior for the cumulative baseline hazard, \(eta0 > 0\)
kappa0shape parameter of gamma process prior for the cumulative baseline hazard, \(kappa0 > 0\)
c0the confidence parameter of gamma process prior for the cumulative baseline hazard, \(c0 > 0\)
r1the shape parameter of the gamma prior for \(\lambda_1^2\)
r2the shape parameter of the gamma prior for \(\lambda_2\)
delta1the rate parameter of the gamma prior for \(\lambda_1^2\)
delta2the rate parameter of the gamma prior for \(\lambda_2\)
sthe set of time partitions for specification of the cumulative baseline hazard function
beta.inithe starting values for \(\beta\)
lambda1Sqthe starting value for \(\lambda_1^2\)
lambda2the starting value for \(\lambda_2\)
sigmaSqthe starting value for \(\sigma^2\)
tauSqthe starting values for \(\tau^2\)
hthe starting values for \(h\)
numBetathe number of components in \(\beta\) to be updated at one iteration
beta.prop.varthe variance of the proposal density for \(\beta\) when rw is set to "TRUE"

References

Lee, K. H., Chakraborty, S., and Sun, J. (2011). Bayesian Variable Selection in Semiparametric Proportional Hazards Model for High Dimensional Survival Data. The International Journal of Biostatistics, Volume 7, Issue 1, Pages 1-32.

Lee, K. H., Chakraborty, S., and Sun, J. (2015). Survival Prediction and Variable Selection with Simultaneous Shrinkage and Grouping Priors. Statistical Analysis and Data Mining, Volume 8, Issue 2, pages 114-127.

Examples

Run this code

if (FALSE) {

# generate some survival data
	
	set.seed(204542)
	
	p = 20
	n = 100
	beta.true <- c(rep(4, 10), rep(0, (p-10)))	
	
	CovX<- diag(0.1, p)
	
	survObj 	<- list()
	survObj$x	<- apply(rmvnorm(n, sigma=CovX, method="chol"), 2, scale)
	
	pred <- as.vector(exp(rowSums(scale(survObj$x, center = FALSE, scale = 1/beta.true))))
	
	t 		<- rexp(n, rate = pred)
	cen		<- runif(n, 0, 8)      
	survObj$t 		<- pmin(t, cen)
	survObj$di 		<- as.numeric(t <= cen)

	priorPara 			<- list()
	priorPara$eta0 		<- 1
	priorPara$kappa0 	<- 1
	priorPara$c0 		<- 2
	priorPara$r1		<- 0.1
	priorPara$r2		<- 1
	priorPara$delta1	<- 0.1
	priorPara$delta2	<- 1
	priorPara$s			<- sort(survObj$t[survObj$di == 1])
	priorPara$s			<- c(priorPara$s, 2*max(survObj$t)
	- max(survObj$t[-which(survObj$t==max(survObj$t))]))
	priorPara$J			<- length(priorPara$s)

	mcmcPara				<- list()
	mcmcPara$numBeta		<- p
	mcmcPara$beta.prop.var	<- 1

	initial				<- list()
	initial$beta.ini	<- rep(0.5, p)
	initial$lambda1Sq	<- 1  
	initial$lambda2		<- 1  
	initial$sigmaSq		<- runif(1, 0.1, 10)
	initial$tauSq		<- rexp(p, rate = initial$lambda1Sq/2)
	initial$h			<- rgamma(priorPara$J, 1, 1)

	rw = FALSE
	num.reps = 20000
	chain = 1
	thin = 5
	save = 5

	fitEN <- psbcEN(survObj, priorPara, initial, rw=FALSE, mcmcPara, 
				num.reps, thin, chain, save)

	vs <- VS(fitEN, X=survObj$x)
    
	}		

Run the code above in your browser using DataLab