# R/bayessurvreg2.R In bayesSurv: Bayesian Survival Regression with Flexible Error and Random Effects Distributions

```################################################
#### AUTHOR:     Arnost Komarek             ####
####             (2005)                     ####
####                                        ####
#### FILE:       bayessurvreg2.R            ####
####                                        ####
#### FUNCTIONS:  bayessurvreg2              ####
################################################

### ======================================
### bayessurvreg2
### ======================================
## Survival regression with normal random effects
## and G-spline error
##
## 29/01/2005: start working on it
## 30/01/2005: more or less finished
##
bayessurvreg2 <- function
(  formula,
random,
formula2,
random2,
data = parent.frame(),
na.action = na.fail,
onlyX = FALSE,
nsimul = list(niter = 10, nthin = 1, nburn = 0, nwrite = 10),
prior,
prior.beta,
prior.b,
init = list(iter = 0),
mcmc.par = list(type.update.a = "slice", k.overrelax.a = 1, k.overrelax.sigma = 1, k.overrelax.scale = 1),
prior2,
prior.beta2,
prior.b2,
init2,
mcmc.par2 = list(type.update.a = "slice", k.overrelax.a = 1, k.overrelax.sigma = 1, k.overrelax.scale = 1),
store = list(a = FALSE, a2 = FALSE, y = FALSE, y2 = FALSE, r = FALSE, r2 = FALSE, b = FALSE, b2 = FALSE),
dir = getwd())
{
thispackage = "bayesSurv"
#thispackage = NULL
store <- bayessurvreg2.checkStore(store)
nsimul <- bayessurvreg.checknsimul(nsimul)

transform = function(t){log(t)}
dtransform = function(t){1/t}
transform2 = function(t){t}
dtransform2 = function(t){1}

## Give a function call to be recorded in a resulting object.
call <- match.call(expand.dots = TRUE)

## Do not transform the response at this moment
doubly <- ifelse(missing(formula2), FALSE, TRUE)
m <- match.call(expand.dots = FALSE)
des <- bayessurvreg.design(m=m, formula=formula, random=random, data=data, transform=transform2, dtransform=transform2)
if (doubly) des2 <- bayessurvreg.design(m=m, formula=formula2, random=random2, data=data, transform=transform2, dtransform=transform2)
else        des2 <- list(nX = 0, n = des\$n, nrandom = 0, randomInt = FALSE)

if (onlyX){
if (doubly) return(list(X=des\$X, X2=des2\$X))
else        return(des\$X)
}

if (!des\$nrandom) store\$b <- FALSE
if (!des2\$nrandom) store\$b2 <- FALSE

## Perform some checks
nobs <- des\$n
if (doubly){
nobs2 <- des2\$n
if (nobs != nobs2) stop("Inconsistent formula and formula2 (different number of observations indicated)")
}

## Priors and inits for beta parameters
if (missing(init))        init <- list()
if (missing(init2))       init2 <- list()

if (missing(prior.beta)) prior.beta <- list()

if (missing(prior.beta2)) prior.beta2 <- list()

## Priors and inits for random effects and their covariance matrix
if (missing(prior.b)) prior.b <- list()
reffdi <- bayessurvreg2.priorb(prior.b, init, des)
init\$b <- attr(reffdi, "init.b")
init\$D <- attr(reffdi, "init.D")
prior.b <- attr(reffdi, "prior.b")

if (missing(prior.b2)) prior.b2 <- list()
reffdi2 <- bayessurvreg2.priorb(prior.b2, init2, des2)
init2\$b <- attr(reffdi2, "init.b")
init2\$D <- attr(reffdi2, "init.D")
prior.b2 <- attr(reffdi2, "prior.b")

## Priors and inits for G-spline, (censored) observations and allocations
if (!doubly){
prior2 <- list()
mcmc.par2 <- list()
}
prinit <- bayessurvreg2.priorInit(prior, init, des, mcmc.par, prior2, init2, des2, mcmc.par2, doubly)
init      <- attr(prinit, "init")
prior     <- attr(prinit, "prior")
mcmc.par  <- attr(prinit, "mcmc.par")
init2     <- attr(prinit, "init2")
prior2    <- attr(prinit, "prior2")
mcmc.par2 <- attr(prinit, "mcmc.par2")

## Compute quantities to determine the space needed to be allocated
##   and numbers of iterations in different phases
if (nsimul\$nburn >= nsimul\$niter) nsimul\$nburn <- nsimul\$niter - 1
if (nsimul\$nburn < 0) nsimul\$nburn <- 0

if (nsimul\$nburn == 0) nruns <- 1
else                   nruns <- 2

nrun <- numeric(2)
nrun[2] <- nsimul\$niter - nsimul\$nburn
nrun[1] <- nsimul\$nburn

nwrite.run <- nrun
nwrite.run[nsimul\$nwrite <= nrun] <- nsimul\$nwrite
max.nwrite <- max(nwrite.run)

## Write headers to files with stored values
bayessurvreg2.writeHeaders(dir, doubly, prinit, store, des, des2)

## Combine similar parameters into one vector
dims <- c(nobs, as.numeric(doubly))
storeV <- c(store\$a, store\$y, store\$r, store\$b, store\$a2, store\$y2, store\$r2, store\$b2)
nsimul.run1 <- c(nrun[1], nsimul\$nthin, nwrite.run[1])
nsimul.run2 <- c(nrun[2], nsimul\$nthin, nwrite.run[2])

cat("Simulation started on                       ", date(), "\n", sep = "")
fit <- .C(C_bayessurvreg2, as.character(dir),
dims = as.integer(dims),
X1 = as.double(if(des\$nX) t(des\$X) else 0),
X2 = as.double(if(des2\$nX) t(des2\$X) else 0),
y1.left = as.double(prinit\$y.left),
y1.right = as.double(prinit\$y.right),
status1 = as.integer(prinit\$status),
t2.left = as.double(prinit\$t2.left),
t2.right = as.double(prinit\$t2.right),
status2 = as.integer(prinit\$status2),
iPML = double(dims[1]),
Ys1 = as.double(prinit\$y),
Ys2 = as.double(prinit\$y2),
r1 = as.integer(prinit\$r),
r2 = as.integer(prinit\$r2),
specif = as.integer(prinit\$specification),
r1.b = integer(1),
r2.b = integer(1),
specif.b = integer(1),
GsplI1 = as.integer(prinit\$Gparmi),
GsplD1 = as.double(prinit\$Gparmd),
GsplI2 = as.integer(prinit\$Gparmi2),
GsplD2 = as.double(prinit\$Gparmd2),
priorbI1 = as.integer(reffdi\$bparmI),
priorbD1 = as.double(reffdi\$bparmD),
priorbI2 = as.integer(reffdi2\$bparmI),
priorbD2 = as.double(reffdi2\$bparmD),
priorCovMatI1 = as.integer(reffdi\$DparmI),
priorCovMatD1 = as.double(reffdi\$DparmD),
priorCovMatI2 = as.integer(reffdi2\$DparmI),
priorCovMatD2 = as.double(reffdi2\$DparmD),
rhob       = as.double(0),
rho.accept = as.integer(0),
rhobI      = as.integer(0),
rhobD      = as.double(0),
mclass.sens.spec = as.double(0),
mclass.logtime   = as.double(0),
mclass.status    = as.integer(0),
mclass.paramI    = as.integer(c(0, 0, 0, rep(0, nobs), 0, 0)),
mclass.paramD    = as.double(c(1, 1, 1, 1)),
iter = as.integer(prinit\$iter),
nsimul = as.integer(nsimul.run1),
store = as.integer(storeV),
version = as.integer(2),
mainSimul = as.integer(0),
err = integer(1),
PACKAGE = thispackage)

if (fit\$err != 0) stop ("Something went wrong during the simulation.")
cat("Burn-up finished on                         ", date(), "   (iteration ", fit\$iter, ")", "\n", sep = "")

## Rewrite sampled values by new files
bayessurvreg2.writeHeaders(dir, doubly, prinit, store, des, des2)

## Main simulation
fit <- .C(C_bayessurvreg2, as.character(dir),
dims = as.integer(dims),
X1 = as.double(fit\$X1),
X2 = as.double(fit\$X2),
y1.left = as.double(fit\$y1.left),
y1.right = as.double(fit\$y1.right),
status1 = as.integer(fit\$status1),
t2.left = as.double(fit\$t2.left),
t2.right = as.double(fit\$t2.right),
status2 = as.integer(fit\$status2),
iPML = double(dims[1]),
Ys1 = as.double(fit\$Ys1),
Ys2 = as.double(fit\$Ys2),
r1 = as.integer(fit\$r1),
r2 = as.integer(fit\$r2),
specif = as.integer(fit\$specif),
r1.b = integer(1),
r2.b = integer(1),
specif.b = integer(1),
GsplI1 = as.integer(fit\$GsplI1),
GsplD1 = as.double(fit\$GsplD1),
GsplI2 = as.integer(fit\$GsplI2),
GsplD2 = as.double(fit\$GsplD2),
priorBetaI1 = as.integer(fit\$priorBetaI1),
priorBetaI2 = as.integer(fit\$priorBetaI2),
priorbI1 = as.integer(fit\$priorbI1),
priorbD1 = as.double(fit\$priorbD1),
priorbI2 = as.integer(fit\$priorbI2),
priorbD2 = as.double(fit\$priorbD2),
priorCovMatI1 = as.integer(fit\$priorCovMatI1),
priorCovMatD1 = as.double(fit\$priorCovMatD1),
priorCovMatI2 = as.integer(fit\$priorCovMatI2),
priorCovMatD2 = as.double(fit\$priorCovMatD2),
rhob       = as.double(0),
rho.accept = as.integer(0),
rhobI      = as.integer(0),
rhobD      = as.double(0),
mclass.sens.spec = as.double(0),
mclass.logtime   = as.double(0),
mclass.status    = as.integer(0),
mclass.paramI    = as.integer(c(0, 0, 0, rep(0, nobs), 0, 0)),
mclass.paramD    = as.double(c(1, 1, 1, 1)),
iter = as.integer(fit\$iter),
nsimul = as.integer(nsimul.run2),
store = as.integer(fit\$store),
version = as.integer(fit\$version),
mainSimul = as.integer(1),
err = integer(1),
PACKAGE = thispackage)

if (fit\$err != 0) stop ("Something went wrong during the simulation.")
cat("Simulation finished on                      ", date(), "   (iteration ", fit\$iter, ")", "\n", sep = "")

toreturn <- fit\$iter
attr(toreturn, "call") <- call

attr(toreturn, "init") <- init
attr(toreturn, "prior") <- prior
attr(toreturn, "prior.beta") <- prior.beta
attr(toreturn, "prior.b") <- prior.b
attr(toreturn, "mcmc.par") <- mcmc.par

if (doubly){
attr(toreturn, "init2") <- init2
attr(toreturn, "prior2") <- prior2
attr(toreturn, "prior.beta2") <- prior.beta2
attr(toreturn, "prior.b2") <- prior.b2
attr(toreturn, "mcmc.par2") <- mcmc.par2
}

class(toreturn) <- "bayessurvreg2"
return(toreturn)
}
```

## Try the bayesSurv package in your browser

Any scripts or data that you put into this service are public.

bayesSurv documentation built on May 2, 2019, 3:26 a.m.