Description Usage Arguments Value Files created Author(s) References Examples
A function to estimate a regression model with bivariate (possibly right, left, interval or doublyintervalcensored) data. In the case of doubly interval censoring, different regression models can be specified for the onset and event times.
The error density of the regression model is specified as a mixture of Bayesian Gsplines (normal densities with equidistant means and constant variance matrices). This function performs an MCMC sampling from the posterior distribution of unknown quantities.
For details, see Kom<c3><a1>rek (2006) and Kom<c3><a1>rek and Lesaffre (2006).
We explain first in more detail a model without doubly censoring. Let T[i,l], i=1,..., N, l=1, 2 be event times for ith cluster and the first and the second unit. The following regression model is assumed:
log(T[i,l]) = beta'x[i,l] + epsilon[i,l], i=1,..., N, l=1,2
where beta is unknown regression parameter vector and x[i,l] is a vector of covariates. The bivariate error terms epsilon[i] = (epsilon[i,1], epsilon[i,2])', i=1,..., N are assumed to be i.i.d. with a~bivariate density g[epsilon](e[1], e[2]). This density is expressed as a~mixture of Bayesian Gsplines (normal densities with equidistant means and constant variance matrices). We distinguish two, theoretically equivalent, specifications.
(epsilon[1],\,epsilon[2])' is distributed as sum[j[1]=K[1]][K[1]] sum[j[2]=K[2]][K[2]] w[j[1],j[2]] N(mu[(j[1],j[2])], diag(sigma[1]^2, sigma[2]^2))
where sigma[1]^2, sigma[2]^2 are unknown basis variances and mu[(j[1],j[2])] = (mu[1,j[1]], mu[2,j[2]])' is an~equidistant grid of knots symmetric around the unknown point (gamma[1], gamma[2])' and related to the unknown basis variances through the relationship
mu[1,j[1]] = gamma[1] + j[1]*delta[1]*sigma[1], j[1]=K[1],..., K[1]
mu[2,j[2]] = gamma[2] + j[2]*delta[2]*sigma[2], j[2]=K[2],..., K[2]
where delta[1], delta[2] are fixed constants, e.g. delta[1]=delta[2]=2/3 (which has a~justification of being close to cubic Bsplines).
(epsilon[1],\,epsilon[2])' is distributed as (alpha[1], alpha[2])' + S (V[1], V[2])'
where (alpha[1], alpha[2])' is an unknown intercept term and S is a diagonal matrix with tau[1] and tau[2] on a diagonal, i.e. tau[1], tau[2] are unknown scale parameters. (V[1], V[2])' is then standardized bivariate error term which is distributed according to the bivariate normal mixture, i.e.
(V[1], V[2])' is distributed as sum[j[1]=K[1]][K[1]] sum[j[2]=K[2]][K[2]] w[j[1],j[2]] N(mu[(j[1],j[2])], diag(sigma[1]^2, sigma[2]^2))
where mu[(j[1],j[2])] = (mu[1,j[1]], mu[2,j[2]])' is an~equidistant grid of fixed knots (means), usually symmetric about the fixed point (gamma[1], gamma[2])' = (0, 0)' and sigma[1]^2, sigma[2]^2 are fixed basis variances. Reasonable values for the numbers of grid points K[1] and K[2] are K[1]=K[2]=15 with the distance between the two knots equal to delta=0.3 and for the basis variances sigma[1]^2=sigma[2]^2=0.2^2.
Personally, I found Specification 2 performing better. In the paper Kom<c3><a1>rek and Lesaffre (2006) only Specification 2 is described.
The mixture weights w[j[1],j[2]], j[1]=K[1],..., K[1], j[2]=K[2],..., K[2] are not estimated directly. To avoid the constraints 0 < w[j[1],j[2]] < 1 and sum[j[1]=K[1]][K[1]]sum[j[2]=K[2]][K[2]]w[j[1],j[2]]=1 transformed weights a[j[1],j[2]], j[1]=K[1],..., K[1], j[2]=K[2],..., K[2] related to the original weights by the logistic transformation:
a[j[1],j[2]] = exp(w[j[1],j[2]])/sum[m[1]]sum[m[2]] exp(w[m[1],m[2]])
are estimated instead.
A~Bayesian model is set up for all unknown parameters. For more details I refer to Kom<c3><a1>rek and Lesaffre (2006) and to Kom<c3><a1>rek (2006).
If there are doublycensored data the model of the same type as above can be specified for both the onset time and the timetoevent.
1 2 3 4 5 6 7 8 9 10 11 12  bayesBisurvreg(formula, formula2, data = parent.frame(),
na.action = na.fail, onlyX = FALSE,
nsimul = list(niter = 10, nthin = 1, nburn = 0, nwrite = 10),
prior, prior.beta, 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, 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),
dir = getwd())

formula 
model formula for the regression. In the case of doublycensored data, this is the model formula for the onset time. Data are assumed to be sorted according to subjects and within subjects according to the types of the events that determine the bivariate survival distribution, i.e. the response vector must be t[1,1],t[1,2],t[2,1],t[2,2],t[3,1],t[3,2],...,t[n,1],t[n,2]. The rows of the design matrix with covariates must be sorted analogically. The lefthand side of the formula must be an object created using

formula2 
model formula for the regression of the timetoevent in
the case of doublycensored data. Ignored otherwise. The same remark as
for 
data 
optional data frame in which to interpret the variables occuring in the formulas. 
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 
prior.beta 
prior specification for the regression parameters,
in the case of doubly censored 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
bayesBisurvreg first with its argument 
init 
an~optional list with initial values for the MCMC related
to the model given by
and
K[1] and all values in the second column must be between
K[2] and K[2]. See argument 
mcmc.par 
a~list specifying how some of the Gspline parameters
related 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.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. 
A list of class bayesBisurvreg
containing an information
concerning the initial values and prior choices.
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:
one column labeled iteration
with
indeces of MCMC iterations to which the stored sampled values
correspond.
columns labeled k
, Mean.1
, Mean.2
,
D.1.1
, D.2.1
, D.2.2
, where
k = number of mixture components that had probability numerically higher than zero;
Mean.1 = E(epsilon[i,1]);
Mean.2 = E(epsilon[i,2]);
D.1.1 = var(epsilon[i,1]);
D.2.1 = cov(epsilon[i,1], epsilon[i,2]);
D.2.2 = var(epsilon[i,2]);
all related to the distribution of the error term from the model given by formula
.
in the case of doublycensored data, the same
structure as mixmoment.sim
, however related to the model
given by formula2
.
sampled mixture weights
w[k[1],k[2]] of mixture components that had
probabilities numerically higher than zero. Related to the model
given by formula
.
in the case of doublycensored data, the same
structure as mweight.sim
, however related to the model
given by formula2
.
indeces k[1], k[2],
k[1] in {K[1], ..., K[1]},
k[2] in {K[2], ..., K[2]}
of mixture components that had probabilities numerically higher
than zero. It corresponds to the weights in
mweight.sim
. Related to the model given by formula
.
in the case of doublycensored data, the same
structure as mmean.sim
, however related to the model
given by formula2
.
characteristics of the sampled Gspline
(distribution of
(epsilon[i,1],
epsilon[i,2])') related to the model given by
formula
. This file together with mixmoment.sim
,
mweight.sim
and mmean.sim
can be used to reconstruct
the Gspline in each MCMC iteration.
The file has columns labeled gamma1
,
gamma2
, sigma1
, sigma2
, delta1
,
delta2
, intercept1
, intercept2
,
scale1
, scale2
. The meaning of the values in these
columns is the following:
gamma1 = the middle knot gamma[1] in the first dimension. If ‘Specification’ is 2, this column usually contains zeros;
gamma2 = the middle knot gamma[2] in the second dimension. If ‘Specification’ is 2, this column usually contains zeros;
sigma1 = basis standard deviation sigma[1] of the Gspline in the first dimension. This column contains a~fixed value if ‘Specification’ is 2;
sigma2 = basis standard deviation sigma[2] of the Gspline in the second dimension. This column contains a~fixed value if ‘Specification’ is 2;
delta1 = distance delta[1] between the two knots of the Gspline in the first dimension. This column contains a~fixed value if ‘Specification’ is 2;
delta2 = distance delta[2] between the two knots of the Gspline in the second dimension. This column contains a~fixed value if ‘Specification’ is 2;
intercept1 = the intercept term alpha[1] of the Gspline in the first dimension. If ‘Specification’ is 1, this column usually contains zeros;
intercept2 = the intercept term alpha[2] of the Gspline in the second dimension. If ‘Specification’ is 1, this column usually contains zeros;
scale1 = the scale parameter tau[1] of the Gspline in the first dimension. If ‘Specification’ is 1, this column usually contains ones;
scale2 = the scale parameter tau[2] of the Gspline in the second dimension. ‘Specification’ is 1, this column usually contains ones.
in the case of doublycensored data, the same
structure as gspline.sim
, however related to the model
given by formula2
.
fully created only if store$a = TRUE
. The
file contains the transformed weights
a[k[1],k[2]],
k[1]=K[1],..., K[1],
k[2]=K[2],..., K[2] of all mixture
components, i.e. also of components that had numerically zero
probabilities.
This file is related to the model given by formula
.
fully created only if store$a2 =
TRUE
and in the case of doublycensored data, the same
structure as mlogweight.sim
, however related to the model
given by formula2
.
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 double indeces
(k[1], k[2]), values from 1 to (2*K[1]+1)*(2*K[2]+1) are stored here. Function
vecr2matr
can be used to transform it back to double
indeces.
fully created only if store$r2 =
TRUE
and in the case of doublycensored data, the same
structure as r.sim
, however related to the model
given by formula2
.
either one column labeled lambda
or two
columns labeled lambda1
and lambda2
. These are the
values of the smoothing parameter(s) lambda
(hyperparameters of the prior distribution of the transformed
mixture weights a[k[1],k[2]]). This file is
related to the model given by formula
.
in the case of doublycensored data, the same
structure as lambda.sim
, however related to the model
given by formula2
.
sampled values of the regression parameters
beta related to the model given by
formula
. The columns are labeled according to the
colnames
of the design matrix.
in the case of doublycensored data, the same
structure as beta.sim
, however related to the model
given by formula2
.
fully created only if store$y = TRUE
. It
contains sampled (augmented) logevent times for all observations
in the data set.
fully created only if store$y2 =
TRUE
and in the case of doublycensored data, the same
structure as Y.sim
, however related to the model
given by formula2
.
columns labeled loglik
, penalty
or penalty1
and
penalty2
, logprw
. This file is related to the model
given by formula
. The columns have the following meaning.
loglik = N(log(2*pi) + log(sigma[1]) + log(sigma[2])) 0.5*sum[i=1][N]( (sigma[1]^2*tau[1]^2)^(1) * (y[i,1]  x[i,1]'beta  alpha[1]  tau[1]*mu[1,r[i,1]])^2 + (sigma[2]^2*tau[2]^2)^(1) * (y[i,2]  x[i,2]'beta  alpha[2]  tau[2]*mu[2,r[i,2]])^2 )
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] log(p((y[i,1], y[i,2])  r[i], beta, Gspline));
penalty1: If prior$neighbor.system
= "uniCAR"
:
the penalty term for the first dimension not multiplied by
lambda1
;
penalty2: If prior$neighbor.system
= "uniCAR"
:
the penalty term for the second dimension not multiplied by
lambda2
;
penalty: If prior$neighbor.system
is different from "uniCAR"
:
the penalty term not multiplied by lambda
;
logprw = 2*N*log(sum[k[1]]sum[k[2]] exp(a[k[1],k[2]])) + sum[k[1]]sum[k[2]] N[k[1],k[2]]*a[k[1],k[2]], where N[k[1],k[2]] is the number of residuals assigned intrinsincally to the (k[1], k[2])th mixture component.
In other words, logprw
is equal to the conditional
logdensity
sum[i=1][N] log(p(r[i]  Gspline weights)).
in the case of doublycensored data, the same
structure as lambda.sim
, however related to the model
given by formula2
.
Arno<c5><a1>t Kom<c3><a1>rek arnost.komarek[AT]mff.cuni.cz
Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling. Applied Statistics, 41, 337  348.
Kom<c3><a1>rek, A. (2006). Accelerated Failure Time Models for Multivariate IntervalCensored Data with Flexible Distributional Assumptions. PhD. Thesis, Katholieke Universiteit Leuven, Faculteit Wetenschappen.
Kom<c3><a1>rek, A. and Lesaffre, E. (2006). Bayesian semiparametric accelerated failure time model for paired doubly intervalcensored data. Statistical Modelling, 6, 3–22.
Neal, R. M. (2003). Slice sampling (with Discussion). The Annals of Statistics, 31, 705  767.
1 2 3 4 5 6 7 8 9 10  ## See the description of R commands for
## the population averaged AFT model
## with the Signal Tandmobiel data,
## analysis described in Komarek and Lesaffre (2006),
##
## R commands available in the documentation
## directory of this package as
##  see extandmobPA.R and
## http://www.karlin.mff.cuni.cz/~komarek/software/bayesSurv/extandmobPA.pdf
##

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