bayessurvreg2 {bayesSurv}  R Documentation 
Clusterspecific accelerated failure time model for multivariate, possibly doublyintervalcensored data. The error distribution is expressed as a penalized univariate normal mixture with high number of components (Gspline). The distribution of the vector of random effects is multivariate normal.
Description
A function to estimate a regression model with possibly clustered (possibly right, left, interval or doublyinterval censored) data. In the case of doublyinterval censoring, different regression models can be specified for the onset and event times.
(Multivariate) random effects, normally distributed and acting as in the linear mixed model, normally distributed, can be included to adjust for clusters.
The error density of the regression model is specified as a mixture of Bayesian Gsplines (normal densities with equidistant means and constant variances). This function performs an MCMC sampling from the posterior distribution of unknown quantities.
For details, see Komárek (2006), and Komárek, Lesaffre and Legrand (2007).
We explain first in more detail a model without doubly censoring.
Let T_{i,l},\; i=1,\dots, N,\; l=1,\dots, n_i
be event times for i
th cluster and the units within that cluster
The following regression model is assumed:
\log(T_{i,l}) = \beta'x_{i,l} + b_i'z_{i,l} + \varepsilon_{i,l},\quad i=1,\dots, N,\;l=1,\dots, n_i
where \beta
is unknown regression parameter vector,
x_{i,l}
is a vector of covariates.
b_i
is a (multivariate) clusterspecific random effect
vector and z_{i,l}
is a vector of covariates for random
effects.
The random effect vectors b_i,\;i=1,\dots, N
are assumed to be i.i.d. with a (multivariate) normal distribution
with the mean \beta_b
and a covariance matrix
D
. Hierarchical centring (see Gelfand, Sahu, Carlin, 1995) is
used. I.e. \beta_b
expresses the average effect of the
covariates included in z_{i,l}
. Note that covariates
included in z_{i,l}
may not be included in the covariate
vector x_{i,l}
. The covariance matrix D
is
assigned an inverse Wishart prior distribution in the next level of hierarchy.
The error terms
\varepsilon_{i,l},\;i=1,\dots, N, l=1,\dots, n_i
are assumed to be i.i.d. with a univariate density
g_{\varepsilon}(e)
. This density is expressed as
a mixture of Bayesian Gsplines (normal densities with equidistant
means and constant variances). We distinguish two,
theoretically equivalent, specifications.
 Specification 1

\varepsilon \sim \sum_{j=K}^{K} w_{j} N(\mu_{j},\,\sigma^2)
where
\sigma^2
is the unknown basis variance and\mu_{j},\;j=K,\dots, K
is an equidistant grid of knots symmetric around the unknown point\gamma
and related to the unknown basis variance through the relationship\mu_{j} = \gamma + j\delta\sigma,\quad j=K,\dots,K,
where
\delta
is fixed constants, e.g.\delta=2/3
(which has a justification of being close to cubic Bsplines).  Specification 2

\varepsilon \sim \alpha + \tau\,V
where
\alpha
is an unknown intercept term and\tau
is an unknown scale parameter.V
is then standardized error term which is distributed according to the univariate normal mixture, i.e.V\sim \sum_{j=K}^{K} w_{j} N(\mu_{j},\,\sigma^2)
where
\mu_{j},\;j=K,\dots, K
is an equidistant grid of fixed knots (means), usually symmetric about the fixed point\gamma=0
and\sigma^2
is fixed basis variance. Reasonable values for the numbers of grid pointsK
isK=15
with the distance between the two knots equal to\delta=0.3
and for the basis variance\sigma^2=0.2^2.
Personally, I found Specification 2 performing better. In the paper Komárek, Lesaffre and Legrand (2007) only Specification 2 is described.
The mixture weights
w_{j},\;j=K,\dots, K
are
not estimated directly. To avoid the constraints
0 < w_{j} < 1
and
\sum_{j=K}^{K}\,w_j = 1
transformed weights a_{j},\;j=K,\dots, K
related to the original weights by the logistic transformation:
a_{j} = \frac{\exp(w_{j})}{\sum_{m}\exp(w_{m})}
are estimated instead.
A Bayesian model is set up for all unknown parameters. For more details I refer to Komárek (2006) and to Komárek, Lesafre, and Legrand (2007).
If there are doublycensored data the model of the same type as above can be specified for both the onset time and the timetoevent.
Usage
bayessurvreg2(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)
Arguments
formula 
model formula for the regression. In the case of doublycensored data, this is the model formula for the onset time. The lefthand side of the In the formula all covariates appearing both in the vector
If
 
random 
formula for the ‘random’ part of the model, i.e. the
part that specifies the covariates If omitted, no random part is included in the model. E.g. to specify the model with a
random intercept, say When some random effects are included the random intercept is added
by default. It can be removed using e.g.  
formula2 
model formula for the regression of the timetoevent in
the case of doublycensored data. Ignored otherwise. The same structure as
for  
random2 
specification of the ‘random’ part of the model for
timetoevent in the case of doublycensored data. Ignored
otherwise. The same structure as for  
data 
optional data frame in which to interpret the variables
occuring in the  
na.action 
the user is discouraged from changing the default
value  
onlyX 
if  
nsimul 
a list giving the number of iterations of the MCMC and other parameters of the simulation.
 
prior 
a list specifying the prior distribution of the Gspline
defining the distribution of the error term in the regression model
given by The item  
prior.b 
a list defining the way in which the random effects
involved in
 
prior.beta 
prior specification for the regression parameters,
in the case of doublycensored data for the regression parameters of
the onset time, i.e. it is related to This should be a list with the following components:
It is recommended to run the function
bayessurvreg2 first with its argument  
init 
an optional list with initial values for the MCMC related
to the model given by
 
mcmc.par 
a list specifying how some of the Gspline parameters
related to the distribution of the error term from In contrast to  
prior2 
a list specifying the prior distribution of the Gspline
defining the distribution of the error term in the regression model
given by  
prior.b2 
prior specification for the parameters related to the
random effects from  
prior.beta2 
prior specification for the regression parameters
of timetoevent in the case of doubly censored data (related to
 
init2 
an optional list with initial values for the MCMC related
to the model given by  
mcmc.par2 
a list specifying how some of the Gspline parameters
related to  
store 
a list of logical values specifying which chains that are not stored by default are to be stored. The list can have the following components.
 
dir 
a string that specifies a directory where all sampled values are to be stored. 
Value
A list of class bayessurvreg2
containing an information
concerning the initial values and prior choices.
Files created
Additionally, the following files with sampled values
are stored in a directory specified by dir
argument of this
function (some of them are created only on request, see store
parameter of this function).
Headers are written to all files created by default and to files asked
by the user via the argument store
. During the burnin, only
every nsimul$nwrite
value is written. After the burnin, all
sampled values are written in files created by default and to files
asked by the user via the argument store
. In the files for
which the corresponding store
component is FALSE
, every
nsimul$nwrite
value is written during the whole MCMC (this
might be useful to restart the MCMC from some specific point).
The following files are created:
 iteration.sim
one column labeled
iteration
with indeces of MCMC iterations to which the stored sampled values correspond. mixmoment.sim
columns labeled
k
,Mean.1
,D.1.1
, wherek = number of mixture components that had probability numerically higher than zero;
Mean.1 =
\mbox{E}(\varepsilon_{i,l})
;D.1.1 =
\mbox{var}(\varepsilon_{i,l})
;all related to the distribution of the error term from the model given by
formula
. mixmoment_2.sim
in the case of doublycensored data, the same structure as
mixmoment.sim
, however related to the model given byformula2
. mweight.sim
sampled mixture weights
w_{k}
of mixture components that had probabilities numerically higher than zero. Related to the model given byformula
. mweight_2.sim
in the case of doublycensored data, the same structure as
mweight.sim
, however related to the model given byformula2
. mmean.sim
indeces
k,
k \in\{K, \dots, K\}
of mixture components that had probabilities numerically higher than zero. It corresponds to the weights inmweight.sim
. Related to the model given byformula
. mmean_2.sim
in the case of doublycensored data, the same structure as
mmean.sim
, however related to the model given byformula2
. gspline.sim
characteristics of the sampled Gspline (distribution of
\varepsilon_{i,l}
) related to the model given byformula
. This file together withmixmoment.sim
,mweight.sim
andmmean.sim
can be used to reconstruct the Gspline in each MCMC iteration.The file has columns labeled
gamma1
,sigma1
,delta1
,intercept1
,scale1
, The meaning of the values in these columns is the following:gamma1 = the middle knot
\gamma
If ‘Specification’ is 2, this column usually contains zeros;sigma1 = basis standard deviation
\sigma
of the Gspline. This column contains a fixed value if ‘Specification’ is 2;delta1 = distance
delta
between the two knots of the Gspline. This column contains a fixed value if ‘Specification’ is 2;intercept1 = the intercept term
\alpha
of the Gspline. If ‘Specification’ is 1, this column usually contains zeros;scale1 = the scale parameter
\tau
of the Gspline. If ‘Specification’ is 1, this column usually contains ones; gspline_2.sim
in the case of doublycensored data, the same structure as
gspline.sim
, however related to the model given byformula2
. mlogweight.sim
fully created only if
store$a = TRUE
. The file contains the transformed weightsa_{k},
k=K,\dots,K
of all mixture components, i.e. also of components that had numerically zero probabilities. This file is related to the error distribution of the model given byformula
. mlogweight_2.sim
fully created only if
store$a2 = TRUE
and in the case of doublycensored data, the same structure asmlogweight.sim
, however related to the error distribution of the model given byformula2
. r.sim
fully created only if
store$r = TRUE
. The file contains the labels of the mixture components into which the residuals are intrinsically assigned. Instead of indeces on the scale\{K,\dots, K\}
values from 1 to(2\,K+1)
are stored here. Functionvecr2matr
can be used to transform it back to indices fromK
toK
. r_2.sim
fully created only if
store$r2 = TRUE
and in the case of doublycensored data, the same structure asr.sim
, however related to the model given byformula2
. lambda.sim
one column labeled
lambda
. These are the values of the smoothing parameter\lambda
(hyperparameters of the prior distribution of the transformed mixture weightsa_{k}
). This file is related to the model given byformula
. lambda_2.sim
in the case of doublycensored data, the same structure as
lambda.sim
, however related to the model given byformula2
. beta.sim
sampled values of the regression parameters, both the fixed effects
\beta
and means of the random effects\beta_b
(except the random intercept which has always the mean equal to zero). This file is related to the model given byformula
. The columns are labeled according to thecolnames
of the design matrix. beta_2.sim
in the case of doublycensored data, the same structure as
beta.sim
, however related to the model given byformula2
. D.sim
sampled values of the covariance matrix
D
of the random effects. The file has1 + 0.5\,q\,(q+1)
columns (q
is the dimension of the random effect vectorb_i
). The first column labeleddet
contains the determinant of the sampled matrix, additional columns labeledD.1.1
,D.2.1
, ...,D.q.1
, ...D.q.q
contain the lower triangle of the sampled matrix. This file is related to the model specified byformula
andrandom
. D_2.sim
in the case of doublycensored data, the same structure as
D.sim
, however related to the model given byformula2
andrandom2
. b.sim
fully created only if
store$b = TRUE
. It contains sampled values of random effects for all clusters in the data set. The file hasq\times N
columns sorted asb_{1,1},\dots,b_{1,q},\dots, b_{N,1},\dots,b_{N,q}
. This file is related to the model given byformula
andrandom
. b_2.sim
fully created only if
store$b2 = TRUE
and in the case of doublycensored data, the same structure asb.sim
, however related to the model given byformula2
andrandom2
. Y.sim
fully created only if
store$y = TRUE
. It contains sampled (augmented) logevent times for all observations in the data set. Y_2.sim
fully created only if
store$y2 = TRUE
and in the case of doublycensored data, the same structure asY.sim
, however related to the model given byformula2
. logposter.sim
columns labeled
loglik
,penalty
, andlogprw
. This file is related to the model given byformula
. The columns have the following meaning.loglik
=
%  (\sum_{i=1}^N\,n_i)\,\Bigl\{\log(\sqrt{2\pi}) + \log(\sigma) \Bigr\} 0.5\sum_{i=1}^N\sum_{l=1}^{n_i} \Bigl\{ (\sigma^2\,\tau^2)^{1}\; (y_{i,l}  x_{i,l}'\beta  z_{i,l}'b_i  \alpha  \tau\mu_{r_{i,l}})^2 \Bigr\}
where
y_{i,l}
denotes (augmented) (i,l)th true logevent time.In other words,
loglik
is equal to the conditional logdensity\sum_{i=1}^N \sum_{l=1}^{n_i}\,\log\Bigl\{p\bigl(y_{i,l}\;\big\;r_{i,l},\,\beta,\,b_i,\,\mbox{Gspline}\bigr)\Bigr\};
penalty: the penalty term
\frac{1}{2}\sum_{k}\Bigl(\Delta\, a_k\Bigr)^2
(not multiplied by
\lambda
);logprw
=
2\,(\sum_i n_i)\,\log\bigl\{\sum_{k}a_{k}\bigr\} + \sum_{k}N_{k}\,a_{k},
whereN_{k}
is the number of residuals assigned intrinsincally to thek
th mixture component.In other words,
logprw
is equal to the conditional logdensity\sum_{i=1}^N\sum_{l=1}^{n_i} \log\bigl\{p(r_{i,l}\;\;\mbox{Gspline weights})\bigr\}.
 logposter_2.sim
in the case of doublycensored data, the same structure as
logposter.sim
, however related to the model given byformula2
.
Author(s)
Arnošt Komárek arnost.komarek@mff.cuni.cz
References
Gelfand, A. E., Sahu, S. K., and Carlin, B. P. (1995). Efficient parametrisations for normal linear mixed models. Biometrika, 82, 479488.
Komárek, A. (2006). Accelerated Failure Time Models for Multivariate IntervalCensored Data with Flexible Distributional Assumptions. PhD. Thesis, Katholieke Universiteit Leuven, Faculteit Wetenschappen.
Komárek, A., Lesaffre, E., and Legrand, C. (2007). Baseline and treatment effect heterogeneity for survival times between centers using a random effects accelerated failure time model with flexible error distribution. Statistics in Medicine, 26, 54575472.
Examples
## See the description of R commands for
## the model with EORTC data,
## analysis described in Komarek, Lesaffre and Legrand (2007).
##
## R commands available in the documentation
## directory of this package
## as exeortc.R and
## https://www2.karlin.mff.cuni.cz/ komarek/software/bayesSurv/exeortc.pdf
##