bayesHistogram {bayesSurv} | R Documentation

## Smoothing of a uni- or bivariate histogram using Bayesian G-splines

### Description

A function to estimate a density of a uni- or bivariate (possibly censored) sample. The density is specified as a mixture of Bayesian G-splines (normal densities with equidistant means and equal variances). This function performs an MCMC sampling from the posterior distribution of unknown quantities in the density specification. Other method functions are available to visualize resulting density estimate.

This function served as a basis for further developed `bayesBisurvreg`, `bayessurvreg2` and `bayessurvreg3` functions. However, in contrast to these functions, `bayesHistogram` does not allow for doubly censoring.

Bivariate case:

Let Y[i,l], i=1,..., N, l=1,2 be observations for the ith cluster and the first and the second unit (dimension). The bivariate observations Y[i] = (Y[i,1], Y[i,2])', i=1,..., N are assumed to be i.i.d. with a~bivariate density g[y](y[1], y[2]). This density is expressed as a~mixture of Bayesian G-splines (normal densities with equidistant means and constant variance matrices). We distinguish two, theoretically equivalent, specifications.

Specification 1

  (Y[1],,Y[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 B-splines).

Specification 2

  (Y[1],,Y[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 observational vector 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.

Univariate case:

It is a~direct simplification of the bivariate case.

### Usage

    bayesHistogram(y1, y2,
       nsimul = list(niter = 10, nthin = 1, nburn = 0, nwrite = 10),
       prior, init = list(iter = 0),
       mcmc.par = list(type.update.a = "slice", k.overrelax.a = 1,
                       k.overrelax.sigma = 1, k.overrelax.scale = 1),
       store = list(a = FALSE, y = FALSE, r = FALSE),
       dir = getwd())

### Arguments

`y1` | response for the first dimension in the form of a survival object created using `Surv`.
`y2` | response for the second dimension in the form of a survival object created using `Surv`. If the response is one-dimensional this item is missing.
`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 that identifies prior hyperparameters and prior choices. See the paper Komarek and Lesaffre (2006) and the PhD. thesis Komarek (2006) for more details. Some prior parameters can be guessed by the function itself. If you want to do so, set such parameters to `NULL`. Set to `NULL` also the parameters that are not needed in your model. specification a~number giving which specification of the model is used. It can be one of the following numbers: 1 with this specification positions of the middle knots gamma[1],...,gamma[q], where q is dimension of the G-spline and basis standard deviations sigma[0,1],...,sigma[0,q] are estimated. At the same time the G-spline intercepts alpha[1],...,alpha[q] and the G-spline scale parameters s[1],...,s[q] are assumed to be fixed (usually, intercepts to zero and scales to 1). The user can specified the fixed quantities in the `init` parameter of this function 2 with this specification, G-spline intercepts alpha[1],...,alpha[q] and the G-spline scale parameters s[1],...,s[q] are estimated at the same time positions of the middle knots gamma[1],...,gamma[q] and basis standard deviations sigma[0,1],...,sigma[0,q] are assumed to be fixed (usually, middle knots to zero ans basis standard deviations to some smaller number like 0.2) The user can specified the fixed quantities in the `init` parameter of this function K specification of the number of knots in each dimension, i.e. `K` is a vector of length equal to the dimension of the data q and K[j], j=1,...,q determines that the subscript k[j] of the knots runs over -K[j],...,0,...,K[j]. A value K[j]=0 is valid as well. There are only some restriction on the minimal value of K[j] with respect to the choice of the neighbor system and possibly the order of the conditional autoregression in the prior of transformed weights (see below). izero subscript k[1]...k[q] of the knot whose transformed weight a[k1...kq] will constantly be equal to zero. This is here for identifiability. To avoid numerical problems it is highly recommended to set `izero=rep(0, q)`. `izero[j]` should be taken from the set -K[j],...,K[j]. neighbor.system identification of the neighboring system for the Markov random field prior of transformed mixture weights a[k1, k2]. This can be substring of one of the following strings: `uniCAR` ``univariate conditional autoregression'': a~prior based on squared differences of given order m (see argument `order`) in each row and column. For univariate smoothing: p(a) propto exp(-lambda/2 * sum[k=-K+m][K] (Delta^m a[k])^2), where Delta^m denotes the difference operator of order m, i.e. Delta^1 a[k] = a[k] - a[k-1] and Delta^m a[k] = Delta^(m-1)a[k] - Delta^(m-1)a[k-1], m >= 2. For bivariate smoothing: p(a) propto exp( -lambda[1]/2 * sum[k[1]=-K[1]][K[1]]sum[k[2]=-K[2]+m][K[2]] (Delta[1]^m a[k[1], k[2]])^2 -lambda[2]/2 * sum[k[2]=-K[2]][K[2]]sum[k[1]=-K[1]+m][K[1]] (Delta[2]^m a[k[1], k[2]])^2), where Delta[l]^m denotes the difference operator of order m acting in the lth margin, e.g. Delta[1]^2 = a[k[1], k[2]] - 2*a[k[1], k[2]-1] + a[k[1], k[2]-2]. The precision parameters lambda[1] and lambda[2] might be forced to be equal (see argument `equal.lambda`.) `eight.neighbors` this prior is based on eight nearest neighbors (i.e. except on edges, each full conditional depends only on eight nearest neighbors) and local quadratic smoothing. It applies only in the case of bivariate smoothing. The prior is then defined as p(a) propto exp(-lambda/2 * sum[k[1]=-K[1]][K[1]-1]sum[k[2]=-K[2]][K[2]-1] (Delta a[k[1],k[2]])^2), where Delta a[k[1],k[2]] = a[k[1],k[2]] - a[k[1]+1, k[2]] - a[k[1], k[2]+1] + a[k[1]+1, k[2]+1]. `twelve.neighbors` !!! THIS FEATURE HAS NOT BEEN IMPLEMENTED YET. !!! order order of the conditional autoregression if `neighbor.system = uniCAR`. Implemented are 1, 2, 3. If `order = 0` and `neighbor.system = uniCAR` then mixture weights are assumed to be fixed and equal to their initial values specified by the `init` parameter (see below). Note that the numbers K[j], j=1,...,q must be all equal to or higher than `order`. equal.lambda `TRUE/FALSE` applicable in the case when a density of bivariate observations is estimated and `neighbor.system = uniCAR`. It specifies whether there is only one common Markov random field precision parameter lambda for all margins (dimensions) or whether each margin (dimension) has its own precision parameter lambda. For all other neighbor systems is `equal.lambda` automatically `TRUE`. prior.lambda specification of the prior distributions for the Markov random field precision parameter(s) lambda (when `equal.lambda = TRUE`) or lambda[1],...,lambda[q] (when `equal.lambda = TRUE`). This is a vector of substring of one of the following strings (one substring for each margin if `equal.lambda = FALSE`, otherwise just one substring): `fixed` the lambda parameter is then assumed to be fixed and equal to its initial values given by `init` (see below). `gamma` a particular lambda parameter has a priori gamma distribution with shape g[j] and rate (inverse scale) h[j] where j=1 if `equal.lambda=TRUE` and j=1,...,q if `equal.lambda=TRUE`. Shape and rate parameters are specified by `shape.lambda`, `rate.lambda` (see below). `sduniform` a particular 1/sqrt(lambda) parameter (i.e.a standard deviation of the Markov random field) has a priori a uniform distribution on the interval (0, S[j]) where j=1 if `equal.lambda=TRUE` and j=1,...,q if `equal.lambda=TRUE`. Upper limit of intervals is specified by `rate.lambda` (see below). prior.gamma specification of the prior distribution for a reference knot (intercept) gamma in each dimension. This is a vector of substrings of one of the following strings (one substring for each margin): `fixed` the gamma parameter is then assumed to be fixed and equal to its initial values given by `init` (see below). `normal` the gamma parameter has a priori a normal distribution with mean and variance given by `mean.gamma` and `var.gamma`. prior.sigma specification of the prior distribution for basis standard deviations of the G-spline in each dimension. This is a vector of substrings of one of the following strings (one substring for each margin): `fixed` the sigma parameter is then assumed to be fixed and equal to its initial values given by `init` (see below). `gamma` a particular sigma^{-2} parameter has a priori gamma distribution with shape zeta[j] and rate (inverse scale) eta[j] where j=1,...,q. Shape and rate parameters are specified by `shape.sigma`, `rate.sigma` (see below). `sduniform` a particular sigma parameter has a priori a uniform distribution on the interval (0, S[j]) where j=1,...,q. Upper limit of intervals is specified by `rate.sigma` (see below). prior.intercept specification of the prior distribution for the intercept terms alpha[1],...,alpha[q] (2nd specification) in each dimension. This is a vector of substrings of one of the following strings (one substring for each margin): `fixed` the intercept parameter is then assumed to be fixed and equal to its initial values given by `init` (see below). `normal` the intercept parameter has a priori a normal distribution with mean and variance given by `mean.intercept` and `var.intercept`. prior.scale specification of the prior distribution for the scale parameter (2nd specification) of the G-spline in each dimension This is a vector of substrings of one of the following strings (one substring for each margin): `fixed` the `scale` parameter is then assumed to be fixed and equal to its initial values given by `init` (see below). `gamma` a particular scale^{-2} parameter has a priori gamma distribution with shape zeta[j] and rate (inverse scale) eta[j] where j=1,...,q. Shape and rate parameters are specified by `shape.scale`, `rate.scale` (see below). `sduniform` a particular scale parameter has a priori a uniform distribution on the interval (0, S[j]) where j=1,...,q. Upper limit of intervals is specified by `rate.scale` (see below). c4delta values of c[1],...,c[q] which serve to compute the distance delta[j] between two consecutive knots in each dimension. The knot mu[j,k], j=1,...,q, k=-K[j],...,K[j] is defined as mu[j,k] = gamma[j] + k*delta[j] with delta[j] = c[j]*sigma[j]. mean.gamma these are means for the normal prior distribution of middle knots gamma[1],...,gamma[q] in each dimension if this prior is normal. For fixed gamma an appropriate element of the vector `mean.gamma` may be whatever. var.gamma these are variances for the normal prior distribution of middle knots gamma[1],...,gamma[q] in each dimension if this prior is normal. For fixed gamma an appropriate element of the vector `var.gamma` may be whatever. shape.lambda these are shape parameters for the gamma prior (if used) of Markov random field precision parameters lambda[1],...,lambda[q] (if `equal.lambda = FALSE`) or lambda[1] (if `equal.lambda = TRUE`). rate.lambda these are rate parameters for the gamma prior (if `prior.lambda = gamma`) of Markov random field precision parameters lambda[1],...,lambda[q] (if `equal.lambda = FALSE`) or lambda[1] (if `equal.lambda = TRUE`) or upper limits of the uniform prior (if `prior.lambda = sduniform`) of Markov random field standard deviation parameters lambda[1]^{-1/2},...,lambda[q]^{-1/2} (if `equal.lambda = FALSE`) or lambda[1]^{-1/2} (if `equal.lambda = TRUE`). shape.sigma these are shape parameters for the gamma prior (if used) of basis inverse variances sigma[1]^{-2},...,sigma[q]^{-2}. rate.sigma these are rate parameters for the gamma prior (if `prior.sigma = gamma`) of basis inverse variances sigma[1]^{-2},...,sigma[q]^{-2} or upper limits of the uniform prior (if `prior.sigma = sduniform`) of basis standard deviations sigma[1],...,sigma[q]. mean.intercept these are means for the normal prior distribution of the G-spline intercepts (2nd specification) alpha[1],...,alpha[q] in each dimension if this prior is normal. For fixed alpha an appropriate element of the vector `mean.intercept` may be whatever. var.intercept these are variances for the normal prior distribution of the G-spline intercepts alpha[1],...,alpha[q] in each dimension if this prior is normal. For fixed alpha an appropriate element of the vector `var.alpha` may be whatever. shape.scale these are shape parameters for the gamma prior (if used) of the G-spline scale parameter (2nd specification) scale[1]^{-2},...,scale[q]^{-2}. rate.scale these are rate parameters for the gamma prior (if `prior.scale = gamma`) of the G-spline inverse variances scale[1]^{-2},...,scale[q]^{-2} or upper limits of the uniform prior (if `prior.scale = sduniform`) of the G-spline scale scale[1],...,scale[q].
`init` | a list of the initial values to start the McMC. Set to `NULL` such parameters that you want the program should itself sample for you or parameters that are not needed in your model. iter the number of the iteration to which the initial values correspond, usually zero. a vector/matrix of initial transformed mixture weights a[k1], k1=-K1,...,K1 if univariate density is estimated; a[k1,k2], k1=-K1,...,K1, k2=-K2,...,K2, if bivariate density is estimated. This initial value can be guessed by the function itself. lambda initial values for Markov random field precision parameter(s) lambda (if `equal.lambda = TRUE`), lambda[1],...,lambda[q] (if `equal.lambda = FALSE`.) gamma initial values for the middle knots in each dimension. If `prior$specification = 2` it is recommended (for easier interpretation of the results) to set `init$gamma` to zero for all dimensions. If `prior$specification = 1` `init$gamma` should be approximately equal to the mean value of the data in each margin. sigma initial values for basis standard deviations in each dimension. If `prior$specification = 2` this should be approximately equal to the range of standardized data (let say 4 + 4) divided by the number of knots in each margin and multiplied by something like 2/3. If `prior$specification = 1` this should be approximately equal to the range of your data divided by the number of knots in each margin and multiplied again by something like 2/3. intercept initial values for the intercept term in each dimension. Note that if `prior$specification = 1` this initial value is always changed to zero for all dimensions. scale initial values for the G-spline scale parameter in each dimension. Note that if `prior$specification = 1` this initial value is always changed to one for all dimensions. y initial values for (possibly unobserved censored) observations. This should be either a vector of length equal to the sample size if the response is univariate or a matrix with as many rows as is the sample size and two columns if the response is bivariate. Be aware that `init$y` must be consistent with data supplied. This initial can be guessed by the function itself. Possible missing values in `init$y` tells the function to guess the initial value. r initial values for labels of components to which the (augmented) observations belong. This initial can be guessed by the function itself. This should be either a vector of length equal to the sample size if the response is univariate or a matrix with as many rows as is the sample size and two columns if the response is bivariate. Values in the first column of this matrix should be between `-prior$K[1]` and `prior$K[1]`, values in the second column of this matrix between `-prior$K[2]` and `prior$K[2]`, e.g. when `init$r[i,1:2] = c(-3, 6)` it means that the ith observation is initially assigned to the component with the mean mu = (mu[1], mu[2])' where mu[1] = mu[1,-3] = gamma[1] -3*c[1]*sigma[1] and mu[2] = mu[1,6] = gamma[2] +6*c[2]*sigma[2].
`mcmc.par` | a list specifying further details of the McMC simulation. There are default values implemented for all components of this list. type.update.a it specifies the McMC method to update transformed mixture weights a. It is a~substring of one of the following strings: slice slice sampler of Neal (2003) is used (default choice); ars.quantile adaptive rejection sampling of Gilks and Wild (1992) is used with starting abscissae equal to 15%, 50% and 85% quantiles of a~piecewise exponential approximation to the full conditional from the previous iteration; ars.mode adaptive rejection sampling of Gilks and Wild (1992) is used with starting abscissae equal to the mode and plus/minus twice approximate standard deviation of the full conditional distribution k.overrelax.a this specifies a frequency of overrelaxed updates of transformed mixture weights a when slice sampler is used. Every kth value is sampled in a usual way (without overrelaxation). If you do not want overrelaxation at all, set `k.overrelax.a` to 1 (default choice). Note that overrelaxation can be only done with the slice sampler (and not with adaptive rejection sampling). k.overrelax.sigma a vector of length equal to the dimension of the G-spline specifying a frequency of overrelaxed updates of basis G-spline variances. If you do not want overrelaxation at all, set all components of `k.overrelax.sigma` to 1 (default choice). k.overrelax.scale a vector of length equal to the dimension of the G-spline specifying a frequency of overrelaxed updates of the G-spline scale parameters (2nd specification). If you do not want overrelaxation at all, set all components of `k.overrelax.scale` to 1 (default choice).
`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[1],k[2]], k[1]=-K[1],..., K[1], k[2]=-K[2],..., K[2], related to the G-spline are stored. y if `TRUE` then augmented log-event times for all observations are stored. r if `TRUE` then labels of mixture components for residuals are stored.
`dir` | a string that specifies a directory where all sampled values are to be stored.

### Value

A list of class `bayesHistogram` 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`. 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`, `Mean.2`, `D.1.1`, `D.2.1`, `D.2.2` in the bivariate case and columns labeled `k`, `Mean.1`, `D.1.1` in the univariate case, where

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

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

  Mean.2 = E(Y[i,2]);

  D.1.1 = var(Y[i,1]);

  D.2.1 = cov(Y[i,1], Y[i,2]);

  D.2.2 = var(Y[i,2]).

mweight.sim
  sampled mixture weights w[k[1],k[2]] of mixture components that had probabilities numerically higher than zero.
mmean.sim
  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`.
gspline.sim
  characteristics of the sampled G-spline (distribution of (Y[i,1], Y[i,2])'). 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`, `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 G-spline in the first dimension. This column contains a~fixed value if ‘Specification’ is 2;

  sigma2 = basis standard deviation sigma[2] of the G-spline in the second dimension. This column contains a~fixed value if ‘Specification’ is 2;

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

  delta2 = distance delta[2] between the two knots of the G-spline in the second dimension. This column contains a~fixed value if ‘Specification’ is 2;

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

  intercept2 = the intercept term alpha[2] of the G-spline in the second dimension. If ‘Specification’ is 1, this column usually contains zeros;

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

  scale2 = the scale parameter tau[2] of the G-spline in the second dimension. ‘Specification’ is 1, this column usually contains ones.

mlogweight.sim
  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.
r.sim
  fully created only if `store$r = TRUE`. The file contains the labels of the mixture components into which the observations 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.
lambda.sim
  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]]).
Y.sim
  fully created only if `store$y = TRUE`. It contains sampled (augmented) log-event times for all observations in the data set.
logposter.sim
  columns labeled `loglik`, `penalty` or `penalty1` and `penalty2`, `logprw`. The columns have the following meaning (the formulas apply for the bivariate case).

  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] - alpha[1] - tau[1]*mu[1,r[i,1]])^2 + (sigma[2]^2*tau[2]^2)^(-1) * (y[i,2] - alpha[2] - tau[2]*mu[2,r[i,2]])^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] log(p((y[i,1], y[i,2]) | r[i], G-spline));

  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 observations assigned intrinsincally to the (k[1], k[2])th mixture component.

  In other words, `logprw` is equal to the conditional log-density sum[i=1][N] log(p(r[i] | G-spline weights)).

### Author(s)

Arnost Komarek arnost.komarek@med.kuleuven.be

### References

Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling. Applied Statistics, 41, 337 - 348.

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`.

Komarek, A. and Lesaffre, E. (2006b). Bayesian accelerated failure time model for paired doubly-interval-censored data. Submitted.
See `Komarek_Lesaffre_2006b.pdf`.

Neal, R. M. (2003). Slice sampling (with Discussion). The Annals of Statistics, 31, 705 - 767.

[Package bayesSurv version 0.4 Index]