bayesHistogram {bayesSurv}

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,\dots,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,\dots,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)' \sim \sum_{j_1=-K_1}^{K_1}\sum_{j_2=-K_2}^{K_2} w_{j_1,j_2} N_2(\mu_{(j_1,j_2)},\,\mbox{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,\quad j_1=-K_1,\dots,K_1,

\mu_{2,j_2} = \gamma_2 + j_2\delta_2\sigma_2,\quad j_2=-K_2,\dots,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)' \sim (\alpha_1,\,\alpha_2)'+ \bold{S}\,(Y_1,\,Y_2)'

where (\alpha_1,\,\alpha_2)' is an unknown intercept term and \bold{S} \mbox{ is a diagonal matrix with } \tau_1 \mbox{ and }\tau_2 \mbox{ 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)'\sim \sum_{j_1=-K_1}^{K_1}\sum_{j_2=-K_2}^{K_2} w_{j_1,j_2} N_2(\mu_{(j_1,j_2)},\,\mbox{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)

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 Komárek and Lesaffre (2008) and the PhD. thesis Komárek (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,\dots,\gamma_q, where q is dimension of the G-spline and basis standard deviations \sigma_{0,1},\dots,\sigma_{0,q} are estimated. At the same time the G-spline intercepts \alpha_1,\dots,\alpha_q and the G-spline scale parameters s_{1},\dots,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,\dots,\alpha_q and the G-spline scale parameters s_{1},\dots,s_{q} are estimated at the same time positions of the middle knots \gamma_1,\dots,\gamma_q and basis standard deviations \sigma_{0,1},\dots,\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,\dots,q determines that the subscript k_j of the knots runs over -K_j,\dots,0,\dots,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\dots k_q of the knot whose transformed weight a_{k_1\dots k_q} 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, \dots, K_j.

neighbor.system

identification of the neighboring system for the Markov random field prior of transformed mixture weights a_{k_1\,k_2}. 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\Bigl\{-\frac{\lambda}{2}\sum_{k=-K+m}^K\bigl(\Delta^m a_{k}\bigr)^2\Bigr\},

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 \geq 2.

For bivariate smoothing:

p(a) \propto \exp\Bigl\{ -\frac{\lambda_1}{2}\sum_{k_1=-K_1}^{K_1}\sum_{k_2=-K_2+m}^{K_2} \bigl(\Delta_1^m a_{k_1,k_2}\bigr)^2 -\frac{\lambda_2}{2}\sum_{k_2=-K_2}^{K_2}\sum_{k_1=-K_1+m}^{K_1} \bigl(\Delta_2^m a_{k_1,k_2}\bigr)^2 \Bigr\},

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} - 2a_{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 \Bigl\{-\frac{\lambda}{2}\sum_{k_1=-K_1}^{K_1-1}\sum_{k_2=-K_2}^{K_2-1} \bigl(\Delta a_{k_1,k_2} \bigr)^2\Bigr\},

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,\dots,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,\dots,\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,\dots,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,\dots,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,\dots,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). Upper limit of intervals is specified by rate.sigma (see below).

prior.intercept

specification of the prior distribution for the intercept terms \alpha_1,\dots,\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,\dots,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). Upper limit of intervals is specified by rate.scale (see below).

c4delta

values of c_1,\dots,c_q which serve to compute the distance \delta_j between two consecutive knots in each dimension. The knot \mu_{j\,k}, j=1,\dots,q, k=-K_j,\dots,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,\dots,\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,\dots,\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,\dots,\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,\dots,\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},\dots,\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},\dots,\sigma_q^{-2}.

rate.sigma

these are rate parameters for the gamma prior (if prior.sigma = gamma) of basis inverse variances \sigma_1^{-2},\dots,\sigma_q^{-2} or upper limits of the uniform prior (if prior.sigma = sduniform) of basis standard deviations \sigma_1,\dots,\sigma_q.

mean.intercept

these are means for the normal prior distribution of the G-spline intercepts (2nd specification) \alpha_1,\dots,\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,\dots,\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},\dots,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},\dots,scale_q^{-2} or upper limits of the uniform prior (if prior.scale = sduniform) of the G-spline scale scale_1,\dots,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_{k_1}, k_1=-K_1,\dots,K_1 if univariate density is estimated; a_{k_1\,k_2}, k_1=-K_1,\dots,K_1, k_2=-K_2,\dots,K_2, 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,\dots,\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 \boldsymbol{\mu}=(\mu_1, \mu_2)' where

\mu_1 = \mu_{1,\,-3} = \gamma_1 -3\,c_1\sigma_1

and

\mu_2 = \mu_{2,\,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,\dots,K_1, k_2=-K_2,\dots,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 = \mbox{E}(Y_{i,1});

Mean.2 = \mbox{E}(Y_{i,2});

D.1.1 = \mbox{var}(Y_{i,1});

D.2.1 = \mbox{cov}(Y_{i,1},\,Y_{i,2});

D.2.2 = \mbox{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, \dots, K_1\}, k_2 \in\{-K_2, \dots, 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,\dots,K_1, k_2=-K_2,\dots,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)\times (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\Bigl\{\log(2\pi) + \log(\sigma_1) + \log(\sigma_2)\Bigr\}- 0.5\sum_{i=1}^N\Bigl\{ (\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 \Bigr\}

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\Bigl\{p\bigl((y_{i,1},\,y_{i,2})\;\big|\;r_{i},\,\mbox{G-spline}\bigr)\Bigr\};

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\bigl\{\sum_{k_1}\sum_{k_2}a_{k_1,\,k_2}\bigr\} + \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\bigl\{p(r_i\;|\;\mbox{G-spline weights})\bigr\}.

Author(s)

Arnošt Komárek [e-mail]

References

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

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

Komárek, A. and Lesaffre, E. (2008). Bayesian accelerated failure time model with multivariate doubly-interval-censored data and flexible distributional assumptions. Journal of the American Statistical Association, 103, 523 - 533.

Komárek, A. and Lesaffre, E. (2006b). Bayesian semi-parametric accelerated failurew time model for paired doubly interval-censored data. Statistical Modelling, 6, 3 - 22.

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