bayessurvreg3 {bayesSurv} | R Documentation

## Cluster-specific accelerated failure time model for multivariate, possibly doubly-interval-censored data. The random intercept can be included in the model formula. Both the error distribution and the distribution of the random intercept is expressed as a~penalized univariate normal mixture with high number of components (G-spline).

### Description

A function to estimate a regression model with possibly clustered (possibly right, left, interval or doubly-interval censored) data. In the case of doubly-interval censoring, different regression models can be specified for the onset and event times.

A~univariate random effect (random intercept) with the distribution expressed as a~penalized normal mixture can be included in the model to adjust for clusters.

The error density of the regression model is specified as a mixture of Bayesian G-splines (normal densities with equidistant means and constant variances). This function performs an MCMC sampling from the posterior distribution of unknown quantities.

For details, see Komarek (2006) and Komarek and Lesaffre (2006).

We explain first in more detail a model without doubly censoring. Let T[i,l], i=1,..., N, l=1,..., n[i] be event times for ith cluster and the units within that cluster The following regression model is assumed:

log(T[i,l]) = beta'x[i,l] + b[i] + epsilon[i,l], i=1,..., N, l=1,..., n[i]

where beta is unknown regression parameter vector, x[i,l] is a vector of covariates. b[i] is a cluster-specific random effect (random intercept).

The random effects b[i], i=1,..., N are assumed to be i.i.d. with a~univariate density g[b](b). The error terms epsilon[i,l], i=1,..., N, l=1,..., n[i] are assumed to be i.i.d. with a~univariate density g[epsilon](e).

Densities g[b] and g[epsilon] are both expressed as a~mixture of Bayesian G-splines (normal densities with equidistant means and constant variances). We distinguish two, theoretically equivalent, specifications.

In the following, the density for epsilon is explicitely described. The density for b is obtained in an analogous manner.

Specification 1

  epsilon is distributed as sum[j=-K][K] w[j] N(mu[j], sigma^2)

  where sigma^2 is the unknown basis variance and mu[j],;j=-K,..., 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, j=K,..., K

  where delta is fixed constants, e.g. delta=2/3 (which has a~justification of being close to cubic B-splines).

  {}

Specification 2

  epsilon[1] is distributed as 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 is distributed as sum[j=-K][K] w[j] N(mu[j], sigma^2)

  where mu[j], j=-K,..., 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 points K is K=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 Komarek and Lesaffre (2006) only Specification 2 is described.

The mixture weights w[j], j=-K,..., 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,..., K related to the original weights by the logistic transformation:

a[j] = 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 Komarek and Lesaffre (2006) (manuscript can be found in the documentation of this package) and to Komarek (2006).

If there are doubly-censored data the model of the same type as above can be specified for both the onset time and the time-to-event.

### Usage

    bayessurvreg3(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,
                       type.update.a.b = "slice", k.overrelax.a.b = 1,
                       k.overrelax.sigma.b = 1, k.overrelax.scale.b = 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,
                       type.update.a.b = "slice", k.overrelax.a.b = 1,
                       k.overrelax.sigma.b = 1, k.overrelax.scale.b = 1),
       store = list(a = FALSE, a2 = FALSE, y = FALSE, y2 = FALSE,
                    r = FALSE, r2 = FALSE, b = FALSE, b2 = FALSE,
                    a.b = FALSE, a.b2 = FALSE, r.b = FALSE, r.b2 = FALSE),
       dir = getwd())

### Arguments

`formula` | model formula for the regression. In the case of doubly-censored data, this is the model formula for the onset time. The left-hand side of the `formula` must be an~object created using `Surv`. Intercept is implicitely included in the model by the estimation of the error distribution. As a~consequence `-1` in the model formula does not have any effect on the model specification. If `random` is used then the formula must contain an identification of clusters in the form `cluster(id)`, where `id` is a name of the variable that determines clusters, e.g.
`Surv(time, event)~gender + cluster(id)`.
`random` | formula for the `random' part of the model. In the case of doubly-censored data, this is the `random` formula for the onset time. With this version of the function only
`random = ~1`
is allowed. If omitted, no random part is included in the model.
`formula2` | model formula for the regression of the time-to-event in the case of doubly-censored data. Ignored otherwise. The same structure as for `formula` applies here.
`random2` | specification of the `random' part of the model for time-to-event in the case of doubly-censored data. Ignored otherwise. The same structure as for `random` applies here.
`data` | optional data frame in which to interpret the variables occuring in the `formula`, `formula2`, `random`, `random2` statements.
`na.action` | the user is discouraged from changing the default value `na.fail`.
`onlyX` | if `TRUE` no MCMC sampling is performed and only the design matrix (matrices) are returned. This can be useful to set up correctly priors for regression parameters in the presence of `factor` covariates.
`nsimul` | a list giving the number of iterations of the MCMC and other parameters of the simulation. niter total number of sampled values after discarding thinned ones, burn-up included; nthin thinning interval; nburn number of sampled values in a burn-up period after discarding thinned values. This value should be smaller than `niter`. If not, `nburn` is set to `niter - 1`. It can be set to zero; nwrite an interval at which information about the number of performed iterations is print on the screen and during the burn-up period an interval with which the sampled values are writen to files;
`prior` | a~list specifying the prior distribution of the G-spline defining the distribution of the error term in the regression model given by `formula` and `random`. See `prior` argument of `bayesHistogram` function for more detail. In this list also ‘Specification’ as described above is specified. The item `prior$neighbor.system` can only be equal to `uniCAR` here.
`prior.b` | a~list specifying the prior distribution of the G-spline defining the distribution of the random intercept in the regression model given by `formula` and `random`. See `prior` argument of `bayesHistogram` function for more detail. In this list also ‘Specification’ as described above is specified. The item `prior.b$neighbor.system` can only be equal to `uniCAR` here.
`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 `formula` and `random`. This should be a~list with the following components: mean.prior a~vector specifying a~prior mean for each `beta` parameter in the model. var.prior a~vector specifying a~prior variance for each `beta` parameter. It is recommended to run the function bayessurvreg3 first with its argument `onlyX` set to `TRUE` to find out how the betas are sorted. They must correspond to a design matrix X taken from `formula`.
`init` | an~optional list with initial values for the MCMC related to the model given by `formula` and `random`. The list can have the following components: iter the number of the iteration to which the initial values correspond, usually zero. beta a~vector of initial values for the regression parameters. It must be sorted in the same way as are the columns in the design matrix. Use `onlyX=TRUE` if you do not know how the columns in the design matrix are created. a a~vector of length 2*K+1 with the initial values of transformed mixture weights for the G-spline defining the distribution of the error term epsilon. lambda initial values for the Markov random fields precision parameter for the G-spline defining the distribution of the error term epsilon. gamma an~initial values for the middle knot gamma for the G-spline defining the distribution of the error term epsilon. If ‘Specification’ is 2, this value will not be changed by the MCMC and it is recommended (for easier interpretation of the results) to set `init$gamma` to zero (default behavior). If ‘Specification’ is 1 `init$gamma` should be approximately equal to the mean value of the residuals. sigma an~initial values of the basis standard deviation sigma for the G-spline defining the distribution of the error term epsilon. If ‘Specification’ is 2, this value will not be changed by the MCMC and it is recommended to set it approximately equal to the range of standardized data (let say 4 + 4) divided by the number of knots and multiplied by something like 2/3. If ‘Specification’ is 1 this should be approximately equal to the range of the residuals divided by the number of knots (2*K+1) and multiplied again by something like 2/3. intercept an~initial values of the intercept term alpha for the G-spline defining the distribution of the error term epsilon. If ‘Specification’ is 1 this value is not changed by the MCMC and the initial value is always changed to zero. scale an~initial value of the scale parameter tau for the G-spline defining the distribution of the error term epsilon. If ‘Specification’ is 1 this value is not changed by the MCMC and the initial value is always changed to one. a.b a~vector of length 2*K+1 with the initial values of transformed mixture weights for the G-spline defining the distribution of the random intercept b. lambda.b initial values for the Markov random fields precision parameter for the G-spline defining the distribution of the random intercept b. gamma.b an~initial values for the middle knot gamma for the G-spline defining the distribution of the random intercept b. Due to identifiability reasons, this value is always changed to zero and is for neither ‘Specification’ updated by the MCMC. sigma.b an~initial values of the basis standard deviation sigma for the G-spline defining the distribution of the random intercept b. If ‘Specification’ is 2, this value will not be changed by the MCMC and it is recommended to set it approximately equal to the range of standardized data (let say 4 + 4) divided by the number of knots and multiplied by something like 2/3. If ‘Specification’ is 1 this should be approximately equal to the range of the residuals divided by the number of knots (2*K+1) and multiplied again by something like 2/3. intercept.b an~initial values of the intercept term alpha for the G-spline defining the distribution of the random intercept b. Due to identifiability reasons, this value is always changed to zero and is for neither ‘Specification’ updated by the MCMC. scale.b an~initial value of the scale parameter tau for the G-spline defining the distribution of the random intercept b. If ‘Specification’ is 1 this value is not changed by the MCMC and the initial value is always changed to one. b a~vector of length N of the initial values of random effects b[i],;i=1,..., N for each cluster. y a~vector of length sum[i=1][N] n[i] with initial values of log-event-times. r a~vector of length sum[i=1][N] n[i] with initial component labels for each residual. All values must be between -K and K. See argument `init` of the function `bayesHistogram` for more details. r.b a~vector of length N with initial component labels for each random intercept. All values must be between -K and K. See argument `init` of the function `bayesHistogram` for more details.
`mcmc.par` | a~list specifying how some of the G-spline parameters related to the distribution of the error term and of the random intercept from `formula` and `random` are to be updated. See `bayesBisurvreg` for more details. Compared to the mcmc.par argument of the function `bayesBisurvreg` additional components related to the G-spline for the random intercept can be present, namely
`type.update.a.b`
`k.overrelax.a.b`
`k.overrelax.sigma.b`
`k.overrelax.scale.b`
`prior2` | a~list specifying the prior distribution of the G-spline defining the distribution of the error term in the regression model given by `formula2` and `random2`. See `prior` argument of `bayesHistogram` function for more detail.
`prior.b2` | prior specification for the parameters related to the random effects from `formula2` and `random2`. This should be a~list with the same structure as `prior.b`.
`prior.beta2` | prior specification for the regression parameters of time-to-event in the case of doubly censored data (related to `formula2` and `random2`). This should be a~list with the same structure as `prior.beta`.
`init2` | an~optional list with initial values for the MCMC related to the model given by `formula2` and `random2`. The list has the same structure as `init`.
`mcmc.par2` | a~list specifying how some of the G-spline parameters related to `formula2` and `random2` are to be updated. The list has the same structure as `mcmc.par`.
`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. a if `TRUE` then all the transformed mixture weights a[k], k=-K,..., K, related to the G-spline defining the error distribution of `formula` are stored. a.b if `TRUE` then all the transformed mixture weights a[k], k=-K,..., K, related to the G-spline defining the distribution of the random intercept from `formula` and `random` are stored. a2 if `TRUE` and there are doubly-censored data then all the transformed mixture weights a[k], k=-K,..., K, related to the G-spline defining the error distribution of `formula2` are stored. a.b2 if `TRUE` then all the transformed mixture weights a[k], k=-K,..., K, related to the G-spline defining the distribution of the random intercept from `formula2` and `random2` are stored. y if `TRUE` then augmented log-event times for all observations related to the `formula` are stored. y2 if `TRUE` then augmented log-event times for all observations related to `formula2` are stored. r if `TRUE` then labels of mixture components for residuals related to `formula` are stored. r.b if `TRUE` then labels of mixture components for random intercepts related to `formula` and `random` are stored. r2 if `TRUE` then labels of mixture components for residuals related to `formula2` are stored. r.b2 if `TRUE` then labels of mixture components for random intercepts related to `formula2` and `random2` are stored. b if `TRUE` then the sampled values of the random interceptss related to `formula` and `random` are stored. b2 if `TRUE` then the sampled values of the random interceptss related to `formula2` and `random2` are stored.
`dir` | a string that specifies a directory where all sampled values are to be stored.

### Value

A list of class `bayessurvreg3` 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 burn-in, only every `nsimul$nwrite` value is written. After the burn-in, 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
  this file is related to the density of the error term from the model given by `formula`.

  Columns labeled `k`, `Mean.1`, `D.1.1`, where

  k = number of mixture components that had probability numerically higher than zero;

  Mean.1 = E(epsilon[i,l]);

  D.1.1 = var(epsilon[i,l]).

mixmoment_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  The same structure as `mixmoment.sim`.

mixmoment_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  The same structure as `mixmoment.sim`.

mixmoment_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  The same structure as `mixmoment.sim`.

mweight.sim
  this file is related to the density of the error term from the model given by `formula`.

  Sampled mixture weights w[k] of mixture components that had probabilities numerically higher than zero.

mweight_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  The same structure as `mweight.sim`.

mweight_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  The same structure as `mweight.sim`.

mweight_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  The same structure as `mweight.sim`.

mmean.sim
  this file is related to the density of the error term from the model given by `formula`.

  Indeces k, k in {-K, ..., K} of mixture components that had probabilities numerically higher than zero. It corresponds to the weights in `mweight.sim`.

mmean_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  The same structure as `mmean.sim`.

mmean_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  The same structure as `mmean.sim`.

mmean_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  The same structure as `mmean.sim`.

gspline.sim
  this file is related to the density of the error term from the model given by `formula`.

  Characteristics of the sampled G-spline. This file together with `mixmoment.sim`, `mweight.sim` and `mmean.sim` can be used to reconstruct the G-spline 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 G-spline. This column contains a~fixed value if ‘Specification’ is 2;

  delta1 = distance delta between the two knots of the G-spline. This column contains a~fixed value if ‘Specification’ is 2;

  intercept1 = the intercept term alpha of the G-spline. If ‘Specification’ is 1, this column usually contains zeros;

  scale1 = the scale parameter tau of the G-spline. If ‘Specification’ is 1, this column usually contains ones;

gspline_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  The same structure as `gspline.sim`.

gspline_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  The same structure as `gspline.sim`.

gspline_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  The same structure as `gspline.sim`.

mlogweight.sim
  this file is related to the density of the error term from the model given by `formula`.

  Fully created only if `store$a = TRUE`. The file contains the transformed weights a[k], k=-K,..., K of all mixture components, i.e. also of components that had numerically zero probabilities.

mlogweight_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  Fully created only if `store$a.b = TRUE`.

  The same structure as `mlogweight.sim`.

mlogweight_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  Fully created only if `store$a2 = TRUE`.

  The same structure as `mlogweight.sim`.

mlogweight_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  Fully created only if `store$a.b2 = TRUE`.

  The same structure as `mlogweight.sim`.

r.sim
  this file is related to the density of the error term from the model given by `formula`.

  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,..., K} values from 1 to (2*K+1) are stored here. Function `vecr2matr` can be used to transform it back to indices from -K to K.

r_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  Fully created only if `store$r.b = TRUE`.

  The same structure as `r.sim`.

r_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  Fully created only if `store$r2 = TRUE`.

  The same structure as `r.sim`.

r_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  Fully created only if `store$r.b2 = TRUE`.

  The same structure as `r.sim`.

lambda.sim
  this file is related to the density of the error term from the model given by `formula`.

  One column labeled `lambda`. These are the values of the smoothing parameterlambda (hyperparameters of the prior distribution of the transformed mixture weights a[k]).

lambda_{}b.sim
  this file is related to the density of the random intercept from the model given by `formula` and `random`.

  The same structure as `lambda.sim`.

lambda_{}2.sim
  in the case of doubly-censored data. This file is related to the density of the error term from the model given by `formula2`.

  The same structure as `lambda.sim`.

lambda_{}b2.sim
  in the case of doubly-censored data. This file is related to the density of the random intercept from the model given by `formula2` and `random2`.

  The same structure as `lambda.sim`.

beta.sim
  this file is related to the model given by `formula`.

  Sampled values of the regression parameters beta.

  The columns are labeled according to the `colnames` of the design matrix.

beta_{}2.sim
  in the case of doubly-censored data, the same structure as `beta.sim`, however related to the model given by `formula2`.
b.sim
  this file is related to the model given by `formula` and `random`.

  Fully created only if `store$b = TRUE`. It contains sampled values of random intercepts for all clusters in the data set. The file has N columns.

b_{}2.sim
  fully created only if `store$b2 = TRUE` and in the case of doubly-censored data, the same structure as `b.sim`, however related to the model given by `formula2` and `random2`.
Y.sim
  this file is related to the model given by `formula`.

  Fully created only if `store$y = TRUE`. It contains sampled (augmented) log-event times for all observations in the data set.

Y_{}2.sim
  fully created only if `store$y2 = TRUE` and in the case of doubly-censored data, the same structure as `Y.sim`, however related to the model given by `formula2`.
logposter.sim
  This file is related to the residuals of the model given by `formula`.

  Columns labeled `loglik`, `penalty`, and `logprw`. The columns have the following meaning.

  loglik = -(sum[i=1][N] n[i]) * (log(sqrt(2*pi)) + log(sigma)) -0.5*sum[i=1][N] sum[l=1][n[i]]( (sigma^2*tau^2)^(-1) * (y[i,l] - x[i,l]'beta - b[i] - alpha - tau*mu[r[i,l]])^2)

  where y[i,l] denotes (augmented) (i,l)th true log-event time.

  In other words, `loglik` is equal to the conditional log-density

  sum[i=1][N] sum[l=1][n[i]] log(p(y[i,l] | r[i,l], beta, b[i], error-G-spline));

  penalty: the penalty term

  -0.5*sum[k] (Delta a[k])^2

  (not multiplied by lambda);

  logprw = -2*(sum[i] n[i])*log(sum[k] exp(a[k])) + sum[k[1]] N[k]*a[k], where N[k] is the number of residuals assigned intrinsincally to the kth mixture component.

  In other words, `logprw` is equal to the conditional log-density

  sum[i=1][N] sum[l=1][n[i]] log(p(r[i,l] | error-G-spline weights)).

logposter_{}b.sim
  This file is related to the random intercepts from the model given by `formula` and `random`.

  Columns labeled `loglik`, `penalty`, and `logprw`. The columns have the following meaning.

  loglik = -N * (log(sqrt(2*pi)) + log(sigma)) -0.5*sum[i=1][N]( (sigma^2*tau^2)^(-1) * (b[i] - alpha - tau*mu[r[i]])^2)

  where b[i] denotes (augmented) ith random intercept.

  In other words, `loglik` is equal to the conditional log-density

  sum[i=1][N] log(p(b[i] | r[i], b-G-spline));

  The columns `penalty` and `logprw` have the analogous meaning as in the case of logposter.sim file.

logposter_{}2.sim
  in the case of doubly-censored data, the same structure as `logposter.sim`, however related to the model given by `formula2`.
logposter_{}b2.sim
  in the case of doubly-censored data, the same structure as `logposter\_{}b.sim`, however related to the model given by `formula2` and `random2`.

### Author(s)

Arnost Komarek arnost.komarek@med.kuleuven.be

### References

Komarek, A. (2006). Accelerated Failure Time Models for Multivariate Interval-Censored Data with Flexible Distributional Assumptions. PhD. Thesis, Katholieke Universiteit Leuven, Faculteit Wetenschappen.

Komarek, A. and Lesaffre, E. (2006). Bayesian accelerated failure time model with multivariate doubly-interval-censored data and flexible distributional assumptions. Submitted.
See `Komarek_Lesaffre_2006.pdf`.

### Examples

    ## See the description of R commands for
    ## the cluster specific 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 tandmobCS.pdf, tandmobCS.R.

[Package bayesSurv version 0.4 Index]