Bayesian inference for recurrent events data using time-dependent frailty

In medical studies, we commonly encounter multiple events data such as recurrent infection or attack times in patients suffering from a given disease. A number of statistical procedures for the analysis of such data use the Cox proportional hazards model, modified to include a random effect term called frailty which summarizes the dependence of recurrent times within a subject. These unobserved random frailty effects capture subject effects that are not explained by the known covariates. They are typically modelled constant over time and are assumed to be independently and identically distributed across subjects. However, in some situations, the subject-specific random frailty may change over time in the same manner as time-dependent covariate effects. This paper presents a time-dependent frailty model for recurrent failure time data in the Bayesian context and estimates it using a Markov chain Monte Carlo method. Our approach is illustrated by a data set relating to patients with chronic granulomatous disease and it is compared to the constant frailty model using the deviance information criterion.     

Nonproportional hazards and unobserved heterogeneity in clustered survival data: When can we tell the difference?

Multivariate survival data are frequently encountered in biomedical applications in the form of clustered failures (or recurrent events data). A popular way of analyzing such data is by using shared frailty models, which assume that the proportional hazards assumption holds conditional on an unobserved cluster-specific random effect. Such models are often incorporated in more complicated joint models in survival analysis. If the random effect distribution has finite expectation, then the conditional proportional hazards assumption does not carry over to the marginal models. It has been shown that, for univariate data, this makes it impossible to distinguish between the presence of unobserved heterogeneity (eg, due to missing covariates) and marginal nonproportional hazards. We show that time-dependent covariate effects may falsely appear as evidence in favor of a frailty model also in the case of clustered failures or recurrent events data, when the cluster size or number of recurrent events is small. When true unobserved heterogeneity is present, the presence of nonproportional hazards leads to overestimating the frailty effect. We show that this phenomenon is somewhat mitigated as the cluster size grows. We carry out a simulation study to assess the behavior of test statistics and estimators for frailty models in such contexts. The gamma, inverse Gaussian, and positive stable shared frailty models are contrasted using a novel software implementation for estimating semiparametric shared frailty models. Two main questions are addressed in the contexts of clustered failures and recurrent events: whether covariates with a time-dependent effect may appear as indication of unobserved heterogeneity and whether the additional presence of unobserved heterogeneity can be detected in this case. Finally, the practical implications are illustrated in a real-world data analysis example.     

Analyzing sickness absence with statistical models for survival data

OBJECTIVES: Sickness absence is the outcome in many epidemiologic studies and is often based on summary measures such as the number of sickness absences per year. In this study the use of modern statistical methods was examined by making better use of the available information. Since sickness absence data deal with events occurring over time, the use of statistical models for survival data has been reviewed, and the use of frailty models has been proposed for the analysis of such data. METHODS: Three methods for analyzing data on sickness absences were compared using a simulation study involving the following: (i) Poisson regression using a single outcome variable (number of sickness absences), (ii) analysis of time to first event using the Cox proportional hazards model, and (iii) frailty models, which are random effects proportional hazards models. Data from a study of the relation between the psychosocial work environment and sickness absence were used to illustrate the results. RESULTS: Standard methods were found to underestimate true effect sizes by approximately one-tenth [method i] and one-third [method ii] and to have lower statistical power than frailty models. CONCLUSIONS: An uncritical use of standard methods may underestimate the effect of work environment exposures or leave predictors of sickness absence undiscovered.     

Generalized survival models for correlated time‐to‐event data

Our aim is to develop a rich and coherent framework for modeling correlated time‐to‐event data, including (1) survival regression models with different links and (2) flexible modeling for time‐dependent and nonlinear effects with rich postestimation. We extend the class of generalized survival models, which expresses a transformed survival in terms of a linear predictor, by incorporating a shared frailty or random effects for correlated survival data. The proposed approach can include parametric or penalized smooth functions for time, time‐dependent effects, nonlinear effects, and their interactions. The maximum (penalized) marginal likelihood method is used to estimate the regression coefficients and the variance for the frailty or random effects. The optimal smoothing parameters for the penalized marginal likelihood estimation can be automatically selected by a likelihood‐based cross‐validation criterion. For models with normal random effects, Gauss‐Hermite quadrature can be used to obtain the cluster‐level marginal likelihoods. The Akaike Information Criterion can be used to compare models and select the link function. We have implemented these methods in the R package rstpm2. Simulating for both small and larger clusters, we find that this approach performs well. Through 2 applications, we demonstrate (1) a comparison of proportional hazards and proportional odds models with random effects for clustered survival data and (2) the estimation of time‐varying effects on the log‐time scale, age‐varying effects for a specific treatment, and two‐dimensional splines for time and age.     

Joint models for efficient estimation in proportional hazards regression models

In survival studies, information lost through censoring can be partially recaptured through repeated measures data which are predictive of survival. In addition, such data may be useful in removing bias in survival estimates, due to censoring which depends upon the repeated measures. Here we investigate joint models for survival T and repeated measurements Y, given a vector of covariates Z. Mixture models indexed as f (T/Z) f (Y/T,Z) are well suited for assessing covariate effects on survival time. Our objective is efficiency gains, using non-parametric models for Y in order to avoid introducing bias by misspecification of the distribution for Y. We model (T/Z) as a piecewise exponential distribution with proportional hazards covariate effect. The component (Y/T,Z) has a multinomial model. The joint likelihood for survival and longitudinal data is maximized, using the EM algorithm. The estimate of covariate effect is compared to the estimate based on the standard proportional hazards model and an alternative joint model based estimate. We demonstrate modest gains in efficiency when using the joint piecewise exponential joint model. In a simulation, the estimated efficiency gain over the standard proportional hazards model is 6.4 per cent. In clinical trial data, the estimated efficiency gain over the standard proportional hazards model is 10.2 per cent.     

Reduced-rank hazard regression for modelling non-proportional hazards

The Cox proportional hazards model is the most common method to analyse survival data. However, the proportional hazards assumption might not hold. The natural extension of the Cox model is to introduce time-varying effects of the covariates. For some covariates such as (surgical)treatment non-proportionality could be expected beforehand. For some other covariates the non-proportionality only becomes apparent if the follow-up is long enough. It is often observed that all covariates show similar decaying effects over time. Such behaviour could be explained by the popular (gamma-) frailty model. However, the (marginal) effects of covariates in frailty models are not easy to interpret. In this paper we propose the reduced-rank model for time-varying effects of covariates. Starting point is a Cox model with p covariates and time-varying effects modelled by q time functions (constant included), leading to a pxq structure matrix that contains the regression coefficients for all covariate by time function interactions. By reducing the rank of this structure matrix a whole range of models is introduced, from the very flexible full-rank model (identical to a Cox model with time-varying effects) to the very rigid rank one model that mimics the structure of a gamma-frailty model, but is easier to interpret. We illustrate these models with an application to ovarian cancer patients.     

Variable selection in subdistribution hazard frailty models with competing risks data

The proportional subdistribution hazards model (i.e. Fine‐Gray model) has been widely used for analyzing univariate competing risks data. Recently, this model has been extended to clustered competing risks data via frailty. To the best of our knowledge, however, there has been no literature on variable selection method for such competing risks frailty models. In this paper, we propose a simple but unified procedure via a penalized h‐likelihood (HL) for variable selection of fixed effects in a general class of subdistribution hazard frailty models, in which random effects may be shared or correlated. We consider three penalty functions, least absolute shrinkage and selection operator (LASSO), smoothly clipped absolute deviation (SCAD) and HL, in our variable selection procedure. We show that the proposed method can be easily implemented using a slight modification to existing h‐likelihood estimation approaches. Numerical studies demonstrate that the proposed procedure using the HL penalty performs well, providing a higher probability of choosing the true model than LASSO and SCAD methods without losing prediction accuracy. The usefulness of the new method is illustrated using two actual datasets from multi‐center clinical trials. Copyright © 2014 John Wiley & Sons, Ltd.     

Proportional hazards model with random effects

We propose a general proportional hazards model with random effects for handling clustered survival data. This generalizes the usual frailty model by allowing a multivariate random effect with arbitrary design matrix in the log relative risk, in a way similar to the modelling of random effects in linear, generalized linear and non-linear mixed models. The distribution of the random effects is generally assumed to be multivariate normal, but other (preferably symmetrical) distributions are also possible. Maximum likelihood estimates of the regression parameters, the variance components and the baseline hazard function are obtained via the EM algorithm. The E-step of the algorithm involves computation of the conditional expectations of functions of the random effects, for which we use Markov chain Monte Carlo (MCMC) methods. Approximate variances of the estimates are computed by Louis' formula, and posterior expectations and variances of the individual random effects can be obtained as a by-product of the estimation. The inference procedure is exemplified on two data sets.     

Approaches in modelling long-term survival: an application to breast cancer

Several modelling techniques have been proposed for non-proportional hazards. In this work we consider different models which can be classified into three wide categories: models with time-varying effects of the covariates; frailty models and cure rate models. We present those different extensions of the proportional hazards model on an application of 2433 breast cancer patients with a long follow-up. We comment on the differences and similarities among the models and evaluate their performance using survival and hazard plots, Brier scores and pseudo-observations.     

Testing for centre effects in multi-centre survival studies: a Monte Carlo comparison of fixed and random effects tests.

The problem of testing for a centre effect in multi-centre studies following a proportional hazards regression analysis is considered. Two approaches to the problem can be used. One fits a proportional hazards model with a fixed covariate included for each centre (except one). The need for a centre specific adjustment is evaluated using either a score, Wald or likelihood ratio test of the hypothesis that all the centre specific effects are equal to zero. An alternative approach is to introduce a random effect or frailty for each centre into the model. Recently, Commenges and Andersen have proposed a score test for this random effects model. By a Monte Carlo study we compare the performance of these two approaches when either the fixed or random effects model holds true. The study shows that for moderate samples the fixed effects tests have nominal levels much higher than specified, but the random effect test performs as expected under the null hypothesis. Under the alternative hypothesis the random effect test has good power to detect relatively small fixed or random centre effects. Also, if the centre effect is ignored the estimator of the main treatment effect may be quite biased and is inconsistent. The tests are illustrated on a retrospective multi-centre study of recovery from bone marrow transplantation.     

Factors associated with readmission to a general hospital in Brazil

The objective of this study was to compare different modeling strategies to identify individual and admissions characteristics associated with readmission to a general hospital. Routine data recorded in the Hospital Information System on all admissions to the Regional Public Hospital of Betim, Minas Gerais State, Brazil, from July 1996 to June 2000 were analyzed. Cox proportional hazards model and variants designed to deal with multiple-events data, like Andersen-Gill (AG), Prentice, Williams and Peterson (PWP), and random effects models were fitted to time between hospital admissions or censoring. For comparison purposes, a Poisson model was fitted to the total number of readmissions, using the same covariates. We analyzed 31,648 admissions of 26,198 patients, including 17,096 adults and 9,102 children. Estimates for the PWP and frailty models were very similar, and both approaches should be fitted and compared. If clinical characteristics are available, the PWP model should be used. Otherwise the random effects model can account for unmeasured differences, particularly some related to severity of the disease. These methodologies can help focus on various related readmission aspects such as diagnostic groups or medical specialties.     

Estimating and testing for center effects in competing risks

The problems of fitting Gaussian frailties proportional hazards models for the subdistribution of a competing risk and of testing for center effects are considered. In the analysis of competing risks data, Fine and Gray proposed a proportional hazards model for the subdistribution to directly assess the effects of covariates on the marginal failure probabilities of a given failure cause. Katsahianbiet al. extended their model to clustered time to event data, by including random center effects or frailties in the subdistribution hazard. We first introduce an alternate estimation procedure to the one proposed by Katsahian et al. This alternate estimation method is based on the penalized partial likelihood approach often used in fitting Gaussian frailty proportional hazards models in the standard survival analysis context, and has the advantage of using standard survival analysis software. Second, four hypothesis tests for the presence of center effects are given and compared via Monte-Carlo simulations. Statistical and numerical considerations lead us to formulate pragmatic guidelines as to which of the four tests is preferable. We also illustrate the proposed methodology with registry data from bone marrow transplantation for acute myeloid leukemia (AML).     

Frailty proportional mean residual life regression for clustered survival data: A hierarchical quasi-likelihood method

Frailty models are widely used to model clustered survival data arising in multicenter clinical studies. In the literature, most existing frailty models are proportional hazards, additive hazards, or accelerated failure time model based. In this paper, we propose a frailty model framework based on mean residual life regression to accommodate intracluster correlation and in the meantime provide easily understand and straightforward interpretation for the effects of prognostic factors on the expectation of the remaining lifetime. To overcome estimation challenges, a novel hierarchical quasi-likelihood approach is developed by making use of the idea of hierarchical likelihood in the construction of the quasi-likelihood function, leading to hierarchical estimating equations. Simulation results show favorable performance of the method regardless of frailty distributions. The utility of the proposed methodology is illustrated by its application to the data from a multi-institutional study of breast cancer.     

Recurrent injury event-time analysis

Public health decision making based on data sources that are characterized by a lack of independence and other complicating factors requires the development of innovative statistical techniques. Studies of injuries in occupational cohorts require methods to account for recurrent injuries to workers over time and the temporary removal of workers from the 'risk set' while recuperating. In this study, the times until injury events are modelled in an occupational cohort of employees in a large power utility company where employees are susceptible to recurrent events. The injury history over a ten-year period is used to compare the hazards of specific jobs, adjusted for age when first hired, and race/ethnicity differences. Subject-specific random effects and multiple event-times are accommodated through the application of frailty models which characterize the dependence of recurrent events over time. The counting process formulation of the proportional hazards regression model is used to estimate the effects of covariates for subjects with discontinuous intervals of risk. In this application, subjects are not at risk of injury during recovery periods or other illness, changes in jobs, or other reasons. Previous applications of proportional hazards regression in frailty models have not needed to account for the changing composition of the risk set which is required to adequately model occupational injury data. Published in 1999 by John Wiley & Sons, Ltd. This article is a US Government work and is in the public domain in the United States.     

Hierarchical likelihood inference on clustered competing risks data

The frailty model, an extension of the proportional hazards model, is often used to model clustered survival data. However, some extension of the ordinary frailty model is required when there exist competing risks within a cluster. Under competing risks, the underlying processes affecting the events of interest and competing events could be different but correlated. In this paper, the hierarchical likelihood method is proposed to infer the cause‐specific hazard frailty model for clustered competing risks data. The hierarchical likelihood incorporates fixed effects as well as random effects into an extended likelihood function, so that the method does not require intensive numerical methods to find the marginal distribution. Simulation studies are performed to assess the behavior of the estimators for the regression coefficients and the correlation structure among the bivariate frailty distribution for competing events. The proposed method is illustrated with a breast cancer dataset. Copyright © 2015 John Wiley & Sons, Ltd.     

Multilevel mixed effects parametric survival models using adaptive Gauss–Hermite quadrature with application to recurrent events and individual participant data meta‐analysis

Multilevel mixed effects survival models are used in the analysis of clustered survival data, such as repeated events, multicenter clinical trials, and individual participant data (IPD) meta‐analyses, to investigate heterogeneity in baseline risk and covariate effects. In this paper, we extend parametric frailty models including the exponential, Weibull and Gompertz proportional hazards (PH) models and the log logistic, log normal, and generalized gamma accelerated failure time models to allow any number of normally distributed random effects. Furthermore, we extend the flexible parametric survival model of Royston and Parmar, modeled on the log‐cumulative hazard scale using restricted cubic splines, to include random effects while also allowing for non‐PH (time‐dependent effects). Maximum likelihood is used to estimate the models utilizing adaptive or nonadaptive Gauss–Hermite quadrature. The methods are evaluated through simulation studies representing clinically plausible scenarios of a multicenter trial and IPD meta‐analysis, showing good performance of the estimation method. The flexible parametric mixed effects model is illustrated using a dataset of patients with kidney disease and repeated times to infection and an IPD meta‐analysis of prognostic factor studies in patients with breast cancer. User‐friendly Stata software is provided to implement the methods. Copyright © 2014 John Wiley & Sons, Ltd.     

Bayesian transformation cure frailty models with multivariate failure time data

We propose a class of transformation cure frailty models to accommodate a survival fraction in multivariate failure time data. Established through a general power transformation, this family of cure frailty models includes the proportional hazards and the proportional odds modeling structures as two special cases. Within the Bayesian paradigm, we obtain the joint posterior distribution and the corresponding full conditional distributions of the model parameters for the implementation of Gibbs sampling. Model selection is based on the conditional predictive ordinate statistic and deviance information criterion. As an illustration, we apply the proposed method to a real data set from dentistry.     

The use of Gaussian quadrature for estimation in frailty proportional hazards models

In this paper, we propose a novel Gaussian quadrature estimation method in various frailty proportional hazards models. We approximate the unspecified baseline hazard by a piecewise constant one, resulting in a parametric model that can be fitted conveniently by Gaussian quadrature tools in standard software such as SAS Proc NLMIXED. We first apply our method to simple frailty models for correlated survival data (e.g. recurrent or clustered failure times), then to joint frailty models for correlated failure times with informative dropout or a dependent terminal event such as death. Simulation studies show that our method compares favorably with the well-received penalized partial likelihood method and the Monte Carlo EM (MCEM) method, for both normal and Gamma frailty models. We apply our method to three real data examples: (1) the time to blindness of both eyes in a diabetic retinopathy study, (2) the joint analysis of recurrent opportunistic diseases in the presence of death for HIV-infected patients, and (3) the joint modeling of local, distant tumor recurrences and patients survival in a soft tissue sarcoma study. The proposed method greatly simplifies the implementation of the (joint) frailty models and makes them much more accessible to general statistical practitioners.     

THE ROLE OF FRAILTY MODELS AND ACCELERATED FAILURE TIME MODELS IN DESCRIBING HETEROGENEITY DUE TO OMITTED COVARIATES

In survival analysis, deviations from proportional hazards may sometimes be explained by unaccounted random heterogeneity, or frailty. This paper recalls the literature on omitted covariates in survival analysis and shows in a case study how unstably frailty models might behave when asked to account for unobserved heterogeneity in standard survival analysis with no replications per heterogeneity unit. Accelerated failure time modelling seems to avoid these difficulties and also to yield easily interpretable results. We propose that it would be advantageous to upgrade the accelerated failure time approach alongside the hazard modelling approach to survival analysis. © 1997 by John Wiley & Sons, Ltd.     

Joint modeling of survival time and longitudinal data with subject-specific changepoints in the covariates

Joint models are frequently used in survival analysis to assess the relationship between time-to-event data and time-dependent covariates, which are measured longitudinally but often with errors. Routinely, a linear mixed-effects model is used to describe the longitudinal data process, while the survival times are assumed to follow the proportional hazards model. However, in some practical situations, individual covariate profiles may contain changepoints. In this article, we assume a two-phase polynomial random effects with subject-specific changepoint model for the longitudinal data process and the proportional hazards model for the survival times. Our main interest is in the estimation of the parameter in the hazards model. We incorporate a smooth transition function into the changepoint model for the longitudinal data and develop the corrected score and conditional score estimators, which do not require any assumption regarding the underlying distribution of the random effects or that of the changepoints. The estimators are shown to be asymptotically equivalent and their finite-sample performance is examined via simulations. The methods are applied to AIDS clinical trial data.