Bmc Medical Research Methodology Open Access a Semiparametric Bayesian Proportional Hazards Model for Interval Censored Data with Frailty Effects

Background: Multivariate analysis of interval censored event data based on classical likelihood methods is notoriously cumbersome. Likelihood inference for models which additionally include random effects are not available at all. Developed algorithms bear problems for practical users like: matrix inversion, slow convergence, no assessment of statistical uncertainty.


Background
Interval-censored survival data occur when the appearance of an event is assessed by means of an examination method that cannot tell the exact time of change in disease status, but only that the change has happened since the last examination. This is in contrast to the standard (naive) thinking that change in status coincides with the time of its first positive examination.
For example, the recurrence of a tumor during the followup of a treated cancer patient is an event which happens between two follow-up (FU) examinations. Often it is not possible to connect this event with an exact time (like time of first symptoms, time of first palpable presence, time of death,...) The information for a patient with recurrence is therefore as follows: it is known that up to a certain time (last FU examination at time t1) the patient is free of a recurrence. The recurrence happened between time t1 and t2 (present FU examination at time t2). This is less informative as the usual situation of right censoring. Non-or semi-parametric methods for interval censored data are not frequently used in clinical research papers. The reason may be that these methods are technically more complicated than standard survival methods based on exact or right-censored times. There is a a rich methodological body of methods and algorithms. But there is no easy to use software package in a popular statistical software environment.
The first frequentist work on interval censored event data is the model of Finkelstein. Finkelstein [1] first introduced a quasi semi-parametric model for regression analysis of interval-censored failure time data based on a finite dimensional maximum likelihood problem. For large datasets the information matrix is sparse and may not be invertible. Therefore, it is often to difficult to derive standard errors of the parameter estimates. Huang [2], Huang and Wellner [3] and Lin, Oakes and Ying [4] present alternative models for semi-parametric regression analysis on interval censored data. Pan [5] proposes the use of the iterative convex minorant algorithm (ICM) to handle the complex likelihood. He avoids the full information matrix and uses a bootstrap procedure to quantify the uncertainty of the inferred regression coefficients. Other approaches to likelihood estimates can be found in the work of Satten [6], or Huber, Solev and Vonta [7]. These methods assume independent observations. Bellamy et al. [8] extend parametric event time models to clustered and interval censored settings by introducing additive frailties to the linear predictor. Frailty comes into play when multiple events are considered for a given unit. Frailty can also count for unobserved covariates. Bellamy et al. implement their algorithm in existing commercial statistical computing software Bellamy et al. [8]. The authors model dependency between multiple events of a patient by frailty: doctor visits during the subject's time in study. Bellamy's idea is easily implemented in a Bayesian framework. WinBUGS [9] allows to implement analyses for interval censored data with frailty in the same model. A parametric approach via Weibull model or Accelerated Failure Time (AFT) model is easily realized. One can accommodate a frailty to the linear predictor part. But, the implementation of semiparametric proportional hazards models in WinBUGS is cumbersome (see Example Leuk in Example Volume 1 of the WinBUGS software).
Bayesian analysis of event data using non-or semi-parametric models started immediately after Cox [10] with work of Ferguson [11] and Kalbfeisch [12]. A summary of the current state of the art is given in Dey, Müller and Sinha [13] and Ibrahim, Chen and Sinha [14].
Many authors discuss a Bayesian approach to interval censored data with different forms of frailty. They demonstrate the advantageous combination of Bayesian inferential methods and MCMC sampling in this specific setting. The basic strategy is formulated by Härkänen, Virtanen and Arjas [15]. They avoid complex likelihoods by Bayesian data augmentation of censored lifetimes. They study a problem of dentistry where the dependence of event times between the teeth of a subject are accounted for by introducing subject-specific frailty parameters. Given such frailty parameters, the tooth lifetimes are inde-pendent. The hazard functions are defined non-parametrically by using piecewise constant functions. Finally, the interval censoring is handled by augmenting the data with unobserved exact event times. The authors also offer software described in Härkänen [16].
Another, more general methodology with similar ingredients has been published by Hennerfeind, Brezger and Fahrmeir [17]. In their work, Hennerfeind et al. do not consider explicitly interval censoring, however this is done by Kneib [18]. Kneib avoids the augmentation approach and handles the likelihood of interval censored observations by numerical integration techniques.
The purpose of the paper is to introduce an algorithm which uses a general software framework for statistical applications http://www.R-project.org and which is easy to apply to a large range of practical applications. We illustrate the theory behind the different modules of the algorithm and describe the specific MCMC strategies used to sample the posterior distribution. We restrict our consideration to specific submodels of Hennerfeind et al. and Kneib which are relevant for specific applications in clinical epidemiology. This is an advantage for many analysts. It bewares of the need to win a full understanding of the more general model classes.
Our approach is structured as follows: The working engine is a Bayesian model for right censored event data which is flexible with respect to the frailty structure underlying the data. We use standard frailty models like log-normal and gamma frailty. Second, we apply a fast data imputation algorithm: a piecewise exponential distributed event times which can easily sampled from piecewise constant hazard functions.
Being Bayesian we have to specify the priori belief in the shape of the baseline hazard which we assume to be smooth. A first order random walk is chosen.
The technical advantage with respect to the data augmentation given by a stepwise baseline hazard function may be offset by its possibly suboptimal fit to the true baseline hazard function. Therefore, we additionally implemented a spline based estimation of the log baseline hazard function. The draw back when performing the data imputation for the interval censored data is the need to linearly approximate the cumulative hazard function.
In Section 'Methods' the methods are described. Subsection 'The basic model' discusses the basic model for rightcensored event data. Subsection 'Sampling procedure' describes the sampling procedure needed for estimating the posterior distributions in the basic model. Subsection 'Extensions of the basic model' introduces the technique of data augmentation and the concept of frailty into the model. In Section 'Results' the results are given. Subsection 'Simulation' shows the results of a simulation and Subsection 'Examples' presents two examples, one of current status event data with frailty and one with right censored data and two frailty terms.

The basic model
This subsection sketches a Bayesian approach to a multivariate fixed effects proportional hazards model for right censored data. The specification of its posterior distribution needs the following ingredients: First, the likelihood of the observed data; second, specific prior distributions for regression parameters, hyperparameters, and baseline hazard; third, a Markov Chain Monte Carlo (MCMC) algorithm which can be used to sample the posterior distribution of the parameters of interest.
The data, based on a sample of size n, consists of the n tri- The likelihood contribution of the i-th single observation is given by where h(s) = ln[ 0 (s)] is the log transformed baseline hazard function. The function h will be modeled as a stepwise constant function as well as a cubic spline. The stepwise approach results piecewise exponential survival distributions. The use of cubic splines makes the estimated baseline hazard function smoother. Both concepts can be formulated as B-splines, see de Boor [19]. Therefore the time axis [0, ) is partitioned into disjoint intervals I k = [ k-1 ,  k ) for k = 1,..., K + 1.
The times  k , k = 1,..., K are chosen to create intervals with comparable information content, i.e. similar number of events. To this end a Kaplan Meier estimator is applied to the right censored data. It is interpolated linearly. Spacing the range of the resulting survival curve into K equal parts and using the inverse function of the modified Kaplan Meier estimator for the interval boundaries gives the times  k , k = 1,..., K. We observed that a large K destabilized the mixing properties of the MCMC chain by resting in a state for long periods. Therefore, we chose K as large as possible by keeping an acceptable mixing of the chain.
The priors for the components of the vector  will be independently normal distributed with mean 0 and a small precision  (Precision is defined as the inverse of the variance:  = ), see Gamerman [20].
The prior for the coefficients of the spline function h, the approximated version of the log-transformed baseline hazard function, will be a first order process which gives prior information on smoothness. This is a Bayesian Pspline approach, see Lang and Brezger [21].
The log baseline hazard function h is approximated by Bspline functions, , where b j ,·(t) denotes a B-spline function of degree j. In our approach j is zero for the step functions or three for the cubic splines. The first order process is defined as h k = h k-1 + k with k Ñ

Sampling procedure
Sampling for the parameter vector Aitkin and Clayton [22] pointed out that the proportional hazards model can be interpreted as a generalized linear model by writing the likelihood function L in the form With a given cumulative baseline hazard  0 (t i ) this is the likelihood function of a Poisson sample where the observations are the survival indicators  i , the link function is the canonical link "log" and there is the offset ln( 0 (t i )).
Gamerman [23] describes how to sample effectively the vector of covariates in generalized linear mixed models in a block updating step. This is a combination of the iterated least squares method (IWLS) with a Metropolis-Hastings sampling.
For the prior a weak informative normal distribution N(0, R -1 ), R =  2 I n is chosen. The start of the iteration is with  =  0 and t = 1. One has to sample * from N(m (t) , where * is accepted with probability ( (t-1) , *) and then  (t) = *, else  (t) =  (t-1) and t is increased by one.

Sampling for the baseline hazard
With the given structure of the log baseline hazard function one has to sample the spline coefficients h from a Gaussian Markov Random Field (GMRF). Here we follow Knorr-Held and Rue [24] and Rue [25] and sample the log-baseline hazard in one step. The posterior of h is Knorr-Held and Rue [24] propose to approximate the exponent by a quadratic form, to use the resulting Gaus-sian random field (multivariate normal distribution) to sample a proposal for h, and to accept or to reject the proposal according to a Metropolis-Hastings step. The cumulative baseline hazard which is part of the exponent can be calculated in a closed form in the case of the degree zero B-splines (step functions). A detailed calculation is given in the appendix [see Additional file 1]. The cumulative baseline hazard in the case of the B-splines of degree three is approximated by the trapezoidal rule because it can not be integrated analytically. The trapezoidal rule results in complex terms which contain the exponent of linear terms of h. The calculation of a good quadratic approximation to this terms allows to derive the multivariate normal distribution.

Sampling for the dispersion parameters
For the dispersion parameters and a flat gamma prior with rate  and shape  is chosen. The full conditional distribution of is again gamma distributed and has rate and shape and the full conditional distribution of has rate and shape .

Extensions of the basic model
Data augmentation is used to impute unobserved event times which creates right censored data from interval censored data. An unit specific random effect or frailty term is introduced to the proportional hazards model to account for potential clustering of event times within a statistical unit.

Chained Data augmentation
The chained data augmentation algorithm applied to interval censored data imputes candidates of possible event times conditional to the model and the data observed. To obtain the posterior (|Data), where  = (, h,  0 ,  1 ), one proceeds iteratively by generating right censored imputations of event times T from the predictive distribution (T|, Data), and calculates (|T, Data), from the augmented data, see Tanner [26].
This suggests that one may draw a value of the parameter vector of interest * from (|T*, Data) where T* is drawn from (T|, Data). In the given case the vector  contains all information on , h,  0 and  1 . To initialize the sampling the interval censored event data is treated as right censored (event times in finite intervals are set to the interval midpoint, event time in infinite intervals are considered as right censored) and the basic model (Subsection 'The basic model') is applied. The imputation is applied to all interval censored events. The times of the right censored events are unchanged. The imputation is a straightforward sampling procedure based on the individ- ual survival functions which can be calculated from  for every event in every unit. The individual survival function has to be conditioned on the interval in which the event happens, see also Bebchuk and Betensky [27].
The cumulative baseline hazard  0 (t) is piecewise linear within the defined intervals for B-splines of degree zero. This implies a piecewise exponential survival distribution and a straight forward sampling of the imputed times. The piecewise linear approximation for B-splines of degree three makes a good proposal for the imputed times. Unfortunately the numerics to calculate an acceptance score relies on an approximation of the numerically not treatable true distribution function. Therefore we skipped the acceptance step.

Potential clustering of event times
Unit-specific random effects are used to handle clustering of event times within statistical units. Cluster specific covariates are introduced which takes the value 1 for every observation corresponding to the relevant cluster and 0 else. It is assumed that observation i belongs to cluster j(i). Different frailty distributions can be chosen.

Log-normal frailty
Log-normal frailty is difficult to handle in frequentist frailty models, see Hougaard [28]. In the Bayesian setting random effects (frailty terms) are treated like regression coefficients In case of log-normal frailty these coefficients have a normal prior with mean 0 and variance ,  ~ N(0, ). A non-informative gamma prior is assumed for .
For the posterior holds The sampling follows the ideas of Knorr-Held and Rue [24] as described in Subsection 'Sampling procedure'. The posterior of   is again gamma distributed.

Gamma frailty
A gamma frailty is proposed by Clayton [29]. The gamma frailty distribution offers technical advantages in the maximum likelihood framework because it allows to express the likelihood as a Laplace transform, see Feller [30]. The cumulative hazard function can be written as The posterior of z is with the gamma prior Ga(, p) again and .
For the gamma prior  = p is chosen such that the expected value of the prior is one. The parameter  is updated by a random walk proposal and a Metropolis-Hastings step.

Simulation
The purpose of the simulation is to show that the theoretical Bayesian framework gives the expected results given a known true model. Let us summarize shortly the essence of the Bayesian paradigm applied to the hazard function.
We define a prior distribution for the log-baseline hazard as a random walk with mean 0 and a specific covariance structure. Within the Bayesian framework the data changes the prior distribution into the posterior distribution. The Bayesian procedure changes the random walk structure of the log-baseline hazard in a two-fold way: 1.) The constant mean 0 is transformed to a specific form of a drift. The drift of the posterior distribution (which is a random walk) of the log-hazard function inferred by our procedure in the simulation study is shown in the accompanying figure and compares the inferred drift to the true log-baseline hazard function used in the data generation.
2.) The quite variable covariance structure of the prior distribution (large variances) allows that the data have a strong influence on shaping the posterior distribution.
Data with a lot of information on the process will remove the variable covariance structure and will create a narrow channel around the true drift where realizations of the random walk can be found which represents the posterior of the log-baseline hazard function. In spite the fact, that the choice of the prior distribution is a random walk with the constant mean 0, the choice of a variable covariance structures allows the estimation procedure to come up with a posterior distribution (again a random walk) whose drift describes sufficient closely the true shape of the log-baseline hazard.
The choice of the mean 0 drift for the prior of the log-baseline hazard is motivated by the idea, that if there is a deviation from a constant hazard the data should produce it in the posterior distribution. Reducing the variance in the covariance structure of the prior distribution has the effect Of course we could start by giving a drift to prior distribution which is close to the supposed true form of the logbaseline hazard. But we rejected the idea of using informative priors. If the covariance structure of the prior is quite variable then the actual drift of the prior distribution has not a strong influence of the drift calculated for the posterior distribution.
Data from a Weibull model is used to validate the proposed MCMC procedure. At the one hand the Weibull model is parametric at the other hand it fits many practically relevant situations. The hazard function is given by (t) = (t/b) a with shape a and scale b. The validation model has four covariates and a gamma frailty. The covariates X 1 and X 2 are binary. X 1 is balanced: 50% of the subjects in our virtual population take value one. The covariate X 2 is unbalanced: 70% of the subjects in our virtual population take value one. The variable X 3 is restricted to a well defined range (uniformly distributed on [-1, 1]). The variable X 4 is standard normal distributed. One thousand observations which are randomly assigned to clusters. The cluster size is random but the number of clusters is fixed: 100, 200, and 500. Each cluster carries a gamma frailty. The frailty is sampled from a gamma distribution with equal rate and shape parameters, Ga(q, q), with q  {0.5, 1, 2}. The mean of the gamma distribution is one.
For each subject the interval [0, t max ) is randomly divided into five intervals by uniformly sampling five random numbers on [0, t max ). If the individual event time is at least t max the event time is taken to be right censored at t max , else the event time is interval censored with the interval it falls in.
In order to make the censoring mechanism independent of the individual survival process t max has to be chosen independent of the individual who has to be censored. Therefore the value of t max is the 0.9-quantile of a specific Weibull distribution with shape a = 0.75 (see above) and as scale the median of the individual scales which are derived from the specific covariate configuration and the random effect. The time axis was divided into 50 intervals using the Kaplan Meier procedure as described in Subsection 'The basic model'.
The chain runs through 50 000 cycles. The cubic spline model requires 70+#clusters parameters to be sampled per cycle. The Raftery-Lewis diagnostic, see Mansmann [31], was calculated for regression and baseline hazard coefficients. The Raftery-Lewis diagnostics indicates a burn in of 40 000 cycles to reach convergence for all parameters. 10 000 additional samples of the parameters were taken and thinned by the factor ten. The remaining 1000 samples were used for the analysis.
The trace of  for one of the scenarios studied can be seen on Figure 1.
The estimates for the regression coefficients  and the frailty parameter q and their 95% credibility intervals are given in Table 1.
The bias for regression coefficients is in general moderate. A pattern for the stronger biased results for the regression coefficients can not be seen in the nine scenarios. The bias for the frailty parameter is in general quite small.
We supposed that increasing the number of clusters (decreasing cluster size) would improve the information content of the data by allowing for more independent observations. Therefore, we expected a clear improvement from 100 clusters to 500 clusters with respect to the precision of estimates. A similar effect was expected from changing the frailty. Reducing the dependency within a cluster should have a similar effect. It was surprising not to see clear trends in the simulation results.
We compared the bias in our strategy to the bias produced by the ICM algorithm proposed by Pan [5], which we provide as R-package intcox. This algorithm is considered as one of the best likelihood based estimation procedures for multivariate interval censored survival data, see Zhang and Jamdhidian [32]. Unfortunately the ICM algorithm is not able to handle frailty. We estimate the parameters by ignoring the clustering but apply clustered bootstrap to calculate their confidence intervals. Table 2  Traces of the regression coefficients ; 100 clusters, rate and shape of the frailty gamma distribution q = 1 Figure 1 Traces of the regression coefficients ; 100 clusters, rate and shape of the frailty gamma distribution q = 1.

Examples
In general flat priors are chosen in the R-package sur-vBayes. The precision of the prior for , the rate and shape of the priors for  0 and  1 and the precision of the random walk for the gamma frailty are set to 0.0001. Possible deviations from these values are indicated below.

Aneurisms
Meisel et al. [33] present data on the shrinkage of aneurisms associated with cerebral arteriovenous malformations (cAVM) after embolization treatment. The time to a shrinkage of the aneurism to below 50% of the baseline volume was of interest. Several patients had multiple aneurisms. Each patient was inspected at a random inspection time. Thus, the data is current status censored, the coarsest form of interval censoring (see example 3.4 in Klein and Moeschberger [34]).
Two covariates were considered: the degree of cAMV occlusion by embolization (dichotomized at 50%, variable mo) and the location of the aneurism, whether at the midline arteries or at other afferent cerebral arteries, vari-able loc. Multiple aneurisms were observed per patient. In this case the aneurisms share the same "environment" and may not behave independently.
The data set is analyzed with the model described in Section 'Methods'. The log baseline hazard is modeled by cubic splines as well as by a step function. There were phases of no sampling in the Metropolis-Hastings steps for the log baseline hazard coefficients. The problems could be resolved by choosing a not too flat prior with rate and shape q = 0.001 (cubic splines) or q = 0.01 (constant splines) for the smoothness parameter . The precision for the random walk for the updating of the gamma frailty parameter q was set to 0.01 in the case of the constant splines. 20 000 samples of the parameters were taken and thinned by the factor ten after a burn in of 80 000. The remaining 2000 samples were used for the analysis. There is still no convergence for the last parameters of the log base line hazard while the results for the regression parameters are stable (see Table 3). Inspecting the data  shows that only 54 of 149 aneurysms showed shrinkage. This may impair precise estimation of the baseline hazard function for larger times (above 3 years).
The estimates and the 95% credibility intervals of the log baseline hazards for both models are shown in Figure 3.
The model with cubic splines gives the smoother fit.
We can compare the Bayesian results with the related maximum likelihood analysis using the Iterated Convex Minorant algorithm as introduced by Pan [5]. The analysis of the Aneurism data set with different likelihood based tools is given in the Vignette accompanying the R-package intcox. The results for the regression parameters are shown in Table 3. The model proposed by Pan [5] considers all observations as independent. This may explain the smaller absolute effect estimates.

Colon cancer
The Tumour Register Munich (TRM) collects the data of all cancer patients in Munich and the surrounding southern Bavaria. Here we use data from patients with colorectal carcinoma which includes 9985 patients without metastases collected from 1988 to 2004. Their age distribution at diagnosis is as follows: 621 were less than 50 years old, 1541 between 50 and 60, 2884 between 60 and 70, 2992 between 70 and 80 years old, and 1947 older Estimated and true log baseline hazards     Table 4.
The hazard of death increases with increasing age, pT and pN category and with residual tumour and is lower in women and if a chemotherapy could be applied.
The posterior means of the frailty coefficients of the clinics and the 95% credibility interval belonging to them are shown in Figure 4.
The credibility intervals of the clinics with high volume are the smallest. There are only two clinics whose credibility intervals do not include the value one. This value indicates that there is no frailty. Clinics with a high volume tent to have a lower frailty (clinics 7, 31) but this is not true for all of them. On the other hand there is no clear pattern for clinics with lower volume. Clinics with low patient volume may have low frailty (clinics 20, 45) as well as higher frailty (clinics 26, 40).
Transforming the frailty into ranks allows to study the posterior distribution of the rank of a specific clinic. These distributions ar shown for some selected clinics in Figure  5.
The rankings are quite stable for clinics with extreme frailties. Clinics with intermediate frailties range in general over all ranks.

Discussion
Interval censored event data is a natural documentation scheme for many observational studies. Often, events of interest are not obvious and can only be detected with special diagnostic procedures during follow-up examinations. This is not only of interest for endpoints like disease free survival in oncological studies but also of a broader clinical interest as demonstrated in the aneurysm example. Interval censoring may also help to avoid special forms of bias in epidemiological studies like recall biassee example 3.4 in Klein and Moeschberger [34]. Including random effects into the model handles unobserved covariates, multiple events, or clustering. Interval censoring and frailty may result in a complex likelihood which is not easy to maximize.
Bellamy et al. [8] showed that within the framework of Weibull models the likelihoods which incorporate interval censoring and simple frailty can be solved by standard software. Formulating the aneurysm example as a Weibull problem in WinBUGS is also possible with a few lines of code: The treatment is more complex if the statistician decides to use a semi-or non-parametric approach for the hazards model. The example leuk in volume I of the examples which come with the WinBUGS software demonstrates the complexity of a code for a semiparametric Cox-model. Its extension for frailty and interval censoring makes things even more complex. Furthermore, the WinBUGS software is quite restricted in the choice of updating schemes for the MCMC procedure which improve mixing effects of the chains. Therefore, this paper proposes a Bayesian concept for the simultaneous treatment of frailty and interval censored event times. Data augmentation reduces the interval censored situation to the right censored case and allows to handle frailty in the framework of a semi-parametric proportional hazards model for right censored data. A MCMC procedure was introduced to estimate the parameters of interest. Gamerman's block updating was used to sample the regression coefficients. A block updating based on ideas of Knorr-Held and Rue [24] and Lang and Brezger [21] allowed a simultaneous sample of the baseline hazard function. The incorporation of frailty into the MCMC scheme was straightforward. Data augmentation could make efficient use of piecewise exponential distributions. This is a consequence from a stepwise log baseline hazard function or a piecewise linear approximation to the cumulative hazard function when the log baseline hazard is modeled by cubic splines.
The algorithm is tailored to a wide class of event data problems which are typical for studies in clinical epidemiology, quality control, prognosis, and epidemiological risk assessment. The ideas on which the algorithm is built are also discussed by other groups like Härkänen et al. [15], Kneib [18], Komarek, Lesaffre, Härkänen, Declerck and Virtanen [36], Hennerfeind et al. [17]. The focus of our paper is not to present novel methodology, but to give a full methodological account on a practical algorithm and the accompanying software.
A software which offers comparable functionality is BITE by Härkänen [16]. It is designed for the analysis of event history data using flexible hierarchical models and Bayesian inference. BITE is a stand-alone software while our package uses the full functionality of the R-environment. The handling of BITE is not straight forward. While sur-vBayes is independent of the machine platform, BITE is Clinic's identification Frailty written for Unix-like operating systems such as Linux. It can also be run on 32-bit Microsoft systems, but a Unixemulator is needed. Its output files have to be processed with PERL scripts to be read into other statistical software like R or CODA which is needed to present the results of the analysis.
Two examples were worked out. The first example fits in the structure of the model presented in Section 'Methods'. The data consist of multiple current status observations on patients. The example was run with cubic and constant P-splines. The cubic P-spline estimate of the baseline hazard shows comparable results to the constant P-spline in Sampling of the ranks of selected clinics the part of the time spectrum with the most events. The regression parameter estimates also remain comparable. We could compare our results with estimates based on maximum likelihood theory. The Bayesian results are more pronounced compared to the maximum likelihood estimates because the frailty effect is included in the model. This example is in line with the typical problems also discussed by Härkänen [16] and observed in applications of BITE.
The second example presents a two level frailty structure which includes individual frailty of patients but also frailty related to single hospitals. The data describe the survival of colon cancer patients collected in a tumor registry (Engel et al. [35]). The focus of the analysis was on the estimates of the frailty of each individual clinic. A second frailty term is needed to adjust for unobserved individual prognostic covariates. The result indicates that there are remarkable clinic effects in the survival of colon cancer patients. The credibility intervals for the frailties of the clinics with the highest volume are the smallest but for all clinics reasonable estimates of the frailties are obtained. The analysis is of interest for the benchmarking of institutions.
A simulation study based on a Weibull model with three different gamma frailties Ga(q, q), q  {0.5, 1, 2}, was used to validate our approach. One thousand events were randomly distributed among 100, 200 or 500 patients (conditioned on the fact that each cluster received at least one event). For the regression parameters we got sensible estimates. The bias for the regression coefficients is in general lower than 20%. No clear pattern for the bias could be detected. A relevant result of the simulation study was the reliable estimation of the frailty parameter which is of primary interest in many problems. The baseline hazard is fitted best for larger frailty values (q = 2).
In general we got the impression that the modeling of interval censored event data needs a very careful analysis. Especially convergence of the sampler chains may be not straightforward in small datasets. The proposed procedure needs data sets of sufficient size. Then it will produce in general reasonable estimates of the parameters. But there is always some bias possible with no clear pattern when it happens. Our concept is implemented into the R-package survBayes and is available at the CRAN depository http:// www.R-project.org.

Conclusion
Occurrence of diseases, their progression and death constitute complex event patterns in time. There are two pending problems: (1) There are interval censored events which can only be observed between two examinations.
(2) Unspecific and not observable effects have to be taken into account by more or less complex random effect structures. We present a practical solution to both challenges by the software bayesSurv. The paper presents the theory which motivates the algorithms, evaluates the algorithm by a simulation study and applies our approach to two relevant clinical examples. The presented software supports the solution of complex analyses for event data in many fields of clinical epidemiology as well as health services research. It helps to build quantitative models for complex observational data and contributes to an improved quantitative modeling of individual disease histories.