Multi-Parameter Regression Survival Modelling with Random Effects

by   Fatima-Zahra Jaouimaa, et al.

We consider a parametric modelling approach for survival data where covariates are allowed to enter the model through multiple distributional parameters, i.e., scale and shape. This is in contrast with the standard convention of having a single covariate-dependent parameter, typically the scale. Taking what is referred to as a multi-parameter regression (MPR) approach to modelling has been shown to produce flexible and robust models with relatively low model complexity cost. However, it is very common to have clustered data arising from survival analysis studies, and this is something that is under developed in the MPR context. The purpose of this article is to extend MPR models to handle multivariate survival data by introducing random effects in both the scale and the shape regression components. We consider a variety of possible dependence structures for these random effects (independent, shared, and correlated), and estimation proceeds using a h-likelihood approach. The performance of our estimation procedure is investigated by a way of an extensive simulation study, and the merits of our modelling approach are illustrated through applications to two real data examples, a lung cancer dataset and a bladder cancer dataset.



There are no comments yet.


page 1

page 2

page 3

page 4


A Flexible Parametric Modelling Framework for Survival Analysis

We introduce a general, flexible, parametric survival modelling framewor...

Semiparametric multi-parameter regression survival modelling

We consider a log-linear model for survival data, where both the locatio...

Yakovlev Promotion Time Cure Model with Local Polynomial Estimation

In modeling survival data with a cure fraction, flexible modeling of cov...

Multi-Parameter Regression Survival Modelling: An Alternative to Proportional Hazards

It is standard practice for covariates to enter a parametric model throu...

A Unifying Framework for Flexible Excess Hazard Modeling with Applications in Cancer Epidemiology

Excess hazard modeling is one of the main tools in population-based canc...

Convergent stochastic algorithm for parameter estimation in frailty models using integrated partial likelihood

Frailty models are often the model of choice for heterogeneous survival ...

A Multi-parameter regression model for interval censored survival data

We develop flexible multi-parameter regression survival models for inter...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Standard survival analysis models such as the proportional hazards (PH) model are known to have a single covariate-dependent parameter, the scale parameter. As a means to extend these models and afford more flexibility in the modelling of survival data, Burke and MacKenzie (2017) developed a multi-parameter regression (MPR) approach for a parametric hazard function. MPR is an approach whereby more than one distributional parameter is allowed to depend on covariates, and this is sometimes referred to as “distributional regression” (Stasinopoulos et al., 2018); see also Rigby and Stasinopoulos (2005) and Stasinopoulos et al. (2007). Using the Weibull MPR model as an example, Burke and Mackenzie’s paper demonstrates the advantages of allowing both the scale and the shape parameters to depend on covariates simultaneously. More recently, Burke et al. (2019) extended a semi-parametric accelerated failure time (AFT) model to multiparameter regression and Burke et al. (2020) explored an MPR parametric survival modelling framework in which the baseline hazard function follows an adapted power generalized Weibull (APGW) distribution.

Standard MPR models rely on the assumption of independent event times, and such an assumption cannot be met in studies with clustered observations. This includes studies where successive or recurrent event times are recorded on each subject; multi-centre studies where survival times of individuals from the same centre may be dependent due to centre-specific conditions, clinical or otherwise; family studies; and matched pair studies. Various methods have been developed to model the lack of independence in clustered data, but perhaps the most standard approach is the one whereby a cluster-specific random effect is introduced into the model. Given these random effects, the data are assumed to be conditionally independent. In the context of survival analysis, such random effects models are commonly referred to as frailty models and since the random effect represents a common effect on all members in a cluster, these types of common risk models are known as shared frailty models (Clayton, 1978; Duchateau and Janssen, 2008; Hougaard, 2012).

Extensions of the MPR models to include random effects have been centred around the “classical” multiplicative frailty model, whereby the frailty term is included in the linear predictor corresponding to the scale parameter. Peng et al. (2020)

extend a Weibull MPR model for interval censored data to include a gamma distributed multiplicative frailty term, but only for univariate frailty. In their paper,

Jones et al. (2020)

extend the MPR PGW and MPR APGW models to handle bivariate data using multiplicative frailty, where the dependence is understood using links to the well known power variance copula

(Goethals et al., 2008; Duchateau and Janssen, 2008; Hougaard, 2012). The attraction of MPR modelling is the flexibility afforded by allowing the hazard shapes to differ, but just as the shape can vary according to covariates, so too could it have random variation, for example, resulting from differences across multiple centres. Hence, in this paper, we go beyond the classical multiplicative model to allow for a wider range of model structures, specifically, models in which a frailty term is included in each distributional parameter. This, we believe, is a more natural way of modelling correlated data using MPR models than multiplicative frailty which relates only to the scale of the hazard. Furthermore, our proposed model accounts for the potential correlation between the frailty terms themselves since scale and shape effects may be correlated.

Parameter estimation in the parametric shared frailty model is commonly achieved through integrating out the frailties from the conditional survival likelihood. The resulting equation is an explicit expression for the marginal likelihood, containing the fixed parameters and no longer the frailties. This marginal likelihood can then be maximized using numerical methods to yield estimates for the fixed effects. Inference about the random effects is not readily available and integrating out the frailties from the joint density typically involves the evaluation of analytically intractable integrals over the random-effect distributions. Numerical integration methods such as the Gauss-Hermite (G-H) quadrature can be used to approximate the value of the integrals, however, when the dimensionality of the integral is high, the number of quadrature points grows exponentially with the number of random effects, and the approximation is sub-optimal. While several methods, such as the Monte Carlo Expectation-Maximisation, Markov Chain Monte Carlo and Gibbs sampling have been used to overcome the issue of intractable integrals, these methods are notoriously known for being computationally intensive. This is especially true when the number of random effects (clusters) is large or when the complex correlation structure among the clustered survival times requires the assumption of multiple frailties

(Vaida and Xu, 2000; Ripatti and Palmgren, 2000; Abrahantes et al., 2007; Duchateau and Janssen, 2008). We adopt a so-called “hierarchical likelihood” (h-likelihood) approach. This approach was originally proposed by Lee and Nelder (1996) for a generalized mixed-effects model but further studied and developed to provide a straightforward, unified framework for various random effect models including frailty models (Ha et al., 2001, 2002; Ha and Lee, 2003; Ha et al., 2017).

In contrast to the standard marginal likelihood approach, the h-likelihood framework treats the random effects or frailties as model parameters, which are then jointly estimated with the fixed parameters and the frailty dispersion parameter(s). Estimates of the fixed and random effects are found by maximising the log-likelihood function conditional on the random effects plus a penalty term whose value depends on how dispersed the random parameters are, i.e., if the random effects have a large dispersion parameter, then the penalty term takes a large value. By treating the random effects as model parameters, the h-likelihood framework avoids the intractable integration needed to calculate the marginal likelihood and provides an efficient estimation procedure. Furthermore, classical analysis of random effects models focuses on the estimation of the fixed parameters and the frailty variance parameter(s), but in many recent applications, estimation of the random effects is also of interest. Such estimates allow for the survivor function for individuals with given characteristics to be estimated, i.e., the cluster specific failure time distribution, and are especially useful in multi-centre studies when the frailty term represents the centres. Estimates of the random effects in such studies provide information about the merits of the different centres in terms of patient survival and are useful for investigating the potential heterogeneity in survival among clusters in order to better understand and interpret the variability in the data (Ha et al., 2016).

The remainder of this article is organized as follows: Section 2 describes the MPR model, its extended version with random effects and the h-likelihood procedure for parameter estimation in the given model. Section 3 presents results of extensive simulation studies. The proposed methods are illustrated on two datasets arising from multi-centre studies and the results are shown in Section 4, followed by a discussion and conclusion in Section 5.

2 The MPR Frailty model

2.1 Model Formulation

To formulate a shared frailty MPR model for survival data, we assume time-to-event data arising from a multi-centre study. In this study, we have centres or clusters, with individuals (patients), . The total sample size is the total number of individuals coming from all centres, i.e., . We define as the survival time for the th individual, , in the th cluster and as the corresponding censored time. The cumulative hazard function for a shared frailty MPR model takes the following parameteric form:


where is the underlying cumulative hazard function with scale parameters and shape parameter . The corresponding hazard function is given by

where is the baseline hazard function. can be one of the commonly used survival distributions (Weibull, Gomperts, or log-logistic), see Table 1.

Table 1: Possible distributions.

The two distributional parameters, and depend on covariates as follows:


is the covariate vector,

and are the corresponding regression coefficient vectors, and denote the scale and shape random effects from the th cluster respectively, and to ensure positivity of the parameters a log link is used. Although we allow the scale and the shape parameters to depend on the same set of covariates, the parameters may or may not have covariates in common depending on the value of the corresponding regression coefficients. The individual cluster-specific effects (scale or shape random effects) are assumed to be independently and identically distributed (i.i.d.) according to some distribution.

In their work, Burke and MacKenzie (2017) found the distributional parameters in an MPR model to be highly correlated. To account for the possibility of the propagation of this correlation to the corresponding frailty terms, we allow for correlation between the two random effects terms, and , by assuming that

follow a bivariate normal distribution such that

Fitting the model which assumes the bivariate normal distribution for the frailty terms avoids the assumption of independence which may be unnecessarily restrictive. Furthermore, this model generalizes various models and simplifications of it by assuming null components or a specific correlation structure include:

  • the proportional hazards (PH) model (Cox, 1972);

  • a multiplicative frailty PH model (Duchateau and Janssen, 2008);

  • a standard MPR model (Burke and MacKenzie, 2017; Burke et al., 2020);

  • a multiplicative scale frailty MPR model (Peng et al., 2020; Jones et al., 2020);

  • a shape frailty MPR model;

  • : we assume that the two frailty terms are independent and fit a model whereby and ;

  • : and , where is a real-valued scaling factor (Ha et al., 2017).

Figure 1 provides a summary of such models. Although, we only consider the normal distribution for the random effects, estimates of the fixed effects are usually robust against violations of this assumption if the censoring rate or frailty variance parameter is not too high (Ha et al., 2001; Ha and Lee, 2003; Ha et al., 2011, 2016).

Figure 1: A schematic diagram of some of the possible models generalized by the MPR frailty model assuming a bivariate normal distribution. Note that .

2.2 Construction of the H-Likelihood

Denoting the observed data by the pairs (), where , the observed survival time for the th individual in the th cluster; and is the censoring indicator, which takes the value 0 for a censored observation and 1 for an event. Under the standard assumptions that, given and , the censored times (s) and the event times (s) are conditionally independent and the censoring is conditionally non-informative, the log-h-likelihood function of the proposed model (1) is given by



is the logarithm of the conditional joint density function for and given the random effects , where is the realisation of and

is the logarithm of the density function of with dispersion parameters (i.e. frailty parameters) , and is the vector of the fixed parameters. For meaningful inferences, it is important to define the h-likelihood on a particular scale of , such that the random effects occur linearly in the linear predictor. Hence, the log link we see in (1). (Lee and Nelder, 1996; Ha et al., 2001; Lee et al., 2017; Ha et al., 2017).

2.3 Estimation Procedure

From the h-likelihood function given in (2), we can derive two likelihoods, namely the marginal likelihood and the restricted likelihood. The marginal likelihood eliminates the random effects, , from h while the restricted likelihood eliminates the fixed effects from the marginal likelihood, i.e., eliminates both the fixed effects and the random effects from h. In theory, the h-likelihood should be used for inference about , the marginal likelihood should be used for inference about and the restricted likelihood for inference on the dispersion parameter(s) (Patterson and Thompson, 1971; Harville, 1977). However, when the marginal likelihood is intractable, Lee and Nelder (1996, 2001) propose using adjusted profile likelihoods, and as approximations to the marginal likelihood and the restricted likelihood respectively. Consider the log-h-likelihood h with nuisance parameters , the adjusted profile likelihood suggested by Lee and Nelder is given by

where and solves . The function produces an adjusted profile likelihood, profiling out the nuisance parameters . The nuisance parameters can be the fixed effects and/or the random effects. is the first-order Laplace approximation to the marginal likelihood, (Lee and Nelder, 2001). Similarly, is the first-order Laplace approximation to the restricted likelihood. Moreover, approximates and therefore, . With the exception of binary data with a small cluster size (i.e., ), it has been found that h and give very similar results when used in the estimation of the fixed effects (Lee et al., 2017). Hence, Lee et al. (2017) recommend the joint maximisation of over the fixed and random effects and refer to this as the uncorrected h-likelihood method. This methods works well for various models (Ha et al., 2002, 2017), including the models we propose.

The h-likelihood function given in (2) is maximized to obtain the maximum h-likelihood estimators (MHLEs) of and . The score functions are given by

is an matrix whose th row is , is an matrix whose th row is , a vector indicating the cluster effect; and are vectors of length such that

and and are vectors of length , such that

The corresponding observed information matrix can be written explicitly as


where and are and model matrices for , and whose th rows are and , respectively. , and are diagonal matrices whose th diagonal elements are given by

where . , and are matrices whose th elements arise from such that

where is a identity matrix. Given , the following system of Newton-Raphson equations can be solved iteratively for the MHLEs of and


where and . For the estimation of the dispersion parameters, , the adjusted profile likelihood, , which eliminates is used. Given , , given in (2) and defined in (3), is defined as follows:


Solving the equations yields the restricted maximum likelihood estimators (REMLEs) of . We opt for a non-linear optimizer implemented in R using the function nlm. The procedure iterates between and

until all the estimates converge. The standard errors for

and can be estimated directly from the inverse of the observed information matrices, and respectively (Ha et al., 2016, 2017).

2.4 Fitting Algorithm

The model estimation algorithm described above can be summarized as follows:

Initialisation: Using 0.01 as the initial values for the scale and shape coefficients, we fit a fixed effects Weibull MPR model. Estimates from this model are used as the initial values for the fixed parameters, , in the mixed effects Weibull MPR model. For the initial values of the random effects, , we use 0.01 and (0.1, 0.1, 0.1) is used for the dispersion parameters, .

Parameter Estimation:

Step 1 Keeping the frailty variance parameters fixed, maximize by iteratively re-solving the system of equations given in (4) to obtain the new estimates .

Step 2 Given the estimates from Step 1, a new estimate is obtained by maximising (5) using nlm.

Iterate between Step 1 and Step 2 until the convergence criterion is met, that is until the maximum absolute difference between the previous and current estimates for and is less than . After convergence is reached, the standard errors are estimated.

From the simulations we have carried out, we have found this algorithm to be time efficient and insensitive to starting values.

3 Simulation Studies

The performance of the proposed methods is evaluated through simulation studies. The Weibull distribution is one of the most commonly used distributions in survival analysis, hence, we chose to generate the survival times from a Weibull MPR frailty model with the following regression parameters

The corresponding covariates, , were generated from an AR(1) process with a correlation coefficient of

. Each variable is marginally standard normal. The corresponding censored times were generated from a uniform distribution with a censoring rate of approximately

or respectively. Following the two real data structures in Section 4, three different cluster sizes, , and two different cluster numbers, , are considered (note that in contrast to this scheme, the real data has varying cluster sizes and this is something we have investigated, see Appendix). The frailty terms are generated from a bivariate normal distribution with the combination of various dispersion parameter values , , . Each simulation scenario was replicated 500 times.

To summarize the simulation results, we compute the mean, the standard deviation (SE) and the average standard errors (SEE) for all the parameters we estimated for each scenario. The results for a censoring rate of

, , and are presented in Table 2. Similar tables of the results from other simulation set-ups can be found in the Appendix. Overall, the h-likelihood estimates of both the fixed parameters and the frailty dispersion parameters perform quite well. The bias in the estimates is reduced as we increase both the cluster size and the number of clusters, this is observed in all the combinations of dispersion parameters and for both censoring percentages considered. The standard errors appear to be underestimated in the smaller sample sizes. This is especially true for the frailty variance parameters, and to a lesser extent, for the fixed effects. As we increase the cluster size and the number of clusters, however, we see that the standard errors shrink for all the parameters, and the SEE converges towards the SE; in any case, the use of the standard errors in hypothesis testing for the frailty variance parameters is recommended against (Maller and Zhou, 2003). Similar but larger bias is observed in the results of scenarios with a censoring (see Appendix).

Mean Mean Mean Mean Mean Mean Mean Mean Mean
(SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE)
() (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE)
True 1 -0.5 0.5 0.5 0.5 -0.5 1 0.5 -0.5
(20, 5) 1.30 -0.58 0.59 0.64 0.50 -0.49 1.18 0.51 -0.41
(0.35) (0.27) (0.25) (0.18) (0.11) (0.11) (0.33) (0.15) (0.39)
(0.32) (0.21) (0.20) (0.15) (0.09) (0.09) (0.20) (0.09) (0.19)

(20, 20)
1.05 -0.52 0.52 0.54 0.50 -0.50 1.00 0.50 -0.50
(0.24) (0.08) (0.08) (0.12) (0.04) (0.04) (0.20) (0.09) (0.23)
(0.24) (0.08) (0.08) (0.12) (0.04) (0.04) (0.16) (0.08) (0.17)

(20, 50)
1.03 -0.51 0.50 0.52 0.50 -0.50 1.01 0.50 -0.49
(0.23) (0.05) (0.05) (0.12) (0.02) (0.02) (0.17) (0.08) (0.19)
(0.23) (0.05) (0.05) (0.12) (0.02) (0.02) (0.16) (0.08) (0.17)

(100, 5)
1.21 -0.56 0.57 0.59 0.49 -0.49 1.07 0.52 -0.46
(0.15) (0.10) (0.10) (0.08) (0.04) (0.04) (0.13) (0.06) (0.16)
(0.13) (0.08) (0.08) (0.07) (0.04) (0.04) (0.08) (0.04) (0.09)

(100, 20)
1.05 -0.51 0.51 0.53 0.50 -0.50 1.01 0.50 -0.50
(0.11) (0.04) (0.04) (0.05) (0.02) (0.02) (0.09) (0.04) (0.10)
(0.11) (0.04) (0.04) (0.05) (0.02) (0.02) (0.07) (0.04) (0.08)

(100, 50)
1.02 -0.50 0.50 0.51 0.50 -0.50 1.00 0.50 -0.50
(0.11) (0.02) (0.02) (0.05) (0.01) (0.01) (0.08) (0.04) (0.08)
(0.10) (0.02) (0.02) (0.05) (0.01) (0.01) (0.07) (0.04) (0.08)
Table 2: Averaged coefficient estimates, standard deviations (SE) and the average standard errors (SEE) for the simulation scenario with dispersion parameters , and and a censoring rate of .

4 Data Analysis

4.1 Modelling Details

We also illustrate the proposed models on two datasets. Though our main focus is the MPR model with BVN frailties, for the purpose of comparison and analsyis we also consider various simplifications of this model. More precisely, we fit Weibull MPR models with the following frailty structures to each of the datasets we consider:

  • BVNF: BVN frailty ()

  • IF: independent frailty ()

  • CF: common frailty ()

  • ScF: scale frailty only (), i.e., the classical multiplicative frailty model

  • ShF: shape frailty only ()

These models are fitted following the procedures described in Sections 2.3 and 2.4, and the standard errors for the estimated parameters are computed as described in Section 2.3. For the selection of the fraily structure that is best supported by the data, we use the Akaike information criterion (AIC). Various extended definitions of the AIC in random-effect models can be formulated based on different likelihood functions (Vaida and Blanchard, 2005; Xu et al., 2009; Ha et al., 2007, 2017). More specifically, we make use of the restricted AIC (rAIC) (Ha et al., 2007) and the conditional AIC (cAIC) (Vaida and Blanchard, 2005; Ha et al., 2017). The rAIC is based on the restricted likelihood approximation , which eliminates from and thus is a function of the frailty parameters only. While the cAIC is based on the conditional joint density function for and given the random effects, . Definitions of the AICs are as follows:

where is the number of dispersion parameters governing the frailty distribution,

, the effective degrees of freedom adjusted for the fixed and random effect estimates, and

(Ha et al., 2017; Lee et al., 2017). The computation of involves the fixed effects, the random effects and the frailty distribution parameters, but in a model with no frailty, is just the number of fixed effects in the model, and hence the cAIC becomes the classical AIC in this case (Ha et al., 2017). After fitting each of the aforementioned models as well as a model with no frailty (NF), we obtain the corresponding rAIC and cAIC values for the purpose of model comparison.

In the Weibull MPR models (or any other parametric MPR model with a scale and shape parameter), the scale coefficients describe the overall scale of the hazard and the shape coefficients describe its evolution over time. A positive scale coefficient indicates an increase in the hazard relative to some reference category and similarly, a positive shape coefficient indicates an increased hazard over time relative to some reference category. While an examination of the and coefficients separately provides some initial understanding of the effect of a variable; it is important to look at the combined information from both coefficients when determining its overall effect, and hence, it is important to look at the hazard ratios. For a binary covariate, , Burke and MacKenzie (2017) show the hazard ratio under the Weibull MPR frailty model is given by

where and are, respectively, the scale and shape coefficients of , and , the covariate vector with set to 0. (Note here that we have dropped the subscripts for notational convenience.) Because depends on the values of the other covariates in the model via the vector , we set them to their empirical modal values. In line with this, we also set the random effect to its modal value of zero.

4.2 Extensive-Stage Small-Cell Lung Cancer

This dataset was collected as part of a randomized, multi-centre study conducted by the Eastern Cooperative Oncology Group (ECOG). The main purpose of the trial was to determine if cyclic alternating combination chemotherapy was superior to cyclic standard chemotherapy in patients with extensive-stage small-cell lung cancer. Patients were randomly assigned to one of two treatment arms, standard chemotherapy (CAV; reference category) or an alternating regimen (CAV-HEM). The dataset includes 579 patients from 31 different institutions, with the number of patients per institution ranging from 1 to 56 and a median of 17 patients per institution. The outcome variable was time (in years) from randomisation until death. The median survival time and maximum follow-up time were 0.86 years and 8.45 years, respectively, and of the 579 study participants only 10 were censored yielding a censoring rate of approximately 1.7%. Besides the survival time, censoring status, institution code and treatment, four other dichotomous variables were included in this dataset, namely (reference category listed first): the presence of bone metastases (no, yes), the presence of liver metastases (no, yes), patient status on entry (confined to bed or chair, ambulatory), and whether there was a weight loss prior to entry (no, yes). More details on the trial and it’s clinical results can be found in Ettinger et al. (1990). This dataset was also previously analysed in Gray (1994) using a fully Bayesian approach, in Vaida and Xu (2000) using a marginal likelihood approach, and in Ha et al. (2016) using a correlated frailty model fitted using a h-likelihood approach.

We fit the models listed in Section 4.1 to this lung cancer dataset and the results are presented in Table 3. To explore the degree of dependence between the two random components, and , we first fit the BVNF model, (i.e., the model which assumes a bivariate normal distribution). The estimate of () indicates a very strong positive correlation between the predicted random components and . This perhaps suggests that the model could be simplified and only one random component is needed. Note that the large value is due to the very small values, and this may be pointing towards a shape frailty only model.

Scale Intercept 0.11 0.12 0.11 0.11 0.12 0.11
(0.13) (0.13) (0.13) (0.13) (0.13) (0.13)
Treatment -0.24 -0.23 -0.24 -0.24 -0.24 -0.24
  CAV-HEM (0.09) (0.09) (0.09) (0.09) (0.09) (0.09)
Bone Metastases 0.23 0.24 0.22 0.22 0.25 0.22
  Yes (0.10) (0.10) (0.10) (0.10) (0.10) (0.10)
Liver Metastases 0.33 0.38 0.39 0.39 0.32 0.39
  Yes (0.10) (0.10) (0.10) (0.10) (0.10) (0.10)
Patient Status -0.58 -0.63 -0.62 -0.62 -0.60 -0.62
  Ambulatory (0.11) (0.11) (0.11) (0.11) (0.11) (0.11)
Weight Loss 0.19 0.21 0.21 0.21 0.19 0.21
  Yes (0.10) (0.10) (0.10) (0.10) ( 0.10) (0.10)
Shape Intercept 0.14 0.23 0.22 0.22 0.16 0.22
(0.08) (0.12) (0.12) (0.12) (0.08) (0.12)
Treatment -0.28 -0.23 -0.24 -0.24 -0.27 -0.24
  CAV-HEM (0.06) (0.07) (0.07) (0.07) (0.06) (0.07)
Bone Metastases 0.04 0.02 0.03 0.03 0.03 0.03
  Yes (0.07) (0.08) (0.08) (0.08) (0.07) (0.08)
Liver Metastases -0.14 -0.11 -0.11 -0.11 -0.13 -0.11
  Yes (0.07) (0.07) (0.07) (0.07) (0.07) (0.07)
Patient Status 0.33 0.38 0.38 0.38 0.33 0.38
  Ambulatory (0.07) (0.09) (0.09) (0.09) (0.08) (0.09)
Weight Loss 0.06 0.02 0.03 0.03 0.05 0.03
  Yes (0.07) (0.08) (0.08) (0.08) (0.07) (0.08)
Frailty       parameters 0.08 0.00 0.00 0.16
(0.01) (0.02) (0.00) (0.03)
0.29 0.31 0.31
(0.04) (0.05) (0.05)
1123.40 1077.74 1079.52 1079.52 1121.35 1079.52
0 3 2 2 1 1
41.88 2.22 2.00 2.00 41.83 0
12.00 30.96 30.89 30.89 19.89 30.89
65.19 0 2.08 2.06 61.09 2.08
  • Note that for the CF model, an estimate for can be found by evaluating . Given the above values .

Table 3: The coefficient estimates, frailty dispersion parameter estimates and estimated standard errors (in brackets) from each model we fit to the lung cancer dataset.

With the exception of the ScF model, all the models containing at least one frailty component have quite similar rAIC values, all much smaller than the rAIC value from a no frailty model, with the model with a frailty term in the shape only, the ShF model, having the lowest value. The ScF model, the more standard multiplicative frailty model, has a similar rAIC value to the model with no frailty components, suggesting that the baseline risk is homogenous across the centres and a frailty term in the scale is not needed. We observe a similar pattern in the cAIC values, in that all the models containing at least one frailty component (with the exception of the ScF model) have similar cAIC values and all much smaller compared to the NF model. Although the BVNF model has the lowest cAIC value, it is still close to those of the IF, CF, and ShF models. Given this, and also the rAIC values, we consider the ShF as the “best” model.

Although we report the standard errors of the frailty variance parameters, one should not use them for testing the hypothesis (Vaida and Xu, 2000). A likelihood ratio test can be carried out instead. Since the value of

in the null hypothesis is at the boundary of the parameter space, the standard approximation of the loglikelihood ratio statistic by a

distribution often leads to over conservative test results (Self and Liang, 1987; Stram and Lee, 1994; Duchateau and Janssen, 2008)

. To correct this bias, a mixture of a chi-square distribution with one and zero degrees of freedom,

, should be used as the approximation of the log-likelihood ratio statistic (Maller and Zhou, 2003)

. The test statistic at the 5% significance level is thus 2.71. The difference in deviance,

between the NF model and the ShF model is 43.878, and hence the scale frailty is significant, i.e., .

Focusing on the results from the model selected by the rAIC, the ShF model; considering the scale parameter results first, all the variables included in the model have significant scale coefficients. The CAV-HEM treatment and the subject being ambulatory on entry have negative scale coefficients and so reduce the hazard of death relative to their respective reference categories, the CAV treatment and subject being confined to bed or chair on entry respectively. The presence of bone metastases, the presence of liver metastases and weight loss prior to study entry are all found to increase the hazard of death relative to their reference categories. Now, considering the shape parameter, only two variables, treatment and patient status on entry have significant coefficients. The CAV-HEM treatment has a negative shape coefficient suggesting that the hazard further decreases over time, relative to the reference category. The positive shape coefficient for the subject being ambulatory on entry suggests that the hazard increases over time, although this variable has an initial effect of reducing the hazard of death, this effect wears off over time, relative to its reference category.

Figure 2

shows the hazard ratio corresponding to each of the variables in our model along with 95% confidence intervals estimated using a parametric bootstrap

(Davison and Hinkley, 1997). The presence of bone metastases and weight loss prior to study entry appear to have a more or less constant hazard over time, this is expected since their corresponding shape effects are quite small and not significant. The presence of liver metastases has a negative effect on the hazard but this effect wears off within the first 2 to 3 years. The hazard ratio for the patient being ambulatory on entry appears to be increasing over time. The effectiveness of the CAV-HEM treatment is only observed after 2 months or so from the treatment start date and the hazard continues to decrease over time relative to the CAV treatment.

The random effects along with their 95% confidence intervals under the ShF model are shown in Figure 3. Note that as the cluster size increases, the confidence bounds around the cluster effect shrink. The biggest changes in centre specific hazard over time can be seen in centres 12, 16 and 20. A positive frailty suggests the hazard is increasing over time relative to the baseline, while a negative one suggests it’s decreasing relative to the baseline.

Figure 2: The hazard ratios with 95% confidence intervals based on results from the ShF model. The modal values were used in the computation of these hazard ratios (Treatment = CAV, Bone Metastases = no, Liver Metastases = no, Patient Status = ambulatory on entry, Weight Loss = yes).
Figure 3: The random effects of the 31 centres with 95% confidence intervals under the ShF model. Centres are sorted in increasing order based on the number of patients.

4.3 Bladder Cancer

This multi-centre dataset was collected as part of the European Organisation for Research and Treatment of Cancer (EORTC) trial 30791 (Sylvester et al., 2006). A total of 410 patients with superficial bladder cancer were included in this dataset. The patients came from 21 different centres and the number of patients per centre varied between 3 and 78 patients with a median of 15 patients per centre. The outcome variable is relapse-free or disease-free interval after transurethral resection, i.e., time from randomisation until cancer relapse. Patients which did not experience a recurrence during the follow-up period were censored at their last date of follow-up. The maximum follow-up time was 10.16 years and of the 410 patients, 204 (approximately 50% of the patients) were right-censored. The two covariates included in this dataset are (reference categories are listed first): a treatment indicator for chemotherapy (no, yes), and a variable representing the prior recurrence (no, yes). This dataset was also previously analysed in Ha et al. (2011) and can be found in the R package frailtyHL (Ha et al., 2012, 2017). As in the previous example, we fit the models listed in Section 4.1 to this dataset and compare the fit using the rAIC and cAIC. The results are presented in Table 4.

Scale Intercept -0.79 -0.71 -0.70 -0.70 -0.70 -0.79
(0.18) (0.19) (0.20) (0.20) (0.20) (0.18)
Chemotherapy -0.72 -0.74 -0.74 -0.74 -0.74 -0.72
  Yes (0.19) (0.19) (0.19) (0.19) (0.19) (0.19)
Prior Recurrence 0.55 0.57 0.57 0.57 0.57 0.55
  Yes (0.17) (0.17) (0.17) (0.17) (0.17) (0.17)
Shape Intercept -0.19 -0.17 -0.19 -0.19 -0.19 -0.19
(0.13) (0.13) (0.13) (0.13) (0.13) (0.13)
Chemotherapy 0.03 0.02 0.03 0.03 0.03 0.03
  Yes (0.13) (0.13) (0.13) (0.13) (0.13) (0.13)
Prior Recurrence -0.01 0.01 0.02 0.02 0.02 -0.01
  Yes (0.12) (0.12) (0.12) (0.12) (0.12) (0.12)
Frailty       parameters 0.22 0.28 0.27 0.28
(0.06) (0.06) (0.06) (0.06)
0.06 0.00 0.03
(0.02) (0.03) (0.03)
946.96 943.75 943.28 943.28 943.28 946.96
0 3 2 2 1 1
1.68 4.47 2.00 2.00 0 3.68
6.00 12.76 13.09 13.11 13.09 6.35
6.46 0.70 0.05 0 0.05 6.45
  • Note that for the CF model, an estimate for can be found by evaluating . Given the above values .

Table 4: The coefficient estimates, frailty dispersion parameter estimates and estimated standard errors (in brackets) from each model we fit to the bladder cancer dataset.

Similar to the previous example, the correlation coefficient between the two random effects, and is suggesting that the model can be simplified. The model with the single frailty parameter in the scale has the lowest rAIC value and therefore is the preferred model, suggesting that there is substantial variation in the baseline risk across the centres and this variation is constant overtime. All the models containing a scale frailty term have similar cAIC values, which are significantly lower than the values for the NF model and the ShF model, suggesting the need for a scale frailty term. Since the difference in cAIC between the models BVNF, IF, CF and ScF is small, we focus on the simplest model in that set, the ScF model.

The difference in deviance, , between the NF model and the ScF model is (), and hence the scale centre effect is significant at the significance level, i.e., . A caterpillar plot of the random effects from the ScF model is presented in Figure 4 and again the centres are sorted by the number of patients. Centre 19 appears to have the only significant effect; and hence, the observed frailty effect is solely due to this one centre having a lower hazard than the other centres included in the study. This, perhaps, explains why the rAIC is only slightly lower than the NF model.

Focusing on the results from the ScF model, both variables included in the model have significant scale coefficients but non-significant shape coefficients. The variable Chemotherapy has a negative coefficient and so the hazard is significantly reduced, time until recurrence is prolonged for patients that received chemotherapy relative to those who did not. Prior Recurrence has a positive coefficient and so having had a recurrence already significantly increases the risk of another recurrence, relative to it being a primary occurrence. The shape coefficients corresponding to the two variables, albeit non-significant, are both positive, suggesting that the effect of chemotherapy wears off over time, while the longer one survives after a recurrence the higher the risk of another recurrence. Plots of the hazard ratios along with their bootstrapped 95% confidence intervals can be seen in Figure 5. As expected, the two variables have approximately constant hazards over time and so perhaps the model can be reduced to a proportional hazards model.

Figure 4: The random effects of the 21 centres with 95% confidence intervals under the ScF model. Centres are sorted in increasing order based on the number of patients.
Figure 5: The hazard ratios with 95% confidence intervals based on the ScF model. The modal values were used in the computation of these hazard ratios (Chemotherapy = yes, Prior Recurrence = no).

5 Discussion and Conclusions

The MPR frailty modelling framework we have proposed not only includes frailty structures not previously explored in the literature, but also generalizes a variety of existing sub-models. Existing literature on MPR frailty models has been limited to multiplicative frailty, and, to the best of our knowledge, a model with correlated frailty in each distributional parameter has not previously been considered. We believe that this is a natural structure in the context of MPR modelling, since it has been shown that estimates of the scale and shape can be quite correlated in practice (Burke and MacKenzie, 2017); hence, it is useful to allow for the possibility that correlation may propagate to the frailty terms.

Although the numerical studies were carried out on a Weibull MPR model, we have developed the model and the estimation procedure in a generic form; the underlying cumulative hazard function can be replaced with that of any other two parameter distribution. In principle, the methods can also be extended to models with more than two distributional parameters, e.g., the power generalized Weibull model of Burke et al. (2020), using a higher order multivariate normal distribution but this is beyond the scope of this paper. The adopted h-likelihood framework provides a computationally inexpensive and straightforward two step procedure to fit our frailty models while avoiding the often intractable integration of the random effects over the frailty distribution. Moreover, the readily available estimates of the frailties allow for the survivor function for individuals with specific characteristics to be estimated, and this is useful in providing information about the merits of the different centres in terms of patient survival in multi-centre studies.

While the proposed MPR framework provides a very flexible approach to modelling correlated survival data at a minimal computational cost, there are various ways in which we can extend it to handle more complicated frailty structures. All of the models that we have considered have a constant frailty variance, perhaps a natural next step for us is to allow the frailty variance to depend on covariates in a similar fashion to that of Peng et al. (2020). Another potential direction worth exploring would be multilevel or nested frailty structures, we can have data on patients, nested within centres or hospitals, with recurrent event times and hence we may need a frailty component for the patients and a frailty component for the centre or hospital. The procedures we present can be straightforwardly extended to include more than one random component for each distributional parameter.


The first author is funded by the Irish Research Council. The second author is supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2020R1F1A1A01056987).


  • J. C. Abrahantes, C. Legrand, T. Burzykowski, P. Janssen, V. Ducrocq, and L. Duchateau (2007) Comparison of different estimation procedures for proportional hazards model with random effects. Computational Statistics & Data Analysis 51 (8), pp. 3913–3930. Cited by: §1.
  • K. Burke, F. Eriksson, and C. Pipper (2019) Semiparametric multiparameter regression survival modeling. Scandinavian Journal of Statistics 47 (2), pp. 555–571. Cited by: §1.
  • K. Burke, M. Jones, and A. Noufaily (2020) A flexible parametric modelling framework for survival analysis. Journal of the Royal Statistical Society: Series C (Applied Statistics) 69 (2), pp. 429–457. Cited by: §1, 3rd item, §5.
  • K. Burke and G. MacKenzie (2017) Multi-parameter regression survival modeling: an alternative to proportional hazards. Biometrics 73, pp. 678–686. Cited by: §1, 3rd item, §2.1, §4.1, §5.
  • D. G. Clayton (1978) A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika 65 (1), pp. 141–151. Cited by: §1.
  • D. R. Cox (1972) Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological) 34 (2), pp. 187–202. Cited by: 1st item.
  • A. C. Davison and D. V. Hinkley (1997) Bootstrap methods and their application. Cambridge University Press. Cited by: §4.2.
  • L. Duchateau and P. Janssen (2008) The frailty model. Springer Science & Business Media. Cited by: §1, §1, §1, 2nd item, §4.2.
  • D. S. Ettinger, D. M. Finkelstein, M. D. Abeloff, J. C. Ruckdeschel, S. C. Aisner, and J. C. Eggleston (1990) A randomized comparison of standard chemotherapy versus alternating chemotherapy and maintenance versus no maintenance therapy for extensive-stage small-cell lung cancer: a phase iii study of the eastern cooperative oncology group.. Journal of Clinical Oncology 8 (2), pp. 230–240. Cited by: §4.2.
  • K. Goethals, P. Janssen, and L. Duchateau (2008) Frailty models and copulas: similarities and differences. Journal of Applied Statistics 35 (9), pp. 1071–1079. Cited by: §1.
  • R. J. Gray (1994) A bayesian analysis of institutional effects in a multicenter cancer clinical trial. Biometrics, pp. 244–253. Cited by: §4.2.
  • I. D. Ha, J. Jeong, and Y. Lee (2017) Statistical modelling of survival data with random effects. Springer. Cited by: §1, 7th item, §2.2, §2.3, §2.3, §4.1, §4.3.
  • I. D. Ha, Y. Lee, and G. MacKenzie (2007) Model selection for multi-component frailty models. Statistics in Medicine 26 (26), pp. 4790–4807. Cited by: §4.1.
  • I. D. Ha, Y. Lee, and J. Song (2001) Hierarchical likelihood approach for frailty models. Biometrika 88 (1), pp. 233–233. Cited by: §1, §2.1, §2.2.
  • I. D. Ha, Y. Lee, and J. Song (2002) Hierarchical-likelihood approach for mixed linear models with censored data. Lifetime Data Analysis 8 (2), pp. 163–176. Cited by: §1, §2.3.
  • I. D. Ha and Y. Lee (2003) Estimating frailty models via poisson hierarchical generalized linear models. Journal of Computational and Graphical Statistics 12 (3), pp. 663–681. Cited by: §1, §2.1.
  • I. D. Ha, M. Noh, and Y. Lee (2012) FrailtyHL: a package for fitting frailty models with hlikelihood. R Journal 4 (2), pp. 28–36. Cited by: §4.3.
  • I. D. Ha, R. Sylvester, C. Legrand, and G. MacKenzie (2011) Frailty modelling for survival data from multi-centre clinical trials. Statistics in Medicine 30 (17), pp. 2144–2159. Cited by: §2.1, §4.3.
  • I. D. Ha, F. Vaida, and Y. Lee (2016) Interval estimation of random effects in proportional hazards models with frailties. Statistical Methods in Medical Research 25 (2), pp. 936–953. Cited by: §1, §2.1, §2.3, §4.2.
  • D. A. Harville (1977) Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association 72 (358), pp. 320–338. Cited by: §2.3.
  • P. Hougaard (2012) Analysis of multivariate survival data. Springer Science & Business Media. Cited by: §1, §1.
  • M. Jones, A. Noufaily, and K. Burke (2020) A bivariate power generalized weibull distribution: a flexible parametric model for survival analysis. Statistical Methods in Medical Research 29 (8), pp. 2295–2306. Cited by: §1, 4th item.
  • Y. Lee, J. A. Nelder, and Y. Pawitan (2017) Generalized linear models with random effects: unified analysis via h-likelihood. Vol. 153, CRC Press. Cited by: §2.2, §2.3, §4.1.
  • Y. Lee and J. A. Nelder (1996) Hierarchical generalized linear models. Journal of the Royal Statistical Society: Series B (Methodological) 58 (4), pp. 619–656. Cited by: §1, §2.2, §2.3.
  • Y. Lee and J. A. Nelder (2001) Hierarchical generalised linear models: a synthesis of generalised linear models, random-effect models and structured dispersions. Biometrika 88 (4), pp. 987–1006. Cited by: §2.3.
  • R. A. Maller and X. Zhou (2003) Testing for individual heterogeneity in parametric models for event history data. Mathematical Methods of Statistics 12 (3), pp. 276–304 (English). External Links: ISSN 1066-5307 Cited by: §3, §4.2.
  • H. D. Patterson and R. Thompson (1971) Recovery of inter-block information when block sizes are unequal. Biometrika 58 (3), pp. 545–554. Cited by: §2.3.
  • D. Peng, G. MacKenzie, and K. Burke (2020) A multiparameter regression model for interval-censored survival data. Statistics in Medicine 39 (14), pp. 1903–1918. Cited by: §1, 4th item, §5.
  • R. A. Rigby and D. M. Stasinopoulos (2005) Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society: Series C (Applied Statistics) 54 (3), pp. 507–554. Cited by: §1.
  • S. Ripatti and J. Palmgren (2000) Estimation of multivariate frailty models using penalized partial likelihood. Biometrics 56 (4), pp. 1016–1022. Cited by: §1.
  • S. G. Self and K. Liang (1987) Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association 82 (398), pp. 605–610. Cited by: §4.2.
  • D. M. Stasinopoulos, R. A. Rigby, et al. (2007) Generalized additive models for location scale and shape (gamlss) in r. Journal of Statistical Software 23 (7), pp. 1–46. Cited by: §1.
  • M. D. Stasinopoulos, R. A. Rigby, and F. D. Bastiani (2018) GAMLSS: a distributional regression approach. Statistical Modelling 18 (3-4), pp. 248–273. Cited by: §1.
  • D. O. Stram and J. W. Lee (1994) Variance components testing in the longitudinal mixed effects model. Biometrics, pp. 1171–1177. Cited by: §4.2.
  • R. J. Sylvester, A. P. van der Meijden, W. Oosterlinck, J. A. Witjes, C. Bouffioux, L. Denis, D. W. Newling, and K. Kurth (2006) Predicting recurrence and progression in individual patients with stage ta t1 bladder cancer using eortc risk tables: a combined analysis of 2596 patients from seven eortc trials. European Urology 49 (3), pp. 466–477. Cited by: §4.3.
  • F. Vaida and S. Blanchard (2005) Conditional akaike information for mixed-effects models. Biometrika 92 (2), pp. 351–370. Cited by: §4.1.
  • F. Vaida and R. Xu (2000) Proportional hazards model with random effects. Statistics in Medicine 19 (24), pp. 3309–3324. Cited by: §1, §4.2, §4.2.
  • R. Xu, F. Vaida, and D. P. Harrington (2009) Using profile likelihood for semiparametric model selection with application to proportional hazards mixed models. Statistica Sinica 19 (2), pp. 819. Cited by: §4.1.

Appendix A Other Simulation Studies

This section displays analogous simulation results to those of Table 2 from the main paper but for different frailty dispersion parameters and/or censoring rate.

Mean Mean Mean Mean Mean Mean Mean Mean Mean
(SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE)
() (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE)
True 1 -0.5 0.5 0.5 0.5 -0.5 2 1 -0.5
(20, 5) 1.37 -0.64 0.64 0.69 0.48 -0.49 2.30 1.00 -0.46
(0.61) (0.28) (0.30) (0.28) (0.09) (0.09) (0.57) (0.20) (0.25)
(0.56) (0.22) (0.22) (0.25) (0.07) (0.07) (0.38) (0.17) (0.18)

(20, 20)
1.05 -0.51 0.51 0.54 0.50 -0.50 2.02 0.99 -0.50
(0.48) (0.09) (0.09) (0.23) (0.03) (0.03) (0.38) (0.18) (0.19)
(0.46) (0.08) (0.08) (0.23) (0.03) (0.03) (0.33) (0.16) (0.16)

(20, 50)
0.99 -0.50 0.51 0.52 0.50 -0.50 1.97 0.98 -0.48
(0.45) (0.05) (0.05) (0.23) (0.02) (0.02) (0.35) (0.16) (0.19)
(0.44) (0.05) (0.05) (0.22) (0.02) (0.02) (0.32) (0.16) (0.17)

(100, 5)
1.27 -0.59 0.58 0.64 0.49 -0.49 2.12 1.03 -0.51
(0.25) (0.11) (0.11) (0.11) (0.04) (0.04) (0.22) (0.09) (0.10)
(0.23) (0.09) (0.09) (0.11) (0.03) (0.03) (0.15) (0.07) (0.08)

(100, 20)
1.04 -0.51 0.52 0.54 0.50 -0.50 2.01 1.00 -0.50
(0.20) (0.04) (0.03) (0.10) (0.01) (0.01) (0.16) (0.07) (0.08)
(0.21) (0.04) (0.04) (0.10) (0.01) (0.01) (0.14) (0.07) (0.08)

(100, 50)
1.03 -0.50 0.50 0.51 0.50 -0.50 2.00 1.00 -0.51
(0.21) (0.03) (0.03) (0.12) (0.01) (0.01) (0.16) (0.08) (0.09)
(0.21) (0.02) (0.02) (0.13) (0.01) (0.01) (0.14) (0.07) (0.07)

Table 5: Averaged coefficient estimates, standard deviations (SE) and the average standard errors (SEE) for the simulation scenario with dispersion parameters , and for censoring rate.
Mean Mean Mean Mean Mean Mean Mean Mean Mean
(SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE)
() (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE)
True 1 -0.5 0.5 0.5 0.5 -0.5 0.5 0.25 -0.5
(20, 5) 1.22 -0.58 0.57 0.62 0.49 -0.50 0.62 0.28 -0.35
(0.26) (0.23) (0.22) (0.13) (0.10) (0.10) (0.25) (0.13) (0.60)
(0.22) (0.19) (0.19) (0.12) (0.09) (0.09) (0.12) (0.06) (0.20)

(20, 20)
1.07 -0.52 0.52 0.53 0.50 -0.50 0.50 0.24 -0.50
(0.13) (0.08) (0.08) (0.08) (0.04) (0.04) (0.13) (0.07) (0.35)
(0.14) (0.08) (0.08) (0.07) (0.04) (0.04) (0.09) (0.04) (0.17)

(20, 50)
1.02 -0.50 0.51 0.51 0.50 -0.50 0.49 0.25 -0.51
(0.12) (0.05) (0.05) (0.06) (0.03) (0.03) (0.09) (0.05) (0.24)
(0.12) (0.05) (0.05) (0.06) (0.03) (0.03) (0.08) (0.04) (0.17)

(100, 5)
1.16 -0.54 0.54 0.56 0.50 -0.49 0.57 0.26 -0.43
(0.10) (0.08) (0.08) (0.06) (0.04) (0.04) (0.10) (0.07) (0.35)
(0.09) (0.08) (0.08) (0.05) (0.04) (0.04) (0.05) (0.02) (0.11)

(100, 20)
1.04 -0.51 0.51 0.52 0.50 -0.50 0.50 0.25 -0.50
(0.06) (0.04) (0.04) (0.03) (0.02) (0.02) (0.05) (0.03) (0.15)
(0.06) (0.03) (0.03) (0.03) (0.02) (0.02) (0.04) (0.02) (0.08)

(100, 50)
1.02 -0.50 0.50 0.51 0.50 -0.50 0.50 0.25 -0.50
(0.05) (0.02) (0.02) (0.03) (0.01) (0.01) (0.04) (0.02) (0.11)
(0.05) (0.02) (0.02) (0.03) (0.01) (0.01) (0.04) (0.02) (0.08)
Table 6: Averaged coefficient estimates, standard deviations (SE) and the average standard errors (SEE) for the simulation scenario with dispersion parameters , and for censoring rate.
Mean Mean Mean Mean Mean Mean Mean Mean Mean
(SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE)
() (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE)
True 1 -0.5 0.5 0.5 0.5 -0.5 1 0.5 0.5

(20, 5)
1.34 -0.62 0.61 0.66 0.48 -0.48 1.21 0.54 0.47
(0.39) (0.29) (0.28) (0.17) (0.12) (0.12) (0.38) (0.17) (0.37)
(0.33) (0.21) (0.21) (0.16) (0.09) (0.09) (0.21) (0.10) (0.16)

(20, 20)
1.08 -0.51 0.51 0.54 0.50 -0.50 1.00 0.50 0.49
(0.25) (0.08) (0.08) (0.13) (0.04) (0.04) (0.19) (0.10) (0.21)
(0.24) (0.08) (0.08) (0.12) (0.04) (0.04) (0.16) (0.08) (0.17)
(20, 50) 1.04 -0.51 0.50 0.51 0.50 -0.50 0.98 0.49 0.47
(0.23) (0.05) (0.05) (0.12) (0.02) (0.03) (0.17) (0.09) (0.19)
(0.22) (0.05) (0.05) (0.11) (0.02) (0.02) (0.16) (0.08) (0.17)
(100, 5) 1.22 -0.56 0.56 0.61 0.48 -0.49 1.09 0.52 0.49
(0.15) (0.10) (0.10) (0.07) (0.04) (0.04) (0.13) (0.063) (0.13)
(0.13) (0.09) (0.09) (0.07) (0.04) (0.04) (0.08) (0.04) (0.08)
(100, 20) 1.06 -0.51 0.51 0.53 0.50 -0.50 1.00 0.50 0.48
(0.11) (0.04) (0.04) (0.06) (0.02) (0.02) (0.08) (0.04) (0.09)
(0.11) (0.04) (0.04) (0.05) (0.02) (0.02) (0.07) (0.04) (0.08)
(100, 50) 1.03 -0.50 0.50 0.51 0.50 -0.50 1.00 0.50 0.49
(0.10) (0.02) (0.02) (0.05) (0.01) (0.01) (0.07) (0.04) (0.08)
(0.10) (0.02) (0.02) (0.05) (0.01) (0.01) (0.07) (0.04) (0.08)
Table 7: Averaged coefficient estimates, standard deviations (SE) and the average standard errors (SEE) for the simulation scenario with dispersion parameters , and for censoring rate.
Mean Mean Mean Mean Mean Mean Mean Mean Mean
(SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE)
() (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE)
True 1 -0.5 0.5 0.5 0.5 -0.5 1 0.5 -0.5
(20) 1.1 -0.51 0.51 0.54 0.50 -0.50 1.03 0.50 -0.47
(0.26) (0.07) (0.07) (0.12) (0.04) (0.04) (0.20) (0.10) (0.24)
(0.25) (0.07) (0.07) (0.12) (0.03) (0.03) (0.17) (0.08) (0.17)
(100) 1.08 -0.51 0.51 0.52 0.50 -0.50 1.00 0.50 -0.50
(0.11) (0.03) (0.03) (0.06) (0.02) (0.01) (0.09) (0.05) (0.10)
(0.11) (0.03) (0.03) (0.06) (0.01) (0.01) (0.07) (0.04) (0.08)
Table 8: Averaged coefficient estimates, standard deviations (SE) and the average standard errors (SEE) for the simulation scenario with dispersion parameters , , , censoring rate and varying cluster sizes ( of the clusters are of size , are of size and the remaining are of size ).
Mean Mean Mean Mean Mean Mean Mean Mean Mean
(SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE) (SE)
() (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE) (SEE)
True 1 -0.5 0.5 0.5 0.5 -0.5 1 0.5 -0.5
(20, 5) 1.59 -0.79 0.79 0.68 0.45 -0.44 1.55 0.61 -0.04
(0.72) (0.66) (0.65) (0.27) (0.20) (0.21) (0.66) (0.27) (0.60)
(0.49) (0.39) (0.38) (0.20) (0.13) (0.13) (0.29) (0.12) (0.19)

(20, 20)
1.11 -0.52 0.51 0.53 0.50 -0.50 0.99 0.50 -0.48
(0.29) (0.16) (0.15) (0.14) (0.06) (0.07) (0.22) (0.12) (0.29)
(0.26) (0.14) (0.14) (0.13) (0.06) (0.06) (0.17) (0.08) (0.18)

(20, 50)
1.02 -0.51 0.51 0.52 0.50 -0.50 0.98 0.50 -0.49
(0.22) (0.08) (0.08) (0.12) (0.04) (0.04) (0.19) (0.09) (0.23)
(0.24) (0.08) (0.08) (0.12) (0.04) (0.04) (0.16) (0.08) (0.17)

(100, 5)
1.38 -0.62 0.63 0.63 0.47 -0.47 1.16 0.56 -0.24
(0.22) (0.23) (0.22) (0.11) (0.08) (0.08)