MEGH: A parametric class of general hazard models for clustered survival data

Abstract

In many applications of survival data analysis, the individuals are treated in different medical centres or belong to different clusters defined by geographical or administrative regions. The analysis of such data requires accounting for between-cluster variability. Ignoring such variability would impose unrealistic assumptions in the analysis and could affect the inference on the statistical models. We develop a novel parametric mixed-effects general hazard (MEGH) model that is particularly suitable for the analysis of clustered survival data. The proposed structure generalises the mixed-effects proportional hazards and mixed-effects accelerated failure time structures, among other structures, which are obtained as special cases of the MEGH structure. We develop a likelihood-based algorithm for parameter estimation in general subclasses of the MEGH model, which is implemented in our R package MEGH. We propose diagnostic tools for assessing the random effects and their distributional assumption in the proposed MEGH model. We investigate the performance of the MEGH model using theoretical and simulation studies, as well as a real data application on leukaemia.

1 Introduction

In the analysis of time-to-event data, it is common to obtain clustered samples. For instance, in medical statistics, one may be interested in modelling the survival times of a sample of patients who receive attention in different health centres with varying characteristics or when the patients live in different geographic or administrative regions. This clustering or grouping information needs to be incorporated in any survival model to avoid under-estimation of the standard errors of the parameter estimators. Moreover, ignoring unobserved heterogeneity, due to clustering, in hazard regression models has an effect on model misspecification as the marginal model obtained by integrating out random effects may not coincide with the model ignoring clustering.^1–3 A simple way to account for clustered survival data consists of fitting a stratified survival model. This is typically done by fitting a proportional hazards (PHs) model for each strata level with different baseline hazards but sharing the same regression coefficients. ⁴ Alternatively, when the number of clusters is small, the cluster indicator can be coded as a categorical variable. The main disadvantage of these approaches is that they do not account for between-cluster variability. Modelling the clustering effect is more formally done by incorporating random effects, which are random variables that capture the hierarchical structure of the data and account for unobserved variability or heterogeneity between clusters. Regression models, including survival models, that incorporate random effects are often referred to as mixed-effects models (we refer the reader to McCulloch et al. ⁵ and the references therein for an extensive review of this kind of models). The most common mixed-effects models for time-to-event data correspond to extensions of common regression survival models. For example, the mixed-effects PHs (MEPHs) model^6,7 is an extension of the semiparametric Cox PHs models. This model is obtained by adding random effects, indexed by the cluster indicator, to the linear predictor function in the Cox model. Another common strategy to capture clustering consists of adding random effects to the linear predictor function in parametric and semiparametric accelerated failure time (AFT) models, leading to the so-called mixed-effects AFT (MEAFT) model (see Rubio and Steel ⁸ and Do Ha and Lee ⁹ for a review on these models). A related class of mixed-effects models is constructed by adding random effects or Gaussian processes to the linear predictor function in hazard-based (parametric and semiparametric) regression models.^10,11 Some examples of these kinds of models are the mixed-effects spline regression models formulated at the cumulative hazard level by Crowther et al.,^12,13 Wood, ¹⁴ or those formulated at the hazard level by Vaida and Xu ¹⁵ and Charvat et al. ¹⁶ In a related vein, Zhou and Hanson ¹⁷ proposed mixed-effects models by adding random effects to PHs, proportional odds and AFT models in a Bayesian semiparametric framework. These classes of models are thus characterised by aiming at capturing the clustering effect via the incorporation of random effects into the linear predictor of different sorts of fixed effects hazard-based regression models.

A natural question about mixed-effects models is the effect of misspecification of the random effects structure on the estimation of quantities of interest. It has been found that misspecifying the random-effects distribution has a negligible effect on the point estimation of the fixed effects parameters of some classes of random effects models. ¹⁸ However, there are more tangible effects of this misspecification on interval estimation of the parameters and predictions^8,19 as well as point and interval estimation of the parameters of the baseline hazard. ²⁰ More crucially, ignoring a significant random effect could affect both parameter estimation and inference. ²¹ Tools for diagnosing misspecification of the random-effects part have been developed for mixed-effects models (without censoring) in recent years (see the authors in^22–24).

We develop a novel mixed-effects general hazard (MEGH) structure which incorporates random effects into a hazard-based regression model^25,26 at the hazard scale and the time scale. The resulting model can capture clustering affecting the hazard scale and the time scale, and contains the MEPH and MEAFT models as particular cases. We focus our attention on two tractable subclasses of mixed-effects models that generalise the MEPH and MEAFT models. The specification of the MEGH model requires modelling a baseline hazard, for which we employ several flexible parametric distributions. We provide diagnostic tools for the random-effects part of the proposed classes of models in the presence of censoring. We conduct extensive simulation studies that reveal that not only misspecifying the random-effects distribution is problematic, but also misspecifying the role of the random effects in the hazard structure has adverse effects.

The remainder of the paper is organised as follows. In Section 2, we describe the proposed MEGH model and discuss two submodels of interest. In Section 3, we present the expressions for the conditional and marginal likelihood functions, and define the marginal maximum likelihood estimator (MMLE), which defines the estimation approach followed in this paper. We also present theoretical results regarding the consistency and asymptotic normality of the MMLE under standard regularity conditions. Section 4 presents a test for the inclusion of random effects, as well as a graphical diagnostic tool to assess the goodness-of-fit of the random-effects distribution. Section 5 presents an extensive simulation study which illustrates the performance of the proposed methodology as well as an exploration of the effects of misspecifying the distribution and the role of the random effects. Section 6 presents an application with real data in the context of cancer survival with a spatial component. We conclude with a general discussion and possible extensions of our work in Section 7.

2 MEGH: Mixed-effects general hazard models

In this section, we describe the proposed MEGH model, and discuss two general subclasess of interest. To this end, let us first introduce some notation. Let $o_{ij}\in {\mathbb{R}}_{+}$ , $i=1,\dots ,r$ indicates the cluster, $j=1,\dots ,n_i$ denotes the individuals, be a sample of times-to-event of interest, and ${\mathbf{x}}_{ij}\in {\mathbb{R}}^p$ be a vector of covariates associated with the $ij$ th individual. Let $c_{ij}\in {\mathbb{R}}_{+}$ be the right-censoring times, and $t_{ij}=min\{o_{ij},c_{ij}\}$ be the observed survival times. Let $d_{ij}=I(o_{ij}<c_{ij})$ be the vital status indicators (1 - dead, 0 - alive), and $n=\sum _{i=1}^r{n}_i$ be the total sample size across $r$ clusters. We propose the MEGHs model, which is defined through the individual hazard function:

$$h(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i)=h_0\left(t_{ij}\exp \{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+{\tilde{u}}_i\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }+u_i\}$$

(1)

where $h_0(\cdot )$ is a baseline hazard, $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_p{)}^{\mathrm{\top }}\in {\mathbb{R}}^p$ , ${{\mathit{\alpha }=({\alpha }_1,\dots ,\alpha }_{\tilde{p}}{)}^{\mathrm{\top }}\in \mathbb{R}}^{\tilde{p}}$ , ${\mathbf{x}}_{ij}\in {\mathbb{R}}^p$ is a vector of covariates affecting the hazard scale, ${\tilde{\mathbf{x}}}_{ij}\in {\mathbb{R}}^{\tilde{p}}$ is a vector of covariates affecting the time scale (typically, ${\tilde{\mathbf{x}}}_{ij}\subset {\mathbf{x}}_{ij}$ ), and $(u_i,{\tilde{u}}_i)\stackrel{iid}{\sim }G$ , where $G$ is a continuous distribution with support on ${\mathbb{R}}^2$ and zero mean. Denote by $\mathbf{X}=\{{\mathbf{x}}_{ij}\}$ , $\tilde{\mathbf{X}}=\{{\tilde{\mathbf{x}}}_{ij}\}$ , the design matrices, and by ${\mathbf{X}}_o,{\tilde{\mathbf{X}}}_o$ the sub-matrices with the rows associated with the uncensored individuals. To avoid collinearity problems, we assume that the matrices associated with uncensored individuals are full column rank. The MEGH model represents an extension of the general hazard (GH) structure by Chen and Jewell, ²⁵ which is obtained by removing the random effects from (1). Clearly, the GH structure contains the PHs model for $\mathit{\alpha }=0$ , the AFT model for $\mathit{\alpha }=\mathit{\beta }$ , and the accelerated hazards (AHs) model for $\mathit{\beta }=0$ . ²⁵ Consequently, we also need to impose that the baseline hazard $h_0$ does not belong to the Weibull family, as in this case the GH structure, and consequently the MEGH structure, becomes non-identifiable. ²⁵ Random effects are typically incorporated at the hazard scale (MEPH) or simultaneously at the hazard and time scales (MEAFT). The MEGH structure allows for incorporating random effects in both scales separately. The MEGH structure can also be seen as a generalisation of shared frailty survival models, which aim at accounting for unobserved heterogeneity between clusters.^27,16 Another appealing feature of the MEGH formulation is that the cumulative hazard function can be written in closed-form as follows

$$H(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i)=H_0\left(t_{ij}\exp \{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+{\tilde{u}}_i\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }-{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+u_i-{\tilde{u}}_i\}$$

The MEGH structure (1) is a very rich mixed hazard structure that contains, as particular cases, a number of hazard models of interest. In this paper, we particularly focus on the following two subclasses.

2.1 MEGH-I: Mixed Structure I

The first subclass, MEGH-I, is defined by the hazard and cumulative hazard functions

$$\begin{array}{rl}{h}_1(t_{ij}\mid {\mathbf{x}}_{ij},u_i)& =h_0\left(t_{ij}\exp \left\{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }\right\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }+u_i\}\\ H_1(t_{ij}\mid {\mathbf{x}}_{ij},u_i)& =H_0\left(t_{ij}\exp \left\{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }\right\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }-{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+u_i\}\end{array}$$

(2)

The mixed hazard structure (2) contains the MEPH model ( $\mathit{\alpha }=0$ ), ⁶ while allowing for the inclusion of time-dependent effects through the parameter $\mathit{\alpha }$ .

2.2 MEGH-II: Mixed Structure II

The second subclass, MEGH-II, is defined by the hazard and cumulative hazard functions

$$\begin{array}{cc}{h}_2(t_{ij}\mid {\mathbf{x}}_{ij},u_i)& =h_0\left(t_{ij}\exp \{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+u_i\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }+u_i\}\\ H_2(t_{ij}\mid {\mathbf{x}}_{ij},u_i)& =H_0\left(t_{ij}\exp \{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+u_i\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }-{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }\}\end{array}$$

(3)

The MEGH-II structure (3) generalises the MEAFT structure ( $\mathit{\alpha }=\mathit{\beta }$ ).

Figure 1 shows the link between the MEGH and some submodels of interest, while the proposed MEGH structure is much richer than illustrated in this figure.

Figure 1.

Diagram of the link between the MEGH and some particular hazard structures of interest. (AH: accelerated hazards; AFT: accelerated failure time; PH: proportional hazards; GH: general hazard; MEPH: mixed-effects proportional hazards; MEAFT: mixed-effects accelerated failure time; MEGH: mixed-effects general hazard).

We will adopt parametric baseline hazards $h_0(\cdot \mid \mathit{\theta })$ for models (2)–(3), based on flexible parametric distributions. This is an important step as the choice of the parametric baseline hazard determines the hazard shapes one can obtain. For instance, the log-normal hazard function is unimodal (up-then-down), while the Gamma hazard function can be increasing, decreasing or flat. There exists other (three-parameter) distributions that can capture the basic shapes of the hazard (increasing, decreasing, unimodal, and bathtub), such as the Exponentiated Weibull, Generalised Gamma and Power Generalised Weibull (PGW). However, it is important to keep in mind that an efficient estimation of the parameters of more flexible distributions typically requires larger sample sizes.^26,28 On the other hand, parametric models have been found to perform well, compared to semiparametric counterparts, in the context of clustered survival data. ²⁹

3 Likelihood function and parameter estimation

We here describe the conditional and marginal likelihood functions, which will be used to calculate the marginal maximum likelihood estimators as well as to develop the diagnostic tools presented in the next section. We discuss the details behind the calculation of the marginal likelihoods and the optimisation methods. Finally, we present a theoretical result on the consistency and asymptotic normality of the marginal maximum likelihood estimators.

3.1 Likelihood and marginal likelihood functions

Let ${\mathbf{t}}_i=(t_{i1},\dots ,t_{i{n}_i})$ be the sample associated with the $i$ th cluster and ${\mathbf{x}}_i=({\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{i{n}_i})$ denote the corresponding matrix of covariates. The cluster-specific marginal likelihood function associated with the $i$ th cluster, after integrating out the random effects $u_i$ and ${\tilde{u}}_i$ , is given by

$${\int }_{{\mathbb{R}}^2}\prod _{j=1}^{n_i}h(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i{)}^{d_{ij}}\exp \{-H(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i)\}dG(u_i,{\tilde{u}}_i)$$

For the case where the random-effects distribution $G$ is assumed to be a multivariate distribution with mean $\mathbf{0}$ and finite variance, with parameters $\mathit{\xi }$ , we can write down the marginal likelihood function of the MEGH model (1) in terms of the hazard and cumulative hazard functions as follows. To simplify notation, let us denote $\mathit{\eta }=(\mathit{\beta },\mathit{\alpha },\mathit{\theta },\mathit{\xi })$ and let $\mathrm{\Gamma }$ be the parameter space. The marginal likelihood function is

$$m(\mathit{\eta })=\prod _{i=1}^r{m}_i(\mathit{\eta })$$

(4)

where $m_i$ represents the cluster-specific marginal likelihood associated with the $i$ th cluster

$$m_i(\mathit{\eta })={\int }_{{\mathbb{R}}^2}\exp \{{\ell }_i(\mathit{\eta },u_i,{\tilde{u}}_i)\}dG(u_i,{\tilde{u}}_i;\mathit{\xi })$$

(5)

and

$${\ell }_i(\mathit{\eta },u_i,{\tilde{u}}_i)=\sum _{j=1}^{n_i}{d}_{ij}\log h(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i)-\sum _{j=1}^{n_i}H(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i)$$

denotes the log-likelihood function conditional on the random effects $u_i$ and ${\tilde{u}}_i$ . The marginal maximum likelihood estimator (MMLE) is then defined as $\hat{\mathit{\eta }}={\mathrm{argmax}}_{\mathrm{\Gamma }}m(\mathit{\eta })$ .

3.2 Computations

To calculate the MMLE $\hat{\mathit{\eta }}$ , we need to evaluate the marginal likelihood function (4), which is a product of integrals that define the cluster-specific marginal likelihoods (5). The integrand in (5) can be small when the number of individuals belonging to that cluster is moderate or large. To deal with this numerical challenge, we scale the integrand by using

$$m_i(\mathit{\eta })=K_i{\int }_{{\mathbb{R}}^2}\exp \{{\ell }_i(\mathit{\eta },u_i,{\tilde{u}}_i)-\log K_i\}dG(u_i,{\tilde{u}}_i)$$

where $\log K_i=\underset{u_i,{\tilde{u}}_i}{max}{\ell }_i(\mathit{\eta },u_i,{\tilde{u}}_i)$ , for fixed values of $\mathit{\eta }$ . This allows for a more stable numerical integration, as the integrand is now bounded between $(0,1]$ (always taking the value $1$ at the maximum). Maximisation over the random effects is simpler as their variance is rarely large. Numerical integration in our empirical examples, which include only one hierarchical level, is done by using the R command integrate with default options. Nonetheless, other methods, such as Laplace approximations or Monte Carlo integration, could also be used. The optimisation step is carried out using standard R routines, such as those implemented in the commands optim and nlminb. The R package MEGH (https://github.com/FJRubio67/MEGH) implements the MEGH-I and MEGH-II models under these specifications.

3.3 Theory

Here, we prove the consistency and asymptotic normality of the MMLE, when $r\to \mathrm{\infty }$ (which implicitly implies $n\to \mathrm{\infty }$ ), under Assumptions A1–A7 given in the online Supplemental Material. We denote by ${\mathit{\eta }}^{\star }=({\mathit{\beta }}^{\star },{\mathit{\alpha }}^{\star },{\mathit{\theta }}^{\star },{\mathit{\xi }}^{\star })$ the unknown true value of the parameters, which is assumed to be an interior point of the parameter space $\mathrm{\Gamma }$ .

Theorem 1

Consider the MEGH model (1). Let $h_0(\cdot \mid \mathit{\theta })$ be a pre-specified parametric baseline hazard as in Section 2. Under Assumptions A1–A7 given in the online Supplemental Material, it follows that

Consistency: $\hat{\mathit{\eta }}\stackrel{P}{\to }{\mathit{\eta }}^{\star }$ as $r\to \mathrm{\infty }$ ,

Asymptotic normality: $\sqrt{r}(\hat{\mathit{\eta }}-{\mathit{\eta }}^{\star })\stackrel{d}{\to }N(\mathbf{0},\mathbf{I}({\mathit{\eta }}^{\star }{)}^{-1})$ as $r\to \mathrm{\infty }$ ,

where $\mathbf{I}({\mathit{\eta }}^{\star })$ is the expectation matrix $\mathbf{I}(\mathit{\eta })=\mathrm{cov}[{\mathrm{\nabla }}_{\mathit{\eta }}\log m_1(\mathit{\eta })]$ evaluated at ${\mathit{\eta }}^{\star }$ .

Regarding the required Assumptions A1–A7, we note that Assumption A1 is a restriction to a compact parameter space, which can be assumed to be a large compact set. Assumption A2 concerns the censoring mechanism, which is assumed to be non-informative. Assumption A3 guarantees identifiability and continuity of the MEGH model. Assumption A4 togheter with Assumptions A5–A7 represent standard integrability and differentiability conditions on the marginal likelihood function, which implicitly impose conditions on the choice of the baseline hazard. Assumption A6 assumes the existence and non-singularity of the information matrix. Verifying such condition in practice for specific choices of the parametric baseline hazard is complicated. However, it has been shown that certain smoothness and identifiability conditions, which are easier to verify in practice, guarantee the non-singularity of this matrix (see Asgharian ³⁰ ). Similar asymptotic results have also been obtained for parametric copula survival models for clustered data in Prenen et al., ²⁹ where parametric models were found to perform better than their semiparametric counterparts.

4 Diagnostic tools for random effects in the MEGH model

As already discussed, the proposed MEGH model (1) incorporates random effects to capture the unknown variability between clusters or groups in the survival data. Since random effects are latent and unobservable variables, their assumed distribution can be subject to misspecification, so it is important to check their distributional assumption to identify potential model misspecification. It is also important to test which random effects are statistically significant and should be included in the model. These two problems are crucial especially for practical use, so we study them for the MEGH model in the following two subsections, respectively.

4.1 Detecting misspecification of the random-effects distribution

We here present a diagnostic tool for detecting misspecification of the random-effects distribution in the MEGH model. The assumed random-effects distribution $G$ is typically a multivariate normal distribution. To check the appropriateness of an assumed random-effects distribution, ³¹ suggested to use the so-called gradient function for models with random effects. Their diagnostic tool based on a gradient function is not directly applicable to the MEGH model (1). We extend the gradient function approach for the MEGH model. For this purpose, we use the key idea of the gradient function, that is, to check if the marginal log-likelihood of the model can be increased considerably by replacing $G$ with another random-effects distribution, say $W$ . If so, then the assumed distribution $G$ may not be adequate for the random effects because another random-effects distribution produces a considerably larger marginal likelihood.

Recall that $\mathit{\eta }=(\mathit{\beta },\mathit{\alpha },\mathit{\theta },\mathit{\xi })$ is the combined vector of the model parameters. We here write the marginal likelihood as $m(G)$ , instead of $m(\mathit{\eta })$ , to emphasise its dependence on the assumed random-effects distribution $G$ . To check if the marginal log-likelihood can be increased considerably by replacing $G$ with another random-effects distribution $W$ , we use the directional derivative of the marginal log-likelihood $\log m(.)$ at $G$ into the direction $W$ as follows

$$\begin{array}{rl}\mathrm{\Phi }(G,W)& =\underset{\epsilon \stackrel{>}{\to }0}{lim}\frac{\log m((1-\epsilon )G+\epsilon W)-\log m(G)}{\epsilon }\\ & =\frac{\mathrm{\partial }\log m((1-\epsilon )G+\epsilon W)}{\mathrm{\partial }\epsilon }{|}_{\epsilon =0}\end{array}$$

Then, there is no better random-effects distribution than $G$ if $\mathrm{\Phi }(G,W)\le 0$ for all $W$ . We can write that

$$\begin{array}{rl}\frac{1}{r}\mathrm{\Phi }(G,W)& =\frac{1}{r}\frac{\mathrm{\partial }\sum _{i=1}^r\log [(1-\epsilon )m_i(G)+\epsilon m_i(W)]}{\mathrm{\partial }\epsilon }{|}_{\epsilon =0}\\ & =\frac{1}{r}\sum _{i=1}^r\frac{m_i(W)-m_i(G)}{m_i(G)}\\ & =\frac{1}{r}\sum _{i=1}^r\frac{{\int }_{{\mathbb{R}}^2}\exp \{{\ell }_i(u,\tilde{u})\}dW(u,\tilde{u})}{m_i(G)}-1\\ & ={\int }_{{\mathbb{R}}^2}[\frac{1}{r}\sum _{i=1}^r\frac{\exp \{{\ell }_i(u,\tilde{u})\}}{m_i(G)}]dW(u,\tilde{u})-1\\ & ={\int }_{{\mathbb{R}}^2}\mathrm{\Delta }(G,(u,\tilde{u}))dW(u,\tilde{u})-1\end{array}$$

in which $\mathrm{\Delta }(G,(u,\tilde{u}))$ is the gradient function for the MEGH model, defined as

$$\mathrm{\Delta }(G,(u,\tilde{u}))=\frac{1}{r}\sum _{i=1}^r\frac{\exp \{{\ell }_i(u,\tilde{u})\}}{m_i(G)},(u,\tilde{u})\in {\mathbb{R}}^2$$

(6)

Hence, we have $\varphi (G,W)\le 0$ for all $W$ , if $\mathrm{\Delta }(G,(u,\tilde{u}))\le 1\mathrm{\forall }(u,\tilde{u})\in {\mathbb{R}}^2$ . The gradient function (6) has useful properties. The calculation of the gradient function does not require specifying an alternative distribution $W$ . Also, since $\mathrm{\Phi }(G,G)=0$ , we must have $\mathrm{\Delta }(G,(u,\tilde{u}))=1$ for all random effects points $(u,\tilde{u})$ in the support of $G$ . As a graphical diagnostic tool, we plot the gradient function versus random effects points $(u,\tilde{u})$ , and if the gradient plot does not exceed $1$ then the assumed random-effects distribution $G$ will be deemed to be adequate for random effects; otherwise, the random-effects distribution is deemed to be misspecified. See ³¹ and ³² who applied a similar approach to generalised linear and nonlinear mixed models. The calculation of the gradient function (6) is straightforward as it only requires the calculation of the marginal and conditional distributions for all $r$ clusters. Note that the dependence of the gradient function on the parameters $\mathit{\eta }=(\mathit{\beta },\mathit{\alpha },\mathit{\theta },\mathit{\xi })$ is suppressed for simplicity of the notation, however, we replace these parameters with their estimates when calculating the gradient function.

4.2 Tests for the need of random effects in the MEGH model

Since random effects are latent and unobservable variables, it is not obvious in practice which random effects must be included in the model. As an initial check, one can use a plot of the survival curves for all clusters to get some understanding of the variability between clusters and the potential need of random effects (see an example in Figure 4). We provide formal tests for verifying which random effects are statistically significant and must be present in the MEGH model. If the test confirms that all the random effects are not significant, then the extended model (1) will not be essential and one could use the reduced models without random effects as in Chen and Jewell ²⁵ and Rubio et al. ²⁶ The test for inclusion or exclusion of random effects can be carried out via a test for zero random effects. From a statistical perspective, to test for zero random effects is equivalent to testing if their variance parameters are equal to $0$ . We want to test whether or not the random effects $u_i$ and ${\tilde{u}}_i$ can be excluded from the MEGH model (1). This hypothesis testing problem can be expressed as follows

$$\{\begin{array}{c}{H}_0:{\sigma }_u^2=0\text{and}{\sigma }_{\tilde{u}}^2=0\\ H_1:{\sigma }_u^2>0\text{or}{\sigma }_{\tilde{u}}^2>0\end{array}$$

(7)

Figure 4.

Leukaemia data: Kaplan-Meier estimators of survival by district.

This is a non-standard testing problem because the null hypothesis places the true variance parameters on the boundary of the parameter space.^33–36 Classical tests such as the likelihood ratio, Wald and score tests suffer from testing on the boundary of the parameter space, because the regularity conditions do not hold under such situations. As a consequence, the usual asymptotic chi-squared distribution of the likelihood ratio or score statistic is incorrect here. The correct asymptotic distribution is well understood to be a mixture of chi-squared distributions when $r\to \mathrm{\infty }$ under some mild conditions.^33,34

The hypotheses $H_0$ and $H_1$ in (7) correspond to Case 7 of Self and Liang. ³³ Thus, under $H_0$ , the correct asymptotic distribution of the likelihood ratio statistic for testing (7) would be $\frac{1}{4}{\chi }_0^2+\frac{1}{2}{\chi }_1+\frac{1}{4}{\chi }_2^2$ , where ${\chi }_0^2$ is a point mass at $0$ , and ${\chi }_1^2$ and ${\chi }_2^2$ are the chi-squared distributions with one and two degrees of freedom, respectively. The $p$ -value of the likelihood ratio test for the hypothesis test (7) is as follows (see also Verbeke and Molenberghs ³⁷ )

$$p=\frac{1}{4}P({\chi }_0\ge R_{\text{obs}})+\frac{1}{2}P({\chi }_1\ge R_{\text{obs}})+\frac{1}{4}P({\chi }_2\ge R_{\text{obs}})$$

(8)

where $R_{\text{obs}}$ denotes the observed likelihood ratio statistic.

For the MEGH model with the mixed hazard structures (2)–(3), since $u_i$ and ${\tilde{u}}_i$ are the same, the hypothesis test (7) collapses to

$$\{\begin{array}{c}{H}_0:{\sigma }_u^2=0\\ H_1:{\sigma }_u^2>0\end{array}$$

The correct asymptotic distribution of the likelihood ratio statistic for this hypothesis test is $\frac{1}{2}{\chi }_0^2+\frac{1}{2}{\chi }_1$ , as it corresponds to Case 5 of Self and Liang. ³³

We will apply the above diagnostic tools to our case study in Section 6. We note that there are also non-asymptotic tests available for testing random effects, which avoid the boundary issue. Bootstrap and permutation tests as well as Bayesian tests do not rely on asymptotic results and hence avoid boundary issues when testing random effects. Those tests have been well studied for linear and generalised linear mixed models^24,36,38; however, such tests may not be easily applied to the case of MEGH model, due to the presence of censored and clustered data which makes bootstrap or permutation challenging. Furthermore, Verbeke and Molenberghs ³⁷ and Claeskens et al. ³⁹ provide asymptotic score tests that can be applied in this context.

5 Simulation studies

In this section, we conduct simulations to evaluate the performance of the proposed MEGH model. In the simulations, we first examine the accuracy of parameter estimation as well as statistical inference with the MEGH model under different scenarios. We investigate how the MEGH model compares with the model ignoring random effects. Such comparisons are crucial in demonstrating the advantages of the proposed MEGH model over the models ignoring random effects. We then assess the impact of misspecifying the mixed structure as well as misspecifying the random-effects distribution. We also examine the finite-sample performance of the diagnostic tests for testing random effects in the MEGH model. Some of the simulation results are reported in the online Supplemental Material for the sake of space.

In the simulations, we study both mixed hazard structures (2)–(3) and consider the following MEGH model, according to our real data application in Section 6,

$$h(t_{ij}\mid {\mathbf{x}}_{ij},u_i,{\tilde{u}}_i)=h_0\left(t_{ij}\exp \{{\alpha }_1{\text{age}}_{ij}+{\tilde{u}}_i\}\right)\exp \{{\beta }_1{\text{age}}_{ij}+{\beta }_2{\text{sex}}_{ij}+{\beta }_3{\text{wbc}}_{ij}+{\beta }_4{\text{tpi}}_{ij}+u_i\}$$

(9)

where the covariates ${\mathbf{x}}_{ij}^{\mathrm{\top }}=({\text{age}}_{ij},{\text{sex}}_{ij},{\text{wbc}}_{ij},{\text{tpi}}_{ij})$ are in accordance with the leukaemia data in Section 6, and $h_0(\cdot )$ is the baseline hazard for which we here consider the PGW baseline hazard with parameters $\mathit{\theta }=(\eta ,\nu ,\delta )$ as in the online Supplemental Material. We also conduct simulations with the log-logistic baseline hazard, which are presented in the Supplemental Material.

Recall that ${\tilde{u}}_i=0$ for the mixed hazard structure MEGH-I, and ${\tilde{u}}_i=u_i$ for the mixed hazard structure MEGH-II. In the simulations, we choose the true parameter values according to the estimates obtained from the leukaemia data application. The number of clusters is 24. We simulate 250 data sets from the above model with each of the hazard structures MEGH-I and MEGH-II, each with the same size as of the real data application (i.e. $n=1043$ ), using our R package MEGH available in the online Supplemental Material. The random effects are generated from a normal distribution with mean $0$ and two standard deviation values of ${\sigma }_u=0.5$ and ${\sigma }_u=1$ . The censoring rate is set to $25\text{\%}$ .

We fit the MEGH model with both mixed hazard structures (2)–(3) separately to the 250 simulated data sets. We also fit the model ignoring random effects to the simulated data sets. We calculate the bias of the parameter estimates from each of these three models. Figures 2 and 3 show the estimation bias for these three models. It can be seen that both MEGH models with mixed hazard structures (2)–(3) produce estimates with considerably smaller bias compared to the model ignoring random effects. The bias of the model ignoring random effects for most parameters is higher due to neglecting the presence of random effects. However, we note that the estimates obtained from models with or without random effects are not generally comparable like with like, as discussed by Lee and Nelder. ⁴⁰ This is because the model with random effects is conditional on random effects, while the model without random effects is marginal (or population-average).

Figure 2.

The bias of the estimates from the three methods: MEGH-I, MEGH-II and MLE, for all the parameters based on $250$ simulation replications when the simulated data are generated from model (9) with the mixed structure I and PGW baseline hazard, and a normal distribution for the generated random effects with ${\sigma }_u=1$ .

Figure 3.

The bias of the estimates from the three methods: MEGH-I, MEGH-II and MLE, for all the parameters based on $250$ simulation replications when the simulated data are generated from model (9) with the mixed structure II and PGW baseline hazard, and a normal distribution for the generated random effects with ${\sigma }_u=1$ .

We also observe from the simulation results presented in Figures 2 and 3 that when the mixed hazard structure is misspecified it leads to relatively higher bias compared to the model with the correct mixed hazard structure. This suggests that the choice of the mixed hazard structure is crucial for achieving a reliable model fit and results. We note that the bias of the MEGH model with the correct mixed hazard structure is very small in the simulations.

Table 1 presents the $95\text{\%}$ confidence intervals for all the regression parameters obtained from the three models considered for the simulation case when the true mixed structure is MEGH-I. The results show that the model with the correct mixed structure MEGH-I produces reliable confidence intervals, but both the model ignoring random effects and the model with the incorrect mixed structure MEGH-II tend to produce inaccurate confidence intervals. In particular, the confidence intervals from the model ignoring random effects mostly do not cover the true parameter values.

Table 1.

The $95\text{\%}$ confidence intervals for all the regression parameters from the three methods: MEGH-I, MEGH-II and MLE, for the case when the true mixed structure is MEGH-I and the random effects are generated from a normal distribution with ${\sigma }_u=1$ .

Parameter	True value	Method	$95\text{\%}$ confidence interval
		MLE	$(0.0880,0.0933)$
$\eta$	$0.20$	MEGH-I	$(0.1999,0.2089)$
		MEGH-II	$(0.2815,0.3017)$
		MLE	$(1.6734,1.7146)$
$\nu$	$1.50$	MEGH-I	$(1.5012,1.5244)$
		MEGH-II	$(1.3492,1.3762)$
		MLE	$(5.4804,5.7141)$
$\delta$	$3.00$	MEGH-I	$(2.9988,3.0882)$
		MEGH-II	$(2.3658,2.4574)$
		MLE	$(0.9812,1.0106)$
$\alpha$	$0.96$	MEGH-I	$(0.9307,0.9643)$
		MEGH-II	$(1.1294,1.1829)$
		MLE	$(0.9847,1.0010)$
${\beta }_1$	$1.00$	MEGH-I	$(0.9951,1.0054)$
		MEGH-II	$(0.9676,0.9782)$
		MLE	$(0.0596,0.0785)$
${\beta }_2$	$0.08$	MEGH-I	$(0.0741,0.0919)$
		MEGH-II	$(0.0713,0.0891)$
		MLE	$(0.1540,0.1633)$
${\beta }_3$	$0.22$	MEGH-I	$(0.2204,0.2289)$
		MEGH-II	$(0.2172,0.2258)$
		MLE	$(0.0641,0.0872)$
${\beta }_4$	$0.10$	MEGH-I	$(0.0914,0.1009)$
		MEGH-II	$(0.0901,0.1001)$

Table 2 shows the average AIC of the model fit across the $250$ simulation replications for the three models considered. The results indicate that the MEGH model with the correct mixed hazard structure has the smallest AIC, while the model ignoring random effects tends to produce a very large AIC. Table 2 also presents the average power of the likelihood ratio test for random effects across the $250$ simulation replications, when the random effects are generated from a normal distribution with ${\sigma }_u=1$ . It can be seen that the test provides a high power with this sample size. Further simulation results on the performance of the test for random effects under different number of clusters, censoring rates and variance values can be found in the Supplemental Material.

Table 2.

Average AIC of the model fit and average power of the test for random effects across $250$ simulation replications, when the random effects are generated from a normal distribution with ${\sigma }_u=1$ .

True structure	Fitted model	AIC	Power
	MLE	$1390.66$	$-$
MEGH-I	MEGH-I	$927.84$	$1.0$
	MEGH-II	$975.79$	$1.0$
	MLE	$1383.18$	$-$
MEGH-II	MEGH-I	$996.35$	$1.0$
	MEGH-II	$962.57$	$1.0$

6 Real data application

In this section, we apply the MEGH model to analyse the data set LeukSurv, which contains information on the survival of $n=1043$ patients of acute myeloid leukaemia in northwest England, recorded between 1982 and 1998. This data set is available in the R package spBayesSurv. ¹⁷ Previous analyses of this data set have suggested that there is evidence of variation in survival across this region, ¹⁷ and we can see that this is also suggested by the nonparametric Kaplan-Meier estimators of the survival curves by district shown in Figure 4. We use information about the survival time (in years), vital status at the end of follow-up (0 – right-censored, 1 – dead), age (in years), sex (0 – female, 1 – male), white blood cell count at diagnosis (wbc, truncated at 500), the Townsend score (tpi, higher values indicates less affluent areas), and administrative district of residence (district, 24 districts). We fit the MEGH models with hazard structures MEGH-I (2) and MEGH-II (3), and with the PGW and log-logistic baseline hazards. For the random-effects distribution, we use the normal, Student- $t$ and two-piece normal distributions (to account for heavier tails than normal and asymmetry). In all fitted models, we consider the time dependent effect ( ${\tilde{\mathbf{x}}}_{ij}$ ) of age (standardised), and the hazard-level effects ( ${\mathbf{x}}_{ij}$ ) of age (standardised), sex, white blood cell count at diagnosis (wbc, standardised) and the Townsend score (tpi, standardised). The variable district is used to define the random effects in the two mixed hazard structures considered. In addition, we fit the models ignoring random effects for comparison. We also assess the need for including random effects using the diagnostic tests in Section 4.2.

The best model selected using AIC (1553.725) is the model with the mixed structure MEGH-I, log-logistic baseline hazard (with log-location and scale parameters $({\theta }_1,{\theta }_2)=(\mu ,\tau )$ ) and normal random effects. The AIC for the model ignoring random effects with log-logistic baseline hazard is 1556.366, which is fairly close to the AIC value for the best model. The reason for this is that the estimate of the variance of random effects is relatively small ( ${\hat{\sigma }}_u=0.144$ ), and the number of clusters ( $r=24$ ) is not large in this data set. However, the p-value for testing $H_0:{\sigma }_u=0$ is $0.0156$ , which suggests that the random effect is significant and must be present in the model. This implies that there is significant between-cluster variability after incorporating the information on the observed covariates (see also the Box and Whisker plot in the online Supplemental Material). Figure 5 shows the maps of aggregated summaries at district level. We notice from Figure 5(c) and (d) that the marginal survival functions (which is obtained as the average of cluster-specific marginal survival functions using Monte Carlo integration) at $t=1,5$ years exhibit some variability by district. Figure 5(a) and (b) shows that the distribution of mean age and mean Townsend score is also quite different for the different districts. Consequently, survival gaps between different districts are explained by a combination of complex factors, including different age and Townsend distribution, a conclusion that could be harder to obtain with fully nonparametric methods that require data stratification.

Figure 5.

Leukaemia data: Maps showing the aggregated summaries for each district: (a) mean age, (b) mean Townsend score, (c) marginal 1-year survival and (d) marginal 5-year survival.

The gradient function plot for the model MEGH-I fitted under normal random effects is presented in Figure 6(a), together with the $95\text{\%}$ confidence bands similar to Verbeke and Molenberghs, ³¹ which indicates that the normality assumption on random effects is valid since the gradient values remain under or close to $1$ . To double check whether the fit can be improved by using a distribution for random effects which allows to capture asymmetry, we calculate the gradient function plot with the two-piece normal distribution for random effects. The gradient plot with the $95\text{\%}$ confidence bands for the two-piece normal distribution, shown in Figure 6(b), is very similar, suggesting that the usual normal distribution is adequate for random effects.

Figure 6.

Leukaemia data: (a) Gradient function plot with the $95\text{\%}$ confidence bands for the MEGH model (2) with normal random effects fitted to the leukaemia data. (b) Gradient function plot with the $95\text{\%}$ confidence bands for the MEGH model (2) with a two-piece normal distribution for the random effects.

7 Discussion

We have introduced a general mixed-effects hazard structure (MEGH), as well as two general subclasses of interest (MEGH -I and MEGH-II). The proposed structures generalise classical survival regression models of interest in practice, such as the MEPH and MEAFT. We have adopted a flexible parametric approach, which allowed us to prove asymptotic results under standard regularity conditions. Our simulation studies show that the estimates of the parameters of the proposed models, based on the marginal likelihood, have good frequentist properties. Another contribution of this work consists of evaluating the impact of misspecifying the role of the random effects and potentially the random-effects distribution. Our simulation study shows that not only misspecifying the distribution of the random effects has an impact on the quality of the inference on the parameters, but also misspecifying the role of the random effects in the hazard regression model, a topic that has received little attention in survival analysis. The routines used to produce our results are implemented in the R package MEGH which is available, along with the real data application, at https://github.com/FJRubio67/MEGH. In our simulations and applications, we have focused on the use of the MEGH-I and MEGH-II structures, as these represent generalisations of the most common mixed models in survival analysis. However, the implementation of the general structure (1), with two random effects, represents a possible extension of this work. Such extension does not only bring computational challenges, as one requires double numerical integration, but also the modelling of the dependence between the two random effects. The latter can be explored by using copulas to capture different types of dependencies, and a variety of parametric marginal distributions. ⁸

The tractability of the MEGH model allows for several extensions, such as accounting for left-censored or interval censored times. A more substantial extension corresponds to the case where the distribution of the random effects reflects a spatial multilevel structure. ⁴¹ Another possible extension consists of modelling the baseline hazard using nonparametric (frequentist or Bayesian) methods. ²⁵ Extensions to the inclusion of more than one hierarchical level would also be of practical interest. This could be done simply by considering multivariate random effects

$$h(t_{ij}\mid {\mathbf{x}}_{ij},{\mathbf{u}}_i,{\tilde{\mathbf{u}}}_i)=h_0\left(t_{ij}\exp \{{\tilde{\mathbf{x}}}_{ij}^{\mathrm{\top }}\mathit{\alpha }+{\tilde{\mathbf{z}}}_{ij}^{\mathrm{\top }}{\tilde{\mathbf{u}}}_i\}\right)\exp \{{\mathbf{x}}_{ij}^{\mathrm{\top }}\mathit{\beta }+{\mathbf{z}}_{ij}^{\mathrm{\top }}{\mathbf{u}}_i\}$$

where ${\mathbf{z}}_{ij}$ and ${\tilde{\mathbf{z}}}_{ij}$ are random effects covariates associated with the hazard and time scales respectively, and $({\mathbf{u}}_i,{\tilde{\mathbf{u}}}_i)\stackrel{iid}{\sim }{G}_{q+\tilde{q}}$ , where $G_{q+\tilde{q}}$ is a continuous distribution with support on ${\mathbb{R}}^{q+\tilde{q}}$ and zero mean. We emphasise that, although interesting, this strategy requires a careful study of the identifiability of parameters. Our formulation remains valid for these extensions, but the main challenges relate to the development of computational methods for addressing those questions.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802221102620 - Supplemental material for MEGH: A parametric class of general hazard models for clustered survival data

Supplemental material, sj-pdf-1-smm-10.1177_09622802221102620 for MEGH: A parametric class of general hazard models for clustered survival data by Francisco Javier Rubio and Reza Drikvandi in Statistical Methods in Medical Research