bayesHistogram: Smoothing of a uni- or bivariate histogram using Bayesian...

In bayesSurv: Bayesian Survival Regression with Flexible Error and Random Effects Distributions

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 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],...,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 \eqn{j=1,…,q}{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 \eqn{j=1,…,q}{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)

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.

bayesSurv documentation built on Dec. 5, 2022, 5:22 p.m.