bayesBisurvreg: Population-averaged accelerated failure time model for bivariate, possibly doubly-interval-censored data. The error distribution is expressed as a~penalized bivariate normal mixture with high number of components (bivariate G-spline).

Description

A function to estimate a regression model with bivariate (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. The error density of the regression model is specified as a mixture of Bayesian G-splines (normal densities with equidistant means and constant variance matrices). This function performs an MCMC sampling from the posterior distribution of unknown quantities.

For details, see Komárek (2006) and Komárek and Lesaffre (2006). We explain first in more detail a model without doubly censoring. Let $T[i,l], i=1,..., N, l=1, 2$ be event times for $i$th cluster and the first and the second unit. The following regression model is assumed: $$\log(T_{i,l}) = \beta'x_{i,l} + \varepsilon_{i,l},\quad i=1,\dots,N,\;l=1,2$$ where $beta$ is unknown regression parameter vector and $x[i,l]$ is a vector of covariates. The bivariate error terms $epsilon[i] = (epsilon[i,1], epsilon[i,2])', i=1,..., N$ are assumed to be i.i.d. with a~bivariate density $g[epsilon](e[1], e[2])$. This density is expressed as a~mixture of Bayesian G-splines (normal densities with equidistant means and constant variance matrices). We distinguish two, theoretically equivalent, specifications.

Specification 1: $$(\varepsilon_1,\,\varepsilon_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

$$(\varepsilon_1,\,\varepsilon_2)' \sim (\alpha_1,\,\alpha_2)'+ \bold{S}\,(V_1,\,V_2)'$$ where $(alpha[1], alpha[2])'$ is an unknown intercept term and $ S is a diagonal matrix with tau[1] and tau[2] on a diagonal,$ i.e. $tau[1], tau[2]$ are unknown scale parameters. $(V[1], V[2])'$ is then standardized bivariate error term which is distributed according to the bivariate normal mixture, i.e. $$(V_1,\,V_2)'\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.$

Personally, I found Specification 2 performing better. In the paper Komárek and Lesaffre (2006) only Specification 2 is described.

The mixture weights $w[j[1],j[2]], j[1]=-K[1],..., K[1], j[2]=-K[2],..., K[2]$ are not estimated directly. To avoid the constraints $0 < w[j[1],j[2]] < 1$ and $sum[j[1]=-K[1]][K[1]]sum[j[2]=-K[2]][K[2]]w[j[1],j[2]]=1$ transformed weights $a[j[1],j[2]], j[1]=-K[1],..., K[1], j[2]=-K[2],..., K[2]$ related to the original weights by the logistic transformation: $$a_{j_1,j_2} = \frac{\exp(w_{j_1,j_2})}{\sum_{m_1}\sum_{m_2}\exp(w_{m_1,m_2})}$$ are estimated instead.

A~Bayesian model is set up for all unknown parameters. For more details I refer to Komárek and Lesaffre (2006) and to Komárek (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

bayesBisurvreg(formula, formula2, data = parent.frame(), na.action = na.fail, onlyX = FALSE, nsimul = list(niter = 10, nthin = 1, nburn = 0, nwrite = 10), prior, prior.beta, init = list(iter = 0), mcmc.par = list(type.update.a = "slice", k.overrelax.a = 1, k.overrelax.sigma = 1, k.overrelax.scale = 1), prior2, prior.beta2, init2, mcmc.par2 = list(type.update.a = "slice", k.overrelax.a = 1, k.overrelax.sigma = 1, k.overrelax.scale = 1), store = list(a = FALSE, a2 = FALSE, y = FALSE, y2 = FALSE, r = FALSE, r2 = FALSE), dir = getwd())

Arguments

formula

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

The left-hand side of the formula must be an object created using Surv.

formula2

model formula for the regression of the time-to-event in the case of doubly-censored data. Ignored otherwise. The same remark as for formula concerning the sort order applies here.

data

optional data frame in which to interpret the variables occuring in the formulas.

na.action

the user is discouraged from changing the default value 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. See prior argument of bayesHistogram function for more detail. In this list also ‘Specification’ as described above is specified.

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. 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 bayesBisurvreg 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. 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~matrix of size $(2*K[1]+1) x (2*K[2]+1)$ with the initial values of transformed mixture weights.
lambda: initial values for the Markov random fields precision parameters. According to the chosen prior for the transformed mixture weights, this is either a~number or a~vector of length 2.
gamma: a~vector of length 2 of initial values for the middle knots $gamma[1], gamma[2]$ in each dimension. 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 for all dimensions (default behavior). If ‘Specification’ is 1 init$gamma should be approximately equal to the mean value of the residuals in each margin.
sigma: a~vector of length~2 of initial values of the basis standard deviations $sigma[1], sigma[2]$. 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 in each margin and multiplied by something like 2/3.
intercept: a~vector of length~2 of initial values of the intercept terms $alpha[1], alpha[2]$. If ‘Specification’ is 1 this value is not changed by the MCMC and the initial value is always changed to zero for both dimensions.
scale: a~vector of length~2 of initial values of the scale parameters $tau[1], tau[2]$.
y: a~matrix with 2 columns and $N$ rows with initial values of log-event-times for each cluster in rows.
r: a~matrix with 2 columns and $N$ rows with initial component labels for each bivariate residual in rows. All values in the first column must be between $-K[1]$

mcmc.par

a~list specifying how some of the G-spline parameters related to formula are to be updated. The list can have the following components (all of them have their default values):

type.update.a

G-spline transformed weights $a$ can be updated using one of the following algorithms:

slice: slice sampler of Neal (2003)
ars.quantile: adaptive rejection sampling of Gilks and Wild (1992) with starting abscissae being quantiles of the envelop at the previous iteration
ars.mode: adaptive rejection sampling of Gilks and Wild (1992) with starting abscissae being the mode plus/minus 3 times estimated standard deviation of the full conditional distribution

Default is slice.

k.overrelax.a

if type.update.a == "slice" some updates are overrelaxed. Then every k.overrelax.ath iteration is not overrelaxed. Default is

k.overrelax.a =
	  1

, i.e. no overrelaxation

k.overrelax.sigma

G-spline basis standard deviations are updated using the slice sampler of Neal (2003). At the same time, overrelaxation can be used. Then every k.overrelax.sigma th update is not overrelaxed. Default is k.overrelax.sigma = 1, i.e. no overrelaxation

k.overrelax.scale

G-spline scales are updated using the slice sampler of Neal (2003). At the same time, overrelaxation can be used. Then every k.overrelax.scale th update is not overrelaxed. Default is k.overrelax.scale = 1, i.e. no overrelaxation

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. See prior argument of bayesHistogram function for more detail.

prior.beta2

prior specification for the regression parameters of time-to-event in the case of doubly censored data (related to formula2). 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. The list has the same structure as init.

mcmc.par2

a~list specifying how some of the G-spline parameters related to formula2 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[1],k[2]],$ $k[1]=-K[1],..., K[1],$ $k[2]=-K[2],..., K[2],$ related to the G-spline of formula are stored.
a2: if TRUE and there are doubly-censored data 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 of formula2 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.
r2: if TRUE then labels of mixture components for residuals related to formula2 are stored.

dir

a string that specifies a directory where all sampled values are to be stored.

Value

A list of class bayesBisurvreg 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: columns labeled k, Mean.1, Mean.2, D.1.1, D.2.1, D.2.2, where k = number of mixture components that had probability numerically higher than zero; Mean.1 = $E(epsilon[i,1])$; Mean.2 = $E(epsilon[i,2])$; D.1.1 = $var(epsilon[i,1])$; D.2.1 = $cov(epsilon[i,1], epsilon[i,2])$; D.2.2 = $var(epsilon[i,2])$; all related to the distribution of the error term from the model given by formula.
mixmoment\_2.sim: in the case of doubly-censored data, the same structure as mixmoment.sim, however related to the model given by formula2.
mweight.sim: sampled mixture weights $w[k[1],k[2]]$ of mixture components that had probabilities numerically higher than zero. Related to the model given by formula.
mweight\_2.sim: in the case of doubly-censored data, the same structure as mweight.sim, however related to the model given by formula2.
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. Related to the model given by formula.
mmean\_2.sim: in the case of doubly-censored data, the same structure as mmean.sim, however related to the model given by formula2.
gspline.sim: characteristics of the sampled G-spline (distribution of $(epsilon[i,1], epsilon[i,2])'$) related to the model given by formula. This file together with mixmoment.sim, mweight.sim and mmean.sim can be used to reconstruct the 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.
gspline\_2.sim: in the case of doubly-censored data, the same structure as gspline.sim, however related to the model given by formula2.
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. This file is related to the model given by formula.
mlogweight\_2.sim: fully created only if store$a2 = TRUE and in the case of doubly-censored data, the same structure as mlogweight.sim, however related to the model given by formula2.
r.sim: fully created only if store$r = TRUE. The file contains the labels of the mixture components into which the residuals are intrinsically assigned. Instead of double indeces $(k[1], k[2])$, values from 1 to $(2*K[1]+1)*(2*K[2]+1)$ are stored here. Function vecr2matr can be used to transform it back to double indeces.
r\_2.sim: fully created only if store$r2 = TRUE and in the case of doubly-censored data, the same structure as r.sim, however related to the model given by formula2.
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]]$). This file is related to the model given by formula.
lambda\_2.sim: in the case of doubly-censored data, the same structure as lambda.sim, however related to the model given by formula2.
beta.sim: sampled values of the regression parameters $beta$ related to the model given by formula. The columns are labeled according to the colnames of the design matrix.
beta\_2.sim: in the case of doubly-censored data, the same structure as beta.sim, however related to the model given by formula2.
Y.sim: 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: columns labeled loglik, penalty or penalty1 and penalty2, logprw. This file is related to the model given by formula. The columns have the following meaning. loglik $=$ $ -N(log(2*pi) + log(sigma[1]) + log(sigma[2])) -0.5*sum[i=1][N]( (sigma[1]^2*tau[1]^2)^(-1) * (y[i,1] - x[i,1]'beta - alpha[1] - tau[1]*mu[1,r[i,1]])^2 + (sigma[2]^2*tau[2]^2)^(-1) * (y[i,2] - x[i,2]'beta - alpha[2] - tau[2]*mu[2,r[i,2]])^2 )$ where $y[i,l]$ denotes (augmented) (i,l)th true 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], beta, 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 residuals 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)).$
logposter\_2.sim: in the case of doubly-censored data, the same structure as lambda.sim, however related to the model given by formula2.

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. (2006). Bayesian semi-parametric accelerated failure 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.

Examples

Run this code

## See the description of R commands for
## the population averaged AFT model
## with the Signal Tandmobiel data,
## analysis described in Komarek and Lesaffre (2006),
##
## R commands available in the documentation
## directory of this package as
## - see ex-tandmobPA.R and
##   http://www.karlin.mff.cuni.cz/~komarek/software/bayesSurv/ex-tandmobPA.pdf
##

Run the code above in your browser using DataLab