dynsurv (version 0.2-2)

bayesCox: Fit Bayesian Cox Model for Interval Censored Survival Data

Description

Fit Bayesian Cox model with time-independent, time-varying or dynamic covariate coefficient. The fit is done within a Gibbs sampling framework. The reversible jump algorithm is employed for the dynamic coefficient model. The baseline hazards are allowed to be either time-varying or dynamic.

Usage

bayesCox(formula, data, grid, out, model=c("TimeIndep", "TimeVarying", "Dynamic"), base.prior=list(), coef.prior=list(), gibbs=list(), control=list())

Arguments

formula

a formula object, with the response on the left of a '~' operator, and the terms on the right. The response must be a survival object as returned by the Surv function.

data

a data.frame in which to interpret the variables named in the formula.

grid

vector of pre-specified time grid points.

out

name of Markov chain Monte Carlo (MCMC) samples output file.

model

model type to fit.

base.prior

list of options for prior of baseline lambda. Use list(type="Gamma", shape=0.1, rate=0.1) for all models; list(type="Const", value=1) for Dynamic model when intercept=TRUE.

coef.prior

list of options for prior of coefficient beta. Use list(type="Normal", mean=0, sd=1) for TimeIndep model; list(type="AR1", sd=1) for TimeVarying and Dynamic models; list(type="HAR1", shape=2, scale=1) for TimeVarying and Dynamic models.

gibbs

list of options for Gibbs sampler.

control

list of general control options.

Value

An object of S3 class bayesCox representing the fit.

Details

To use default hyper parameters in the specification of either base.prior or coef.prior, one only has to supply the name of the prior, e.g., list(type="Gamma"), list(type="HAR1").

The gibbs argument is a list of components:

iter:

number of iterations, default 3000.

burn:

number of burning, default 500.

thin:

number of thinning, default 1.

verbose:

a logical value, default TRUE. If TRUE, print the iteration.

nReport:

print frequency, default 100.

The control argument is a list of components:

intercept:

a logical value, default FALSE. If TRUE, the model will estimate the intercept, which is the log of baseline hazards. If TRUE, please remember to turn off the direct estimation of baseline hazards, i.e., base.prior=list(type="Const").

a0:

multiplier for initial variance in time-varying or dynamic models, default 100.

eps0:

size of auxiliary uniform latent variable in dynamic model, default 1.

References

X. Wang, M.-H. Chen, and J. Yan (2011). Bayesian dynamic regression models for interval censored survival data. Under review.

D. Sinha, M.-H. Chen, and S.K. Ghosh (1999). Bayesian analysis and model selection for interval-censored survival data. Biometrics 55(2), 585--590.

See Also

coef.bayesCox, jump.bayesCox, nu.bayesCox, plotCoef, plotJumpTrace, and plotNu.

Examples

## Not run: 
# ################################################################################
# # Load one of the following two data sets
# ################################################################################
# 
# # breast cancer data
# data(bcos) ## load bcos and bcos.grid
# mydata <- bcos
# mygrid <- bcos.grid
# myformula <- Surv(left, right, type="interval2") ~ trt
# 
# # tooth data
# # data(tooth) ## load tooth and tooth.grid
# # mydata <- tooth
# # mygrid <- tooth.grid
# # myformula <- Surv(left, right, type="interval2") ~ dmf + sex
# 
# ################################################################################
# # Fit Bayesian Cox models
# ################################################################################
# 
# # Fit time-independent coefficient model
# fit0 <- bayesCox(myformula, mydata, mygrid, out="tiCox.txt",
#                  model="TimeIndep",
#                  base.prior=list(type="Gamma", shape=0.1, rate=0.1),
#                  coef.prior=list(type="Normal", mean=0, sd=1),
#                  gibbs=list(iter=100, burn=20, thin=1, verbose=TRUE, nReport=5))
# plotCoef(coef(fit0))
# 
# # Fit time-varying coefficient model
# fit1 <- bayesCox(myformula, mydata, mygrid, out="tvCox.txt",
#                  model="TimeVarying",
#                  base.prior=list(type="Gamma", shape=0.1, rate=0.1),
#                  coef.prior=list(type="AR1", sd=1),
#                  gibbs=list(iter=100, burn=20, thin=1, verbose=TRUE, nReport=5))
# plotCoef(coef(fit1))
# 
# # Fit dynamic coefficient model with time-varying baseline hazards
# fit2 <- bayesCox(myformula, mydata, mygrid, out="dynCox1.txt",
#                  model="Dynamic",
#                  base.prior=list(type="Gamma", shape=0.1, rate=0.1),
#                  coef.prior=list(type="HAR1", shape=2, scale=1),
#                  gibbs=list(iter=100, burn=20, thin=1, verbose=TRUE, nReport=5))
# plotCoef(coef(fit2))
# plotJumpTrace(jump(fit2))
# plotJumpHist(jump(fit2))
# plotNu(nu(fit2))
# 
# # Plot the coefficient estimates from three models together
# plotCoef(rbind(coef(fit0), coef(fit1), coef(fit2)))
# 
# # Fit dynamic coefficient model with dynamic hazards (in log scales)
# fit3 <- bayesCox(myformula, mydata, mygrid, out="dynCox2.txt",
#                  model="Dynamic",
#                  base.prior=list(type="Const"),
#                  coef.prior=list(type="HAR1", shape=2, scale=1),
#                  gibbs=list(iter=100, burn=20, thin=1, verbose=TRUE, nReport=5),
#                  control=list(intercept=TRUE))
# plotCoef(coef(fit3))
# plotJumpTrace(jump(fit3))
# plotJumpHist(jump(fit3))
# plotNu(nu(fit3))
# ## End(Not run)