bayesCox {dynsurv} | R Documentation |
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 = NULL,
out = NULL,
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 function |
data |
A data.frame in which to interpret the variables named in the
|
grid |
Vector of pre-specified 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 |
An optional character string specifying the name of Markov chain
Monte Carlo (MCMC) samples output file. By default, the MCMC samples
will be output to a temporary directory set by |
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. |
... |
Other arguments that are for futher extension. |
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
. IfTRUE
, print the iteration;- nReport:
Print frequency, default 100.
The control
argument is a list of components:
- intercept:
A logical value, default
FALSE
. IfTRUE
, the model will estimate the intercept, which is the log of baseline hazards. IfTRUE
, 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.
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 time-invariant 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.
- Part 1:
The first
k
numbers represent the jump size of baseline hazard function at each time grid. If we take the column mean of the firstk
columns of the output file, we will get the same numbers withobj$est$lambda
, whereobj
is thebayesCox
object returned by the function.- Part 2:
The sequence from
(k + 1) to (k + ck * p)
represent the coefficients of covariates at the time grid. The firstk
numbers in the sequence are the coefficients for the first covariate at the time grid; The secondk
numbers' sub-sequence are the coefficients for the second covariate and so on.- Part 3:
The sequence from
(k + ck * p + 1)
to(k + ck * p + p)
represents the sampled latent variance of coefficients.- Part 4:
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' sub-sequence is for the first covariate, the secondk
numbers' sub-sequence is for the second covariate, and so on.
For the model with time-independent coefficients, the output file only
has Part 1 and Part 2 in each row; For time-varying 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.
Value
An object of S3 class bayesCox
representing the fit.
References
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 interval-censored survival data. In: D. Chen, J. Sun, and K. Peace, Interval-censored time-to-event 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 interval-censored survival data. Biometrics 55(2), 585–590.
See Also
coef.bayesCox
, jump.bayesCox
,
nu.bayesCox
, plotCoef
,
plotJumpTrace
, plotNu
,
survCurve
, survDiff
, and
plotSurv
.
Examples
## Not run:
library(dynsurv)
set.seed(1216)
## breast cancer data
data(bcos)
mydata <- bcos
myformula <- Surv(left, right, type = "interval2") ~ trt
### Fit time-independent coefficient model
fit0 <- bayesCox(myformula, mydata, 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 = FALSE))
## plot coefficient estimates
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 time-varying coefficient model
fit1 <- bayesCox(myformula, mydata, 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 time-varying baseline hazards
fit2 <- bayesCox(myformula, mydata, 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, 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)
## extract MCMC samples
mcmc_list <- bayesCoxMcmc(fit3, part = "coef")
posterior_coef <- mcmc_list$coef
## posterior probabilities of hazard ratio of RadChem (vs. Rad)
## greater than 1 at time 10
posterior_coef[covariate == "trtRadChem" & time == 10, mean(exp(coef) > 1)]
## End(Not run)