Multivariate generalized linear frailty models for clustered competing risk data

Multivariate generalized linear frailty models for clustered competing risk data Ren Teranishi1, Kyoji Furukawa2, Takeshi Emura3 1Kurume University Graduate School of Medicine, Kurume, Japan 2Biostatistics Center, Kurume University, Kurume, Japan 3School of Informatics and Data Science, Hiroshima University, Hiroshima, Japan Abstract Clustered competing risk data occur when individuals within clusters are subject to multiple mutually exclusive event types, inducing dependence both across events and within clusters. We propose a joint modeling framework for such data based on a multivariate generalized linear frailty model. Cause-specific hazards are specified in a piecewise exponential form with shared clusterlevel random effects to account for unobserved heterogeneity and between-event correlation. Using the equivalence between the piecewise exponential and Poisson regression likelihoods, estimation is performed under the generalized linear mixed model (GLMM) framework, allowing implementation with standard mixed-model software. Simulation studies show that the proposed method yields nearly unbiased and efficient estimation across a wide range of correlation structures, whereas conventional univariate and Cox-type frailty models exhibit bias or instability under moderate dependence. Application to multicenter clinical trial data illustrates the practical utility and interpretability of the proposed model. The approach offers a flexible and extensible framework for modeling clustered survival data with competing risks. Key Words: Competing risks, cluster survival data, generalized linear mixed models, Poisson regression 1. Introduction In survival analysis, clustered data frequently arise when individuals are grouped within higherlevel units such as families, hospitals, or clinical centers, resulting in correlated event times within clusters (1). Standard survival models that assume independence between individuals may yield biased estimates and underestimated variances (2). The presence of competing risks further complicates the analysis, since individuals may experience one of several mutually exclusive event types that preclude the occurrence of the other types of events (3, 4). Ignoring the dependence induced by clustering or the presence of competing risks can result in misleading inference and poor predictive performance (5). Therefore, flexible statistical models that simultaneously account for both clustering and competing risks are needed for valid analysis in biomedical and epidemiological research. Several statistical approaches have been proposed for analyzing clustered competing risk data. A common approach is to use marginal models, such as extensions of the Cox proportional hazards model to incorporate random effects for within-cluster correlation (1). Another approach introduces random effects, or frailties, to explicitly capture unobserved heterogeneity across clusters. Frailty models are particularly appealing because they allow for direct modeling of the dependence structure and can be extended to multivariate settings to accommodate multiple competing event types simultaneously (6). In the competing risks framework, two main modeling perspectives are often considered: cause-specific hazard models, which separately model the hazard of each event type, and subdistribution hazard models, which focus on the cumulative incidence of events (7). Incorporating frailties into these frameworks provides a flexible means of accounting for both clustering and competing risks, but this requires careful formulation to handle the joint dependence between event types as well as the correlation that arise from clustering and competing events (2, 5). In this study, we propose a joint modeling framework for the analysis of clustered competing risk data that simultaneously accommodates within-cluster dependence and the presence of multiple event types. Specifically, the hazard for each competing risk event is assumed to follow a piecewise exponential form, with a shared frailty introduced to capture unobserved heterogeneity among individuals within the same cluster. This shared frailty induces correlation not only across individuals but also across competing events within a cluster, providing a coherent multivariate structure. By exploiting the equivalence between the piecewise exponential model and the Poisson regression likelihood, the proposed approach can be conveniently implemented within the generalized linear mixed model (GLMM) framework. This connection allows the use of standard mixed-model methodology and software, while offering flexibility to extend the model to multivariate settings and to incorporate covariate effects on different types of event hazards. The remainder of this paper is structured as follows. Section 2 provides an overview of clustered competing risk data and outlines the key statistical challenges arising from within-cluster dependence and multiple event types. Section 3 presents the proposed multivariate generalized linear frailty model, detailing its formulation, likelihood construction, and estimation procedures under the GLMM framework. Section 4 presents a simulation study conducted to assess the bias, efficiency, and robustness of the proposed method relative to existing approaches. Section 5 illustrates the practical application of the method using data from a multicenter clinical trial on bladder cancer, demonstrating its interpretability and computational feasibility. Finally, Section 6 discusses the implications of the findings, summarizes the contributions, and highlights potential avenues for future methodological development. 2. Clustered competing risk data 2.1. Clustered competing risk data In many survival analyses, event times are correlated because individuals are grouped into clusters. Such clustering can arise when, for example, subjects in the same group, such as hospital or family, share unmeasured factors affecting the risk of disease of interest. Presence of such unmeasured shared factors can introduce correlation among the survival times of individuals within the same cluster if they are unobserved or not included in the analysis model(8). Ignoring such within-cluster dependence can result in underestimated variances and overly optimistic conclusions about the precision of covariate effects(9). Competing risks arise when individuals are simultaneously at risk of experiencing one of several mutually exclusive events, where the occurrence of one event precludes the occurrence of the others. For example, in studies of cancer patients, death due to the primary cancer and death due to other causes (such as cardiovascular disease) represent competing events. Ignoring competing risks and treating occurrence of alternative events as censored can bias estimates of the hazard for the event of interest and misleading interpretation of survival and covariate effects. When clustering and competing risks occur together, the analysis becomes more complex because both within-cluster correlation and the presence of multiple event types must be addressed simultaneously(8). For instance, in multicenter cancer studies, patients treated at the same hospital creates clustering, while different causes of death—such as cancer progression or treatment-related complications—constitute competing risks. Individuals within the same cluster may share unmeasured factors that influence the risks of all competing events, such as institutional practices and clinical characteristics of a hospital. This induces dependence not only in overall survival times but also across different types of events. Ignoring this structure can lead to biased estimation of cause-specific or subdistribution hazards and incorrect standard errors, ultimately compromising statistical inference. Proper modeling must therefore account for both the competing nature of events and the correlation induced by clustering to obtain valid and interpretable results. 2.2. Statistical challenges in clustered competing risk data The analysis of clustered competing risk data presents methodological challenges. The first arises from the dependence structure induced by clustering. Individuals within the same cluster—such as patients treated at the same hospital or members of the same family—often share unmeasured factors influencing their risks for multiple event types. This shared heterogeneity creates correlations both among individuals and between competing events for the same individual. Ignoring such dependence leads to underestimated variances, inflated type I error rates, and biased inference(9). A second challenge concerns the distinction between marginal and conditional inference. Marginal approaches, such as cause-specific Cox models with robust variance estimators, provide population-averaged effects but cannot explicitly capture cluster-level heterogeneity. In contrast, frailty-based conditional models directly incorporate latent random effects, allowing inference on cluster-specific risks. In clustered competing risks, the choice between these frameworks affects both interpretation and model adequacy(8). A third issue is bias due to model misspecification. When clustering or competing risks are ignored, regression coefficients and hazard functions may be distorted, leading to misleading conclusions about treatment or prognostic effects. Covariates may also exert different influences on distinct event types, and their effects may vary across clusters, further complicating analysis. Building on these considerations, we propose a unified joint modeling framework for clustered competing risk data that simultaneously accounts for within-cluster dependence and cross-event correlation. By formulating cause-specific hazards under a piecewise exponential model and exploiting their equivalence to Poisson regression likelihoods, the framework can be implemented within the generalized linear mixed model (GLMM) structure. This approach offers both conceptual clarity and computational convenience, allowing flexible modeling of complex survival data using standard statistical software. 3. Methods This section introduces the proposed multivariate generalized linear frailty model for analyzing clustered competing risk data. This approach exploits the equivalence between piecewise exponential survival models and Poisson regression likelihoods, allowing estimation within the generalized linear mixed model (GLMM) framework. We first describe the discretization of survival times and its connection to Poisson regression, then extend the model to jointly analyze multiple competing event types with cluster-level random effects. 3.1. Discretized survival data analysis using Poisson regression Let 𝑇𝑖 denote the survival time and 𝑑𝑖 the event indicator for individual 𝑖 ∈{1,…,𝑛}. To model the hazard flexibly, the follow-up period is divided into 𝐽 disjoint intervals [𝜏𝑗−1,𝜏𝑗),𝑗=1,..,𝐽. Assuming a constant hazard in each interval yields the piecewise exponential hazard model; ℎ𝑖(𝑡)=𝜆𝑖𝑗 =exp(𝛼𝑗+𝒙𝒊 ′𝜷) ∀𝑡∈[𝜏𝑗−1,𝜏𝑗)， 𝑗=1,…,𝐽 where 𝒙𝑖 ′=(𝑥𝑖1,…,𝑥𝑖𝑝) is a vector of covariates, 𝜷′=(𝛽1,…,𝛽𝑝) is the regression coefficients, and 𝛼𝑗 is the baseline log-hazard for interval 𝑗. Under this formulation, each contributes one record for each interval during which he/she is at risk. Let 𝑡𝑖𝑗 denote the observed duration at risk and 𝑑𝑖𝑗 the event indicator for subject 𝑖 in interval 𝑗. Given the observed data {𝑑𝑖𝑗,𝑡𝑖𝑗;;𝑖=1,…,𝑛,𝑗=1,…,𝐽}, the log-likelihood for parameters is log𝐿=∑∑{𝑑𝑖𝑗log𝜆𝑖𝑗 −𝜆𝑖𝑗𝑡𝑖𝑗} 𝐽 𝑗=1 𝑛 𝑖=1 . This expression is essentially equivalent to the log-likelihood of a Poisson regression model for the count 𝑑𝑖𝑗 with mean 𝜆𝑖𝑗𝑡𝑖𝑗. This equivalence provides an intuitive and computationally convenient approach to model survival data within the standard generalized linear model (GLM) framework. Figure 1 illustrates a simple example of constructing a discrete survival dataset. Suppose we have survival data from 25 patients (Figure 1a). By dividing the follow-up period into time intervals defined with cut points at 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0, the data can be expanded and aggregated accordingly (Figure 1b). A Poisson regression model can then be fitted to estimate the hazard function in a piecewise constant form (dashed lines in Figure 1c). Alternatively, a smooth parametric function can be fitted using the midpoints (red curve in Figure 1c). This intuitive representation of the hazard is a key advantage of this approach. 𝑑 𝑖 𝑡𝑖 𝑖 1 2.8 1 1 2.1 2 0 3.0 3 1 0.9 4 0 3.0 5 1 2.7 6 1 2.6 7 1 2.3 8 1 1.7 9 1 1.3 10 0 3.0 11 1 1.9 12 1 2.5 13 1 0.6 14 0 3.0 15 1 2.7 16 1 0.8 17 1 1.4 18 1 2.3 19 1 1.2 20 0 3.0 21 1 1.3 22 1 0.4 23 1 1.7 24 0 3.0 25 (a) Survival data (n=25) (b) Discrete Survival data (c) Fitting hazard models 0 1 2 3 Hazard rate Total events Total time Interval 𝜏𝑗−1,𝜏𝑗 𝑗 0.081 112.4[0,0.5) 1 0.265 311.3 [0.5,1.0) 2 0.412 49.7 [1.0,1.5) 3 0.385 37.8 [1.5,2.0) 4 0.645 46.2 [2.0,2.5) 5 1.053 43.8 [2.5,3.0) 6 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 time hazard rate 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 time hazard rate 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 time hazard rate Piecewise-constant midpoints Parametric event Figure 1. A simple illustration of constructing a discrete survival dataset and fitting hazard models. 3.2. Extension to multi-endpoints model with cluster To jointly account for clustering and competing risks, the above formulation is extended to a multivariate setting. Let 𝑇𝑖𝑘 denote the survival time to event type 𝑘 and 𝑑𝑖𝑘 the event indicator for individual 𝑖 in cluster 𝑐(𝑖), 𝑖∈{1,…,𝑛}, 𝑘∈{1,…,𝐾}, 𝑐(𝑖)∈{1,…,𝐶}. For each event type 𝑘, a piecewise-exponential hazard model is specified as logℎ𝑘𝑖(𝑡)=𝛼𝑘𝑗 +𝒙𝑖 ′𝜷𝑘+𝜂𝑘𝑢𝑐(𝑖) ∀𝑡∈[𝜏𝑗−1,𝜏𝑗)， 𝑗∈{1,…,𝐽}, where 𝑢1,…,𝑢𝐶 ~𝑖𝑖𝑑 𝑁(0,1) are cluster-level random effects. Parameters {𝜂1,…,𝜂𝐾} represent event-specific random effects variation, reflecting the strength of within-cluster correlation for each competing event. For simplicity, a common set of covariates 𝒙𝑖 ′=(𝑥𝑖1,…,𝑥𝑖𝑝) is assumed across event types, although in practice different events may have distinct covariate associations. For the bivariate case (𝐾=2), the model may be described as logℎ1𝑖(𝑡)=𝜆1𝑖𝑗 =𝛼1𝑗 +𝒙𝑖 ′𝜷1+𝜂1𝑢𝑐(𝑖) logℎ2𝑖(𝑡)=𝜆2𝑖𝑗 =𝛼2𝑗 +𝒙𝑖 ′𝜷2+𝜂2𝑢𝑐(𝑖) ∀𝑡∈[𝜏𝑗−1,𝜏𝑗)， 𝑗∈{1,…,𝐽}. Given the observed data {𝑑𝑘𝑖𝑗,𝑡𝑘𝑖𝑗;𝑖=1,…,𝑛,𝑗=1,…,𝐽,𝑘=1,2},𝑡he likelihood for estimating parameters 𝜽={𝜶,𝜷,𝜼}, 𝜶={(𝛼1𝑗,𝛼2𝑗);𝑗=1,…𝐽},𝜷=(𝜷1,𝜷2), 𝜼=(𝜂1,𝜂2) is log𝐿= ∫∏∏(𝜆1𝑖𝑗)𝑑1𝑖𝑗 exp(−𝜆1𝑖𝑗𝑡1𝑖𝑗) 𝐽 𝑗=1 𝑛 𝑖=1 ×(𝜆2𝑖𝑗)𝑑2𝑖𝑗 exp(−𝜆2𝑖𝑗𝑡2𝑖𝑗)𝑓(𝑢) ∞ −∞ 𝑑𝑢 where 𝑓(𝑢) denotes the probability density of the frailty distribution. This unified formulation allows simultaneous modeling of cause-specific hazards while properly accounting for within-cluster dependence and cross-event correlation. The model can be estimated using standard GLMM software—for example, the glmer() function in the lme4 package in R (10)— facilitating implementation in practical biomedical studies. 3.3. Estimation and inference Because the marginal likelihood does not have a closed-form solution, estimation relies on numerical integration, such as adaptive Gaussian quadrature or the Laplace approximation, both available in standard GLMM routines. Parameter estimation proceeds by maximizing the approximated log-likelihood using iterative optimization algorithms (e.g., Newton–Raphson). Standard errors for the regression coefficients 𝜷 and variance components 𝜼 are obtained from the inverse of the observed information matrix or, for robustness, via sandwich variance estimation. Wald-type confidence intervals and hypothesis tests are then constructed under asymptotic normality. The proposed framework provides consistent and efficient inference for covariate effects while properly accounting for within-cluster correlation and event-type dependence. Model adequacy can be assessed using likelihood-based criteria such as the Akaike Information Criterion (AIC) or likelihood ratio tests for the inclusion of frailty terms. 4. Simulation study A simulation study was conducted to evaluate the performance of the proposed multivariate generalized linear frailty model under clustered competing risk settings. The primary objectives were to assess estimation accuracy, robustness to within-cluster dependence, and comparative performance relative to existing Cox-based methods. 4.1. Simulation model Data were generated from a bivariate clustered survival model with two competing event types. The hazard functions for the two endpoints are assumed to have the following forms: logℎ1𝑖(𝑡1)=logℎ01(𝑡1)+𝛼𝑥𝑖+𝜎𝑢𝑐(𝑖) logℎ2𝑖(𝑡2)=logℎ02(𝑡2)+𝛽𝑥𝑖+𝜂𝜎𝑢𝑐(𝑖) 𝑖=1,…,𝑛, where 𝑛=50 or 200. A common Weibull baseline hazard was separately assumed for the two endpoints: ℎ01(𝑡)=ℎ02(𝑡)=𝜆𝛾𝑡𝛾−1 where 𝜆=0.001, 𝛾=3. There was a binary covariate (𝑥) with a true relative risk of 1.5. A half of the subjects were assigned 𝑥=1 and the rest 𝑥 =0. The corresponding log-relative risk parameter for the second event (𝛽) was the main estimand of interest. The number of clusters was assumed to be 20, and cluster-level random effects 𝑢1,…,𝑢20 were generated independently from the standard normal distribution. Two parameters, 𝜎(= 0, 0.1,…,or 2.0) and 𝜂(=0.5,.1.25, or 1.5), were assumed to control the within-cluster correlations for the two competing events. Specifically, 𝜎 controls the magnitude of cluster-level heterogeneity for event 1, while 𝜂 determines the relative heterogeneity for event 2 relative to the heterogeneity for event 2, determining the correlation strength between the competing events within each cluster (if 𝜎=0, the two endpoints are independent). Censoring times were independently drawn from a uniform distribution on (0,10), resulting in an average censoring rate of approximately 20%. 4.2. Analysis models Each simulated dataset was analyzed using the following models: 1. Univariate models ignoring the competing events: • Cox proportional hazards regression model • Piecewise exponential regression with Poisson likelihood formulation 2. Bivariate models accounting for clustering and competing risks: • Cox frailty model with shared gamma frailty across endpoints • Proposed multivariate generalized linear frailty model using the Poisson likelihood formulation For each method, estimates of 𝛽 were obtained, and performance was evaluated in terms of empirical bias and mean squared error (MSE) across the 500 simulation replications. 4.3. Results (a) Results for 𝑛=200 Figure 2a shows the empirical bias for estimation of 𝛽 across varying levels of frailty variation, determined by 𝛽 and 𝜂. When 𝜎=0, the competing events were independent, and all methods yielded negligible bias. As 𝜎 increased, univariate Cox and piecewise exponential models (PWE), which ignore the competing event, exhibited substantial downward bias. The bivariate Cox frailty model reduced this bias but showed upward bias when 𝜂=0.5, corresponding to weaker correlation for the primary event compared to the competing event. In contrast, the multivariate Poisson frailty model (Proposed) remained essentially unbiased across all parameter settings. Table 1a summarizes the results for MSE. As 𝜎 increased, estimation precision declined for all models, whereas the proposed approach consistently achieved the lowest MSE, demonstrating superior efficiency and stability. (b) Results for 𝑛=50 Figure 2b and Table 1b present results with a smaller sample size (𝑛=50), where differences among models became more pronounced. The Cox frailty model frequently encountered convergence issues and unstable estimation of frailty variance, leading to inflated bias and MSE. In contrast, the proposed Poisson frailty model maintained good estimation accuracy and efficiency even with limited number of clusters and events. The model provided stable estimates of 𝛽 and variance components, demonstrating robustness in small-sample settings. 1 𝜂=0.5 𝜂=1.25 𝜂=1.5 (a) 𝑛=200 (b) 𝑛=50 Figure 2. Empirical bias in estimating 𝛽 with four models: Cox proportional hazards (Cox), piecewise exponential (PWE), Cox frailty, and the proposed multivariate Poisson frailty model (Proposed) under varying degrees of within-cluster heterogeneity and cross-event correlation (determined by 𝜎 and 𝜂). Table 1. Empirical mean squared error (MSE) of 𝛽 with four models: Cox proportional hazards (Cox), piecewise exponential (PWE), Cox frailty, and the proposed multivariate Poisson frailty model (Proposed) under varying degrees of within-cluster heterogeneity and cross-event correlation (determined by 𝜎 and 𝜂). σ0 0.5 1 1.5 2 Cox 0.052 0.047 0.058 0.062 0.077 Cox frailty 0.053 0.049 0.070 0.074 0.080 PWE 0.049 0.047 0.056 0.061 0.076 Proposed 0.050 0.047 0.058 0.057 0.063 σ0 0.5 1 1.5 2 Cox 0.056 0.058 0.073 0.101 0.104 Cox frailty 0.057 0.056 0.052 0.055 0.056 PWE 0.054 0.059 0.075 0.104 0.109 Proposed 0.056 0.057 0.053 0.059 0.058 σ0 0.5 1 1.5 2 Cox 0.056 0.061 0.080 0.103 0.126 Cox frailty 0.056 0.053 0.051 0.045 0.047 PWE 0.054 0.062 0.082 0.108 0.132 Proposed 0.056 0.055 0.056 0.047 0.048 1 σ0 0.5 1 1.5 2 Cox 0.265 0.254 0.250 0.271 1.029 Cox frailty 0.179 0.266 0.607 1.180 1.798 PWE 0.246 0.240 0.235 0.260 0.269 Proposed 0.269 0.256 0.264 0.276 0.296 σ0 0.5 1 1.5 2 Cox 0.252 0.282 0.259 0.271 0.281 Cox frailty 0.179 0.413 1.165 2.459 3.921 PWE 0.229 0.271 0.256 0.276 0.280 Proposed 0.239 0.274 0.250 0.285 0.269 σ0 0.5 1 1.5 2 Cox 0.253 0.218 0.263 0.332 0.318 Cox frailty 0.185 0.508 1.467 2.770 4.958 PWE 0.239 0.207 0.259 0.320 0.320 Proposed 0.257 0.233 0.266 0.326 0.255 (a) 𝑛=200 (b) 𝑛=50 𝜂=0.5 𝜂=1.25 𝜂=1.5 𝜂=0.5 𝜂=1.25 𝜂=1.5 5. Real data application To illustrate the utility of the proposed multivariate generalized linear frailty model, we analyzed data from a clinical trial of bladder cancer patients. The aim was to evaluate the effects of chemotherapy and age on two clinically relevant competing events — recurrence of bladder cancer and death before recurrence — while accounting for potential correlation among patients treated at the same medical facility. 5.1. Dataset The data were obtained from the EORTC 30791 trial (11), which enrolled 392 patients diagnosed with stage Ta–T1 bladder cancer across 21 participating medical facilities. We shall use the subset of the dataset as considered in Section 1.2.4 of Ha et al. (12), which is available as the “bladder” dataset from the R package frailtyHL (13). The study followed patients prospectively from treatment initiation until the occurrence of either the first recurrence of bladder cancer or death before recurrence, whichever came first. These two types of events were treated as competing risks, since the occurrence of one precluded observation of the other. • Event 1: Death before recurrence (81 cases; 21%). • Event 2: First recurrence of bladder cancer (200 cases; 51%). A total of 111 patients (28.3%) were censored during follow-up due to withdrawal or the end of follow-up without an event. Cluster sizes varied markedly, ranging from 3 to 178 patients per facility, reflecting differences in institutional enrollment rates and patient loads. Two covariates were included in the analysis: • Chemotherapy (CHEMO): coded as 1 for patients receiving chemotherapy and 0 otherwise. • Age group: coded as 1 for patients aged ≥ 66 years, and 0 otherwise. The primary objective was to quantify the effect of chemotherapy on the risk of recurrence and the effect of age on mortality before recurrence, while accounting for between-facility heterogeneity through shared random effects. 5.2. Models We analyzed the data using the proposed multivariate generalized linear frailty model under the piecewise exponential formulation, implemented as a generalized linear mixed model (GLMM). Separate submodels of a common form were specified for each the two competing event type, allowing event-specific regression parameters and baseline hazards. Let 𝑇𝑖𝑘 denote the observed time to event for patient 𝑖 in facility 𝑐(𝑖), and let 𝑘 ∈{1,2} index the competing events (1=death before recurrence and 2=first recurrence). The cause-specific hazard functions were modeled as Event-1(death): logℎ1𝑖(𝑡1)=logℎ01(𝑡1)+𝛼1𝐴𝐺𝐸𝑖+𝛽1𝐶𝐻𝐸𝑀𝑂𝑖+𝜂1𝑢𝑐(𝑖) Event-2(recurrence): logℎ2𝑖(𝑡2)=logℎ02(𝑡2)+𝛼2𝐴𝐺𝐸𝑖+𝛽2𝐶𝐻𝐸𝑀𝑂𝑖+𝜂2𝑢𝑐(𝑖) where 𝑢𝑐(𝑖) is a frailty term for the facility that subject 𝑖 belongs to. The frailties were assumed to follow independent normal distributions, 𝑢1,…,𝑢𝐶~𝑁(0,1), where 𝐶=21 is the total number of facilities. The baseline hazards for each event type were modeled flexibly using cubic splines with degrees of freedom (df) of 3, which is expected to allow smooth, data-driven estimation of timedependent hazard patterns: logℎ0𝑘(𝑡)=𝜂𝑘0 +∑𝜉𝑘𝑟𝐵𝑟(𝑡(𝑗)) df 𝑟=1 ,𝑘=1,2, for 𝑡∈[𝜏𝑗−1,𝜏𝑗), where 𝑡(𝑗) denotes the midpoint of the 𝑗-th interval. Estimation was performed using the Poisson likelihood representation of the piecewise exponential model within the GLMM framework. Computation was implemented by the glmer() function in the R package lme4, which uses the Laplace approximation to the marginal likelihood. All models converged successfully without numerical instability. For comparison, Cox frailty models with normally distributed shared frailties were fitted using the same covariates and cluster structure. This allowed direct assessment of the proposed model against a conventional semiparametric approach. 5.3. Results Parameter estimates and their corresponding 95% confidence intervals for both endpoints (death before recurrence and recurrence) are summarized in Table 2. Across all models, results were highly consistent. For Event 2 (recurrence of bladder cancer), chemotherapy significantly reduced the risk of recurrence. The estimated hazard ratio for chemotherapy ranged between 0.51 and 0.53 across models, corresponding to a 47-49% reduction in the hazard of recurrence compared to the nonchemotherapy group. The estimated treatment effect remained robust across both the proposed Poisson frailty model and the Cox frailty model, suggesting that the choice of modeling framework had little influence on inference regarding treatment efficacy. For Event 1 (death before recurrence), age was the dominant prognostic factor. Patients aged 66 years or older exhibited approximately 1.8–2.0 times higher hazard of death relative to younger patients, after adjusting for treatment and cluster effects. The effect of chemotherapy on death was not statistically significant, consistent with the trial’s clinical findings that the treatment primarily reduces recurrence risk rather than mortality. The estimated variance of the shared frailty terms (𝜂1, 𝜂2) were modest but non-zero, indicating moderate between-facility heterogeneity for both endpoints. The inclusion of frailty improved model fit as measured by the Akaike Information Criterion (AIC), confirming the presence of facility-level clustering effects. Figure 3 displays the estimated baseline hazards for both events obtained from the spline-based baseline hazard model (degrees of freedom = 3). The estimated hazard for recurrence (Event 2) peaked within the first two years after treatment and gradually declined thereafter, while the hazard for death (Event 1) remained relatively stable over time. These temporal patterns are consistent with the natural history of early-stage bladder cancer, where recurrence risk is highest shortly after initial treatment, while mortality risk accumulates more slowly. Overall, the proposed model provided stable estimates, efficient computation, and interpretable parameter estimates. No convergence issues were observed, and results were in close agreement with those obtained from the Cox frailty models. Table 2. Parameter estimates from the multivariate frailty analysis of the EORTC 30791 bladder cancer trial. Estimated regression coefficients, hazard ratios (HR), and 95% confidence intervals (CI) for age and chemotherapy effects under two competing events: (1) death before recurrence and (2) first recurrence of bladder cancer. Figure 3. Estimated baseline hazard functions for recurrence (Event 1) and death before recurrence (Event 2) in the EORTC 30791 trial. Spline-based baseline hazard estimates (degrees of freedom = 3) obtained from the proposed multivariate Poisson frailty model. 6. Discussion and summary This study introduced a unified statistical framework for analyzing clustered competing risk data by extending the piecewise exponential survival model into a multivariate generalized linear frailty model. The key innovation lies in exploiting the equivalence between the piecewise exponential likelihood and Poisson regression, thereby embedding complex survival structures within the generalized linear mixed model (GLMM) framework. This formulation offers a flexible and computationally efficient means of jointly modeling within-cluster dependence and cross-event correlation using standard mixed-model software. The proposed model provides several methodological and practical advantages over existing approaches. First, by incorporating shared random effects across competing risks, it explicitly