Description Usage Arguments Details Value References See Also Examples
Fit Bayesian Cox model with timeindependent, timevarying 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 timevarying or dynamic.
1 2 3 4 
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 function 
data 
A data.frame in which to interpret the variables named in the

grid 
Vector of prespecified time grid points for model fitting.
It will be automatically set up from data if it is left unspecified
in the function call. By default, it consists of all the unique
finite endpoints (rounded to two significant numbers) of the
censoring intervals after time zero. The 
out 
Name of Markov chain Monte Carlo (MCMC) samples output file. Each row contains one MCMC sample information. The file is needed for those functions further summarizing estimation results in this package. See Section Details for details. 
model 
Model type to fit. Available options are 
base.prior 
List of options for prior of baseline lambda. Use

coef.prior 
List of options for prior of coefficient beta. Use

gibbs 
List of options for Gibbs sampler. 
control 
List of general control options. 
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:
Number of iterations, default 3000;
Number of burning, default 500;
Number of thinning, default 1;
A logical value, default TRUE
. If
TRUE
, print the iteration;
Print frequency, default 100.
The control
argument is a list of components:
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")
Multiplier for initial variance in timevarying or dynamic models, default 100;
Size of auxiliary uniform latent variable in dynamic model, default 1.
For users interested in extracting MCMC sampling information from the output files, the detail of the output files is presented as follows: Let k denote the number of time points (excluding time zero) specified in grid, ck equal 1 for model with timeinvariant coefficients; ck equal k otherwise, and p denote the number of covariates. Then the each sample saved in each row consists of the following possible parts.
The first k numbers represent the jump size of
baseline hazard function at each time grid. If we take the column mean
of the first k columns of the output file, we will get the same
numbers with obj$est$lambda
, where obj
is the bayesCox
object returned by the function.
The sequence from (k + 1) to (k + ck * p) represent the coefficients of covariates at the time grid. The first k numbers in the sequence are the coefficients for the first covariate at the time grid; The second k numbers' subsequence are the coefficients for the second covariate and so on.
The sequence from (k + ck * p + 1) to (k + ck * p + p) represents the sampled latent variance of coefficients.
The sequence from (k + ck * p + p + 1) to (k + 2 * ck * p + p) represents the indicator of whether there is a jump of the covariate coefficients at the time grid. Similar with Part 2, the first k numbers' subsequence is for the first covariate, the second k numbers' subsequence is for the second covariate, and so on.
For the model with timeindependent coefficients, the output file only
has Part 1 and Part 2 in each row; For timevarying coefficient model,
the output file has Part 1, 2, and 3; The output file for the dynamic
model has all the four parts. Note that the dynamic baseline hazard will
be taken as one covariate. So p needs being replaced with
(p + 1) for model with dynamic baseline hazard rate.
No function in the package actually needs the Part 1 from the output file
now; The Part 2 is used by function coef
and survCurve
;
The Part 3 is needed by function nu
; Function jump
extracts
the Part 4.
An object of S3 class bayesCox
representing the fit.
X. Wang, M.H. Chen, and J. Yan (2013). Bayesian dynamic regression models for interval censored survival data with application to children dental health. Lifetime data analysis, 19(3), 297–316.
X. Wang, X. Sinha, J. Yan, and M.H. Chen (2014). Bayesian inference of intervalcensored survival data. In: D. Chen, J. Sun, and K. Peace, Intervalcensored timetoevent data: Methods and applications, 167–195.
X. Wang, M.H. Chen, and J. Yan (2011). Bayesian dynamic regression models for interval censored survival data. Technical Report 13, Department of Statistics, University of Connecticut.
D. Sinha, M.H. Chen, and S.K. Ghosh (1999). Bayesian analysis and model selection for intervalcensored survival data. Biometrics 55(2), 585–590.
coef.bayesCox
, jump.bayesCox
,
nu.bayesCox
, plotCoef
,
plotJumpTrace
, plotNu
,
survCurve
, survDiff
, and
plotSurv
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81  ## Not run:
library(dynsurv)
set.seed(1216)
############################################################################
### Attach one of the following two data sets
############################################################################
### breast cancer data
data(bcos) # attach bcos and bcos.grid
mydata < bcos
mygrid < bcos.grid
myformula < Surv(left, right, type = "interval2") ~ trt
### tooth data
## data(tooth) # attach tooth and tooth.grid
## mydata < tooth
## mygrid < tooth.grid
## myformula < Surv(left, rightInf, type = "interval2") ~ dmf + sex
############################################################################
### Fit Bayesian Cox models
############################################################################
## Fit timeindependent coefficient model
fit0 < bayesCox(myformula, mydata, 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 = 20))
plotCoef(coef(fit0, level = 0.9))
## Plot the estimated survival function for given new data
newDat < data.frame(trt = c("Rad", "RadChem"))
row.names(newDat) < unique(newDat$trt)
plotSurv(survCurve(fit0, newDat), conf.int = TRUE)
## Fit timevarying coefficient model
fit1 < bayesCox(myformula, mydata, out = "tvCox.txt", model = "TimeVary",
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 = 20))
plotCoef(coef(fit1))
plotSurv(survCurve(fit1), conf.int = TRUE)
## Fit dynamic coefficient model with timevarying baseline hazards
fit2 < bayesCox(myformula, mydata, 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 = 20))
plotCoef(coef(fit2))
plotJumpTrace(jump(fit2))
plotJumpHist(jump(fit2))
plotNu(nu(fit2))
plotSurv(survCurve(fit2), conf.int = TRUE)
## 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, 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 = 20),
control = list(intercept = TRUE))
plotCoef(coef(fit3))
plotJumpTrace(jump(fit3))
plotJumpHist(jump(fit3))
plotNu(nu(fit3))
plotSurv(survCurve(fit3), conf.int = TRUE)
## Plot the estimated survival function and the difference
plotSurv(survCurve(fit3, newdata = newDat, type = "survival"),
legendName = "Treatment", conf.int = TRUE)
plotSurv(survDiff(fit3, newdata = newDat, type = "survival"),
legendName = "Treatment", conf.int = TRUE, smooth = TRUE)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.