Generalized Gamma Frailty and Normal Random Effects Model for Repeated Time-to-Event Data

1. Introduction

Clustered time-to-event data are commonly encountered in contemporary medical and statistical studies. It is important to accommodate data hierarchies that can result from repeated measurements of survival outcomes on the same subject [2]. In addition, overdispersion [3] as well as censored observations are most likely to be found in the data. Therefore, it is crucial to build a suitable model to capture the aforementioned features present in the data.

Models that simultaneously deal with issues of overdispersion and data hierarchy in survival analysis is not so common. Molenberghs et al. [4] proposed a general framework to model non-Gaussian outcomes, accommodating overdispersion and data hierarchies through two separate sets of random effects. The case of time-to-event data is one specific application under their broad class of models. Molenberghs et al. [1] extended and presented a model for repeated, overdispersed time-to-event data that are also subject to censoring. These authors [1] combined conjugate random effect for overdispersion with generalized linear mixed model (GLMM [5,6]) in which normal random effects are embedded within the linear predictor for clustering. Two estimation methods-a partial marginalization approach to full maximum likelihood and a pairwise-likelihood version of pseudo-likelihood- were presented and compared, based on a limited simulation study as well as two real data analysis. Molenberghs et al. [1] chose gamma distribution to be the frailty random effects to capture overdispersion. This convenient choice was motivated through the concept of conjugacy. However, convenience does not assure a good fit to the observed data. To overcome this limitation, we propose here a new model that is flexible to provide adequate fit to the repeated and overdispersed time-to-event data, while also accommodating right censoring. The frailty random effect in this paper is assumed to follow a generalized gamma distribution, previously studied by Balakrishnan and Peng [7]. The generalized gamma distribution has one more parameter than the gamma distribution, and is therefore more flexible and less parametric [7]. Moreover, this distribution also includes the well-known gamma, lognormal and Weibull models all as special cases.

The rest of this paper is organized as follows. Key ingredients for our model formulation and the new combined model are the subjects of Section 2. Section 3 discusses an estimation method for the proposed model. In Section 4, an extensive simulation study is conducted to investigate the performance of the proposed estimation method and to compare the performance of the proposed model with the existing one [1]. The illustration of the proposed model with two motivating data sets is made in Section 5. Finally, some concluding remarks are made in Section 6.

2. Generalized Gamma Frailty and Normal Random Effects Model

2.1. Generalized Gamma Distribution

The generalized gamma (GG) distribution has been applied in many areas of statistical applications, due mainly to its flexible form. A random variable $Y>0$ is said to follow the generalized gamma distribution if its density function is as follows:

$$f(y;q,\sigma ,\lambda )=\left\{\begin{array}{cc}\left|q\right|{\left(q^{-2}\right)}^{q^{-2}}{\left(\lambda y\right)}^{q^{-2}(q/\sigma )}exp[-q^{-2}{\left(\lambda y\right)}^{q/\sigma }]/[\Gamma \left(q^{-2}\right)\sigma y],\hfill & q\ne 0\hfill \\ {\left(\sqrt{2\pi }\sigma y\right)}^{-1}exp\{-{[log\left(\lambda y\right)]}^2/\left(2{\sigma }^2\right)\},\hfill & q=0\hfill \end{array}\right.,$$

(1)

where -$\infty <q<\infty$ and $\sigma >0$ are shape parameters and $\lambda >0$ is a scale parameter. It was first developed by Stacy [8] and a review of the generalized gamma distribution is available in [9]. This distribution was studied by Prentice [10], wherein a parameterization was proposed that allowed the maximum likelihood estimates to be computed using standard algorithms. The application of the generalized gamma distribution in an accelerated failure time model was addressed by Lawless [11]. Yamaguchi [12], Peng et al. [13] and Balakrishnan and Pal [14] investigated the use of the distribution in the context of cure rate model. Shin et al. [15] proposed a statistical model for speech signals based on the generalized gamma distribution. Cox et al. [16] provided a taxonomy of hazard functions for the generalized gamma distribution, while Balakrishnan and Peng [7] proposed a frailty model using the generalized gamma distribution as the frailty distribution.

Generalized gamma distribution is a broad family of distributions that includes exponential, gamma, and Weibull distributions as subfamilies, and lognormal as a limiting distribution. For example, the distribution reduces to the gamma distribution when ${\displaystyle \frac{q}{\sigma }}=1$; it is the lognormal distribution when $q=0$ (the limit of $f(y;q,\sigma ,\lambda )$ when $q\to 0$); and it is the Weibull distribution when $q=1$. As the generalized gamma distribution can offer a considerable amount of flexibility, in this way, it can be used to capture more features of the data that may be missed when any of its special cases is used instead. The mean of the generalized gamma distribution given in (2), exists only when $q>{\displaystyle \frac{1}{\sigma }}$, and is given by

$${\displaystyle \frac{\Gamma (q^{-2}+\frac{\sigma }{q})}{\Gamma \left(q^{-2}\right){\left(q^{-2}\right)}^{\sigma /q}\lambda }}.$$

(2)

If the mean is set to one, then

$$\lambda ={\displaystyle \frac{\Gamma (q^{-2}+\frac{\sigma }{q})}{\Gamma \left(q^{-2}\right){\left(q^{-2}\right)}^{\sigma /q}}}.$$

(3)

We can then substitute (3) to (1), and also get the variance as

$${\displaystyle \frac{\Gamma (q^{-2}+2\frac{\sigma }{q})\Gamma \left(q^{-2}\right)}{{\Gamma }^2(q^{-2}+\frac{\sigma }{q})}}-1.$$

(4)

2.2. Weibull-Generalized Gamma-Normal Model

A random variable $\mathit{Y}$ follows an exponential family of distributions if its density is of the form

$$\begin{array}{c}\hfill f\left(y\right)=f\left(y\right|\eta ,\varphi )=exp\{{\varphi }^{-1}[y\eta -\psi \left(\eta \right)]+c(y,\varphi )\},\end{array}$$

(5)

where $\eta$ is termed as the natural parameter, $\varphi$ is called the dispersion parameter, and $\psi (·)$ and $\mathit{c}(·)$ are some known functions. In the regression context, the model needs to incorporate measured covariates to explain variability among outcome values, and this is referred to as a generalized linear model. $Y_1,...,Y_N$ is a set of independent outcomes, and $x_1,...,x_N$ represent the corresponding p-dimensional vectors of covariate values. Then, $Y_i$’s are assumed to have densities $f\left(y_i\right|{\eta }_i,\varphi )$ that belong to the exponential family in (5), while ${\mu }_i$’s are modeled as functions of the covariates. It is assumed that ${\mu }_i=h\left(x_i^{\prime }\xi \right)$, with $\xi$ is a vector of p fixed unknown regression coefficients and $h^{-1}(·)$ is called the link function. In most applications, the natural link function is used, i.e., $h(·)={\psi }^{\prime }(·)$ with ${\mu }_i=h\left(x_i^{\prime }\xi \right)={\psi }^{\prime }\left(x_i^{\prime }\xi \right)={\psi }^{\prime }\left({\eta }_i\right)$. This is equivalent to assuming ${\eta }_i=x_i^{\prime }\xi$. As overdispersion is often present in time-to-event data, we can use a two-stage approach to construct an overdispersed model. Consider a distribution for the outcome, given a random effect, $f\left(y_i\right|{\theta }_i)$, and assume a distribution for the random effect, $f\left({\theta }_i\right)$. The marginal distribution is then produced as

$$\begin{array}{c}\hfill f\left(y_i\right)=\int f\left(y_i\right|{\theta }_i)f\left({\theta }_i\right)d{\theta }_i.\end{array}$$

(6)

If $f\left(y_i\right|{\theta }_i)$ follows a Weibull distribution, and $f\left({\theta }_i\right)$ follows a gamma distribution, then we obtain a Weibull-Gamma model. The choice of gamma distribution for the random effect is motivated through the concept of conjugacy, so that the marginal distribution has a closed-form expression. Although gamma distribution is convenient to use, the convenience, however, does not assure that the fit will be good. The generalized gamma distribution has one more parameter than the gamma distribution, and so can be more flexible enabling better modelling of the data. This can be generalized to clustered data in survival analysis, where $Y_{ij}$ can be used to denote the $j^{th}$ outcome value in cluster i, for $i=1,2,3,...,N$ and $j=1,2,...,n_i$. On the other hand, it is possible to include normal random effects in the linear predictor of the generalized linear model, giving rise to the family known as generalized linear mixed model. Conditional on q-dimensional random effects ${\mathbf{b}}_{\mathbf{i}}\sim N(0,D)$, the outcomes $Y_{ij}$ are independent with exponential-family densities of the form

$$\begin{array}{c}\hfill f_i\left(y_{ij}\right|{\mathbf{b}}_{\mathbf{i}},\xi ,\varphi )=exp\{{\varphi }^{-1}[y_{ij}{\lambda }_{ij}-\psi \left({\lambda }_{ij}\right)]+c(y_{ij},\varphi )\},\end{array}$$

(7)

with

$$\begin{array}{c}\hfill \eta \left[{\psi }^{\prime }\left({\lambda }_{ij}\right)\right]=\eta \left({\mu }_{ij}\right)=\eta \left[E\left(Y_{ij}\right|{\mathbf{b}}_{\mathbf{i}},\xi )\right]=x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}},\end{array}$$

(8)

where $x_{ij}$ and $z_{ij}$ are p-dimensional and q-dimensional vectors of known covariate values, $\eta (·)$ is a known link function, $\xi$ is a p-dimensional vector of unknown fixed regression coefficients, and $\varphi$ is a scale (overdispersion) parameter. $f\left({\mathbf{b}}_{\mathbf{i}}\right|D)$ is the density of $N(0,D)$ distribution for the random effects ${\mathit{b}}_{\mathit{i}}$. Now, combination of overdispersion and normal random effect leads to the following combined model

$$\begin{array}{c}\hfill f_i\left(y_{ij}\right|{\mathbf{b}}_{\mathbf{i}},\xi ,{\theta }_{ij},\varphi )=exp\{{\varphi }^{-1}[y_{ij}{\lambda }_{ij}-\psi \left({\lambda }_{ij}\right)]+c(y_{ij},\varphi )\}.\end{array}$$

(9)

The conditional mean is $E\left(Y_{ij}\right|{\mathbf{b}}_{\mathbf{i}},\xi ,{\theta }_{ij})={\mu }_{ij}^c={\theta }_{ij}{k}_{ij}$, where ${\theta }_{ij}\sim G_{ij}({\nu }_{ij},{\sigma }_{ij}^2)$ and $k_{ij}=g(x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}})$. The relationship between mean and natural parameter now is ${\lambda }_{ij}=h\left({\mu }_{ij}^c\right)=h\left({\theta }_{ij}{k}_{ij}\right)$. Random effects ${\theta }_{ij}$ can capture overdispersion, and $k_{ij}$ can be considered as the generalized linear mixed model component that can capture the in-between cluster effect. In the present work, we consider the general Weibull model for repeated measures, with both generalized gamma and normal random effects, of the following hierarchical form:

$$\begin{array}{c}\hfill f\left(y_i\right|{\theta }_i,{\mathbf{b}}_{\mathbf{i}})=\prod _{j=1}^{n_i}\lambda \rho {\theta }_{ij}{y}_{ij}^{\rho -1}exp(x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}})exp(-\lambda y_{ij}^{\rho }{\theta }_{ij}{e}^{x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}}}),\end{array}$$

(10)

$$\begin{array}{c}\hfill f\left({\mathbf{b}}_{\mathbf{i}}\right)={\displaystyle \frac{1}{{\left(2\pi \right)}^{q/2}{\left|D\right|}^{1/2}}}{e}^{-{\displaystyle \frac{1}{2}}{{\mathbf{b}}_{\mathbf{i}}}^{\prime }{D}^{-1}{\mathbf{b}}_{\mathbf{i}}},\end{array}$$

(11)

$$\begin{array}{c}\hfill f\left({\theta }_i\right)=\prod _{j=1}^{n_i}{\displaystyle \frac{\left|q\right|{\left(q^{-2}\right)}^{q^{-2}}{\left(p{\theta }_{ij}\right)}^{q^{-2}(q/\sigma )}exp[-q^{-2}{\left(p{\theta }_{ij}\right)}^{q/\sigma }]}{[\Gamma \left(q^{-2}\right)\sigma {\theta }_{ij}]}}.\end{array}$$

(12)

To make the parameters in the model identifiable, the mean of the generalized gamma distribution is usually set to be one. As the mean of the generalized gamma distribution is ${\displaystyle \frac{\Gamma (q^{-2}+\frac{\sigma }{q})}{\Gamma \left(q^{-2}\right){\left(q^{-2}\right)}^{\sigma /q}p}}$, with the mean set to one, we then have $p={\displaystyle \frac{\Gamma (q^{-2}+\frac{\sigma }{q})}{\Gamma \left(q^{-2}\right){\left(q^{-2}\right)}^{\sigma /q}}}$.

3. Likelihood Function and Estimation Method

In a frailty model, one of the most common methods used for estimating the unknown parameters of the model is the maximum likelihood method, as used, for example, in Balakrishnan and Peng [7]. In this paper, we make use of the marginal likelihood approach as a direct way for finding the maximum likelihood estimates of the unknown parameters in the combined generalized gamma and normal random effect model. However, one of the difficulties in the parameter estimation is that the likelihood function involves many integrals which render it to be intractable. For this reason, we can approximate the integral in the likelihood function directly and then maximize it to obtain the maximum likelihood estimates of the parameters of interest. The likelihood contribution of subject i (cluster i) is

$$f_i\left(y_i\right|\vartheta ,D,{\vartheta }_i,{\Sigma }_i)=\int \prod _{j=1}^{n_i}{f}_{ij}\left(y_{ij}\right|\vartheta ,{\mathit{b}}_{\mathit{i}},{\theta }_i)f\left({\mathit{b}}_{\mathit{i}}\right|D)f\left({\theta }_i\right|{\vartheta }_i,{\Sigma }_i)d{\mathit{b}}_{\mathit{i}}d{\theta }_i,$$

(13)

while the likelihood function is given by

$$L(\vartheta ,D,{\vartheta }_i,{\Sigma }_i)=\prod _{i=1}^N\int \prod _{j=1}^{n_i}{f}_{ij}\left(y_{ij}\right|\vartheta ,{\mathit{b}}_{\mathit{i}},{\theta }_i)f\left({\mathit{b}}_{\mathit{i}}\right|D)f\left({\theta }_i\right|{\vartheta }_i,{\Sigma }_i)d{\mathit{b}}_{\mathit{i}}d{\theta }_i.$$

(14)

A numerical method, with the use of Monte Carlo simulation, is implemented to approximate the integral $I_i$ in the likelihood function, where

$$I_i=\int \prod _{j=1}^{n_i}{f}_{ij}\left(y_{ij}\right|\vartheta ,{\mathit{b}}_{\mathit{i}},{\theta }_i)f\left({\mathit{b}}_{\mathit{i}}\right|D)f\left({\theta }_i\right|{\vartheta }_i,{\Sigma }_i)d{\mathit{b}}_{\mathit{i}}d{\theta }_i,$$

(15)

and it is approximated as

$$\begin{array}{cc}\hfill I_i& =E_{{\theta }_i}\left[E_{{\mathit{b}}_{\mathit{i}}}\left[\prod _{j=1}^{n_i}{f}_{ij}\left(y_{ij}\right|\vartheta ,{\mathit{b}}_{\mathit{i}},{\theta }_i)\right]\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{N_1\times N_2}}\sum _{k_1=1}^{N_1}\sum _{k_2=1}^{N_2}\prod _{j=1}^{n_i}{f}_{ij}\left(y_{ij}\right|\vartheta ,{\mathit{b}}_{\left({\mathit{k}}_{\mathit{1}}\right)},{\theta }_{\left(k_2\right)}),\hfill \end{array}$$

(16)

where ${\mathit{b}}_{\left({\mathit{k}}_{\mathit{1}}\right)}$, $k_1=1,2,...,N_1$, and ${\theta }_{\left(k_2\right)}$, $k_2=1,2,...,N_2$, are realizations of normal random variables and generalized gamma random variables, respectively. Now, to accommodate right-censoring in the data, for each j, we integrate the conditional distribution over the time interval [$C_{ij},+\infty$]:

$$S\left(C_{ij}\right|{\theta }_{ij},{\mathit{b}}_{\mathit{i}})={\int }_{C_{ij}}^{\infty }{f}_{ij}\left(y_{ij}\right|{\mathit{b}}_{\mathit{i}},{\theta }_i)d{y}_{ij}.$$

(17)

For the case of Weibull-GG-Normal model, we have

$$\begin{array}{cc}\hfill S\left(C_{ij}\right|{\theta }_{ij},{\mathit{b}}_{\mathit{i}})& ={\int }_{C_{ij}}^{\infty }\lambda \rho {\theta }_{ij}{y}_{ij}^{\rho -1}exp(x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}})exp(-\lambda y_{ij}^{\rho }{\theta }_{ij}{e}^{x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}}})d{y}_{ij}\hfill \\ \hfill & =exp(-\lambda e^{x_{ij}^{\prime }\xi +z_{ij}^{\prime }{\mathbf{b}}_{\mathbf{i}}}{\theta }_{ij}{C}_{ij}^{\rho }).\hfill \end{array}$$

(18)

The maximum likelihood estimates of the parameters are then obtained numerically through quasi-Newton-Raphson algorithm. The standard errors are also obtained numerically from this approximation.

4. Simulation Study and Model Discrimination

An extensive simulation study has been conducted in order to evaluate the performance of the proposed model and the associated estimation method. In the simulation study, three different cluster sizes have been considered to examine the impact of sample size on the proposed model: $n=100$, $n=200$, and $n=400$. Each cluster can be treated as one subject, and the number of observations for one subject was then generated randomly from a normal distribution $N(\mu =8,{\sigma }^2=4)$. 10%, 25% and 50% censoring rates have also been introduced through a Bernoulli distribution (with $\pi =0.9$, $0.75$, and $0.5$). True values of the treatment effects used in the simulation study are chosen to be ${ϵ}_0=2$ and ${ϵ}_1=0.1$. The parameter in the normal distribution to model the cluster effect has been set as $\sqrt{d}=0.5$. Gamma, Weibull and lognormal distribution with mean 1 and variance 0.5 have been used to model over-dispersion in the simulated data sets. We then generated 500 random samples from each scenario, and then fitted the three special distributions along with the generalized gamma distribution of the proposed model to these simulated data sets.

Table 1 presents the bias and mean square errors (MSEs) ($\times 100$) of estimated treatment effect, normal random effect and frailty variance when the data are simulated with 10% censoring rate. Looking at the estimated treatment effect, all models perform well in estimation, even when the frailty distributions differ from the correct one. It is interesting to note that the MSEs are the smallest when Gamma frailty model is being fitted; also, in most cases, the MSEs are the second smallest when the generalized gamma frailty is being fitted. When the number of clusters increase from 100 to 400, we notice a trend of the MSEs tending to become smaller. When examining the estimated normal random effect, models with Weibull and generalized gamma frailty have smaller bias than models with Gamma and Lognormal frailty. The bias and MSEs are the smallest when the generalized gamma model is being fitted. There is also a trend that as the number of clusters increases, the bias and MSEs tend to become smaller. There are noticeable differences in estimating the variance of the frailty distribution among the model considered. When the fitted model matches the true model, the bias and MSEs are quite small. However, it is evident that the lognormal frailty model tends to have very large bias while estimating the variance except for the case when the true frailty distribution is lognormal. This shows that Weibull-Lognormal-Normal Model is not robust.

Model discrimination is motivated by the fact that the generalized gamma distribution encompasses some commonly used distributions as special cases. Choosing the scale and shape parameters from this broad family of distributions suitably, we can adequately fit an appropriate model to a data set. The model discrimination study will allow us to investigate how often a true model gets selected and others get rejected, considering the generality and flexibility of the generalized gamma distribution. The model discrimination is carried out here based on an information-based criterion. In our model selection, we use Akaike’s information criterion (AIC) as the criterion to discriminate among the candidate models. AIC is given by $-2\widehat{l}+2p$, where $\widehat{l}$ is the maximized log-likelihood value and p is the number of model parameters to be estimated. It is to be noted that the model with the minimum value of AIC is the model that best describes the data.

The number of times that the correct model gets chosen and the incorrect models get selected have been computed.

Table 2 presents the selection rates based on AIC for the cases with 10% censoring rate. From Table 2, it can be seen that the selection rate of the correct model ranges from 31.2% to 40.0% if the true distribution is gamma; from 54.4% to 67.8% if the true distribution is lognormal; and from 49.6% to 62.2% if the true distribution is Weibull. It can also be seen that as the sample size increases, the selection rates for the correct model also increases, as expect. We observe that the generalized gamma always has the highest selection rate, which is also close to the correct selection rates. For example, when the true model is gamma and the sample size is 100, the selection rate of generalized gamma is 47%. It tells us that it contains all 31.2% of the time the correct gamma is selected, and the remaining 15.8% of the time some other member of generalized gamma distribution is selected.

If one is interested in testing the validity of one model within the generalized gamma family, likelihood ratio test can be performed. But that will only test for the validity of one model, and would not lead to the problem of selecting the best model.

5. Illustrative Real-Life Data Analyses

In this section, we will illustrate the performance of the proposed model with two real data sets from the survival analysis literature.

5.1. Asthma Study

Asthma is claimed to be the most common chronic condition in children globally, and it is associated with accelerated loss of lung function [17]. The data of interest are the recurrent asthma attacks in children, studied by Duchateau and Janssen [2]. In this study, a clinical trial is conducted with young children with higher risk of developing asthma to examine the effectiveness of a new application of an existing anti-allergic drug. There are 232 such children in the study and followed up to 600 days; they are randomized to either the treatment (new application of the drug) or placebo group (standard application of the drug). A child can develop more than one asthma attack, so intermittent events are ordered in time and clustered within one child. A calendar time format is chosen to represent the data, where the time at-risk for a particular event is the time from the end of the previous event (asthma attack) to the start of the next event (start of the next asthma attack). During the course of the study, the start and end dates of different at-risk periods are recorded for one patient, and these at-risk periods are separated either by an asthmatic event or by a period in which the patient is not under observation. The goal of the study is then to investigate whether the new application of the anti-allergic drug has an effect on reducing the frequency and length of the asthma attacks in children. In Table 3, the data for the first two patients are presented, for example.

In this section, various frailty distributions under the proposed generalized gamma frailty and normal random effects framework are fitted for the time-to-event data introduced above. Our aim is to determine the impact of different frailty distributions on estimating the treatment effects, and to find a model that provides the best fit to these data.

We first consider an exponential model by setting $\rho =1$ in the model of the form in (18). The predictor has a form of ${\kappa }_{ij}={\xi }_0+{\xi }_1{T}_i+b_i$, where $T_i$ takes on the value 1 if the patient is in the treatment group and 0 if in the control group, and $b_i\sim N(0,d)$. We fit the data with the generalized gamma frailty and normal random effect model and its three special cases: gamma, lognormal and Weibull frailty and normal random effect models. The obtained results are reported in Table 4. Wald tests have been carried out to assess the estimated treatment effects and the corresponding results are presented in Table 5. It is seen that the treatment effect ${\xi }_1$ are reasonably close to each other. The generalized gamma, which contains lognormal, has the smallest AIC. This suggests that using generalized gamma can aviod the use of the wrong model as the generalized gamma could direct us effectively to the correct model. From Table 5, the p-values appear to be around 0.14, which suggests that the treatment effect is identifiable in all of the four fitted models.

Next, we take censoring into consideration when fitting all four models to the asthma data. The corresponding results are presented Table 6. In both Table 4 and Table 6, it is important to note that the standard errors of the two parameters, ${\xi }_0$ and $\lambda$, are quite large. This could be the result of overdispersion being present in the data. In order to address this issue, it would make sense to set the scale parameter $\lambda$ to be 1. We then re-fit all four models, and the corresponding results are shown in Table 7 and Table 8, for data without and with censoring, respectively.

Next, we take censoring into consideration when fitting all four models to the asthma data. The corresponding results are presented Table 6. In both Table 4 and Table 6, it is important to note that the standard errors of the two parameters, ${\xi }_0$ and $\lambda$, are quite large. This could be the result of overdispersion being present in the data. In order to address this issue, it would make sense to set the scale parameter $\lambda$ to be 1. We then re-fit all four models, and the corresponding results are shown in Table 7 and Table 8, for data without and with censoring, respectively.

From Table 7 and Table 8, it is evident that the standard errors are now lower and reasonable. Wald tests have again been performed to test the significance of the estimated treatment effects from all of the models based on the results in Table 7 and Table 8. Wald tests suggest that the treatment effects are still identifiable in models both with and without censoring. From Table 9, we observe that estimated treatment effects and the standard deviation of the normal random effects are similar in the four models.

5.2. Bladder Study

The model proposed here is further illustrated by the bladder dataset, which is available from survival package in R [18,19]. The study was conducted by the Veterans Administration Cooperative Urological Rsearch Group. This study involved patients with superficial bladder tumors at the start of the trial. Following the transurethral removal of tumors, the patients were randomly treated by placebo, thiotepa or pyridoxine. Most patients experienced repeated tumor recurrences and new tumors were removed at each follow-up. The recurrence times of tumor were recorded. The subset of 85 patients who were assigned to either thiotepa or placebo has been used for analysis here. The exponential model taking $\rho =1$ with generalized gamma frailty and its special cases are fitted with the data. The corresponding results are presented in Table 10. It is observed that the model with generalized gamma frailty has the smallest AIC, which suggests once again the usefulness of the generalized gamma frailty in picking the best model in a flexible and efficient way.

6. Discussion and Conclusions

Building on the work of Molenberghs et al. [1], in this article, we have extended the combined model of gamma frailty and normal random-effects for repeated survival data by using generalized gamma frailty distribution. Particular attention has been given to Weibull models for overdispersed, clustered time-to-event outcomes, with generalized gamma and normal random effects. Censoring has also been accommodated in the data. Model selection and goodness-of-fit are critical issues in frailty models for time-to-event data. Due to the latent nature of frailties in modelling, it is often challenging to determine an appropriate frailty distribution for a given data set. Balakrishnan and Peng [7] proposed generalized gamma distribution in a frailty model context thereby significantly improving the fit of the frailty model to the observed data. The generalized gamma distribution is a flexible frailty distribution with two parameters besides the scale parameter, and it includes gamma, lognormal and Weibull as special cases. The integration of the generalized gamma distribution to the combined model can therefore provide more flexibility for modelling clustered time-to-event data.

The combined Weibull-gamma-normal model proposed by Molenberghs et al. [1] enjoys the so-called strong conjugacy property, and closed-form expressions for the marginal distribution can be derived conveniently and then made use of subsequently in the maximum likelihood estimation. The computational challenge faced in considering the generalized gamma distribution as the frailty distribution is handled conveniently by using Monte Carlo approximation for the integrals involved in the likelihood function. The performance of the proposed model has been evaluated through an extensive simulation study. The estimates of treatment effect and normal random effect coefficients are seen to be accurate. As the sample size increases, the bias and MSEs of the estimates decrease. In addition, we observe that a lower censoring rate leads to a more accurate estimation of the parameters in the proposed model. A model discrimination is also carried out between Weibull-generalized gamma-normal model and its three special cases: Weibull-gamma-normal, Weibull-lognormal-normal, and Weibull-Weibull-normal model. Data sets are simulated from the four models, and the four candidate models are then fitted to each data set. The model selection has been done based on AIC or the value of the maximized log-likelihood function. From our analysis, we observe that the selection rates of the correct model increase as the sample size increases. Finally, a recurrent asthma attack in children data set and a recurrent bladder data set are used to illustrate the proposed Weibull-generalized gamma-normal random effect model. The proposed model is sufficient and flexible enough to provide the best fit for the recurrent asthma data set and the recurrent bladder data set.

In our model formulation, the generalized gamma and normal random effects capture overdispersion and between-subject association. Although the proposed model is flexible and widely applicable for univariate time-to-event outcomes, it is possible to generalize to bivariate or multilevel repeated time-to-event settings. Furthermore, regarding censoring, we have confined the analysis to right-censored outcomes. However, the methoodology can be extended with no trouble to left-censoring and interval censoring. All computations have been implemented in R. Readers who are interested in the proposed model can request a copy of all datasets and software codes from the authors.