Bias reduction for semi-competing risks frailty model with rare events: application to a chronic kidney disease cohort study in South Korea

In a semi-competing risks model in which a terminal event censors a non-terminal event but not vice versa, the conventional method can predict clinical outcomes by maximizing likelihood estimation. However, this method can produce unreliable or biased estimators when the number of events in the datasets is small. Specifically, parameter estimates may converge to infinity, or their standard errors can be very large. Moreover, terminal and non-terminal event times may be correlated, which can account for the frailty term. Here, we adapt the penalized likelihood with Firth's correction method for gamma frailty models with semi-competing risks data to reduce the bias caused by rare events. The proposed method is evaluated in terms of relative bias, mean squared error, standard error, and standard deviation compared to the conventional methods through simulation studies. The results of the proposed method are stable and robust even when data contain only a few events with the misspecification of the baseline hazard function. We also illustrate a real example with a multi-centre, patient-based cohort study to identify risk factors for chronic kidney disease progression or adverse clinical outcomes. This study will provide a better understanding of semi-competing risk data in which the number of specific diseases or events of interest is rare.

Penalized variable selection in multi-parameter regression survival modeling

Standard survival models such as the proportional hazards model contain a single regression component, corresponding to the scale of the hazard. In contrast, we consider the so-called "multi-parameter regression" approach whereby covariates enter the model through multiple distributional parameters simultaneously, for example, scale and shape parameters. This approach has previously been shown to achieve flexibility with relatively low model complexity. However, beyond a stepwise type selection method, variable selection methods are underdeveloped in the multi-parameter regression survival modeling setting. Therefore, we propose penalized multi-parameter regression estimation procedures using the following penalties: least absolute shrinkage and selection operator, smoothly clipped absolute deviation, and adaptive least absolute shrinkage and selection operator. We compare these procedures using extensive simulation studies and an application to data from an observational lung cancer study; the Weibull multi-parameter regression model is used throughout as a running example.

Penalized variable selection for cause-specific hazard frailty models with clustered competing-risks data

Competing risks data usually arise when an occurrence of an event precludes other types of events from being observed. Such data are often encountered in a clustered clinical study such as a multi-center clinical trial. For the clustered competing-risks data which are correlated within a cluster, competing-risks models allowing for frailty terms have been recently studied. To the best of our knowledge, however, there is no literature on variable selection methods for cause-specific hazard frailty models. In this article, we propose a variable selection procedure for fixed effects in cause-specific competing risks frailty models using a penalized h-likelihood (HL). Here, we study three penalty functions, LASSO, SCAD, and HL. Simulation studies demonstrate that the proposed procedure using the HL penalty works well, providing a higher probability of choosing the true model than LASSO and SCAD methods without losing prediction accuracy. The proposed method is illustrated by using two kinds of clustered competing-risks cancer data sets.

Semiparametric estimation for nonparametric frailty models using nonparametric maximum likelihood approach

A consequence of using a parametric frailty model with nonparametric baseline hazard for analyzing clustered time-to-event data is that its regression coefficient estimates could be sensitive to the underlying frailty distribution. Recently, there has been a proposal for specifying both the baseline hazard and the frailty distribution nonparametrically, and estimating the unknown parameters by the maximum penalized likelihood method. Instead, in this paper, we propose the nonparametric maximum likelihood method for a general class of nonparametric frailty models, i.e. models where the frailty distribution is completely unspecified but the baseline hazard can be either parametric or nonparametric. The implementation of the estimation procedure can be based on a combination of either the Broyden-Fletcher-Goldfarb-Shanno or expectation-maximization algorithm and the constrained Newton algorithm with multiple support point inclusion. Simulation studies to investigate the performance of estimation of a regression coefficient by several different model-fitting methods were conducted. The simulation results show that our proposed regression coefficient estimator generally gives a reasonable bias reduction when the number of clusters is increased under various frailty distributions. Our proposed method is also illustrated with two data examples.

Flexible parametric copula modeling approaches for clustered survival data

Copula-based survival regression models, which consist of a copula function and marginal distribution (i.e., marginal survival function), have been widely used for analyzing clustered multivariate survival data. Archimedean copula functions are useful for modeling such dependence. For the likelihood inference, one-stage and two-stage estimation methods have been usually used. The two-stage procedure can give inefficient estimation results because of separate estimation of the marginal and copula's dependence parameters. The more efficient one-stage procedure has been mainly developed under a restrictive parametric assumption of marginal distribution due to complexity of the full likelihood with unknown marginal baseline hazard functions. In this paper, we propose a flexible parametric Archimedean copula modeling approach using a one-stage likelihood procedure. In order to reduce the complexity of the full likelihood, the unknown marginal baseline hazards are modeled based on a cubic M-spline basis function that does not require a specific parametric form. Simulation results demonstrate that the proposed one-stage estimation method gives a consistent estimator and also provides more efficient results over existing one- and two-stage methods. The new method is illustrated with three clinical data sets. The Appendix provides an R function so that the proposed method becomes directly accessible to interested readers.

Frailty modelling approaches for semi-competing risks data

In the semi-competing risks situation where only a terminal event censors a non-terminal event, observed event times can be correlated. Recently, frailty models with an arbitrary baseline hazard have been studied for the analysis of such semi-competing risks data. However, their maximum likelihood estimator can be substantially biased in the finite samples. In this paper, we propose effective modifications to reduce such bias using the hierarchical likelihood. We also investigate the relationship between marginal and hierarchical likelihood approaches. Simulation results are provided to validate performance of the proposed method. The proposed method is illustrated through analysis of semi-competing risks data from a breast cancer study.

Comparison of the marginal hazard model and the sub-distribution hazard model for competing risks under an assumed copula

For the analysis of competing risks data, three different types of hazard functions have been considered in the literature, namely the cause-specific hazard, the sub-distribution hazard, and the marginal hazard function. Accordingly, medical researchers can fit three different types of the Cox model to estimate the effect of covariates on each of the hazard function. While the relationship between the cause-specific hazard and the sub-distribution hazard has been extensively studied, the relationship to the marginal hazard function has not yet been analyzed due to the difficulties related to non-identifiability. In this paper, we adopt an assumed copula model to deal with the model identifiability issue, making it possible to establish a relationship between the sub-distribution hazard and the marginal hazard function. We then compare the two methods of fitting the Cox model to competing risks data. We also extend our comparative analysis to clustered competing risks data that are frequently used in medical studies. To facilitate the numerical comparison, we implement the computing algorithm for marginal Cox regression with clustered competing risks data in the R joint.Cox package and check its performance via simulations. For illustration, we analyze two survival datasets from lung cancer and bladder cancer patients.

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.

Penalized variable selection for accelerated failure time models with random effects

Accelerated failure time (AFT) models allowing for random effects are linear mixed models under the log-transformation of survival time with censoring and describe dependence in correlated survival data. It is well known that the AFT models are useful alternatives to frailty models. To the best of our knowledge, however, there is no literature on variable selection methods for such AFT models. In this paper, we propose a simple but unified variable-selection procedure of fixed effects in the AFT random-effect models using penalized h-likelihood (HL). We consider four penalty functions (ie, least absolute shrinkage and selection operator (LASSO), adaptive LASSO, smoothly clipped absolute deviation (SCAD), and HL). We show that the proposed method can be easily implemented via a slight modification to existing h-likelihood estimation procedures. We thus demonstrate that the proposed method can also be easily extended to AFT models with multilevel (or nested) structures. Simulation studies also show that the procedure using the adaptive LASSO, SCAD, or HL penalty performs well. In particular, we find via the simulation results that the variable selection method with HL penalty provides a higher probability of choosing the true model than other three methods. The usefulness of the new method is illustrated using two actual datasets from multicenter clinical trials.

A Methodological Perspective on the Longitudinal Cognitive Change after Stroke

BACKGROUND/AIMS

Most studies of poststroke cognitive impairment (PSCI) have analyzed cognitive levels at specific time points rather than their changes over time. Furthermore, they seldom consider correlations between cognitive domains. We aimed to investigate the effects of these methodological considerations on determining significant PSCI predictors in a longitudinal stroke cohort.

METHODS

In patients who underwent neuropsychological tests at least twice after stroke, we adopted a multilevel hierarchical mixed-effects model with domain-specific cognitive changes and a multivariate model for multiple outcomes to reflect their correlations.

RESULTS

We enrolled 375 patients (median follow-up of 34.1 months). Known predictors of PSCI were generally associated with cognitive levels; however, most of the statistical significances disappeared when cognitive changes were set as outcomes, except age for memory, prior stroke and baseline cognition for executive/attention domain, and baseline cognition for visuospatial function. The multivariate analysis which considered multiple outcomes simultaneously further altered these associations.

CONCLUSIONS

This study shows that defining outcomes as changes over time and reflecting correlations between outcomes may affect the identification of predictors of PSCI.

H-likelihood approach for joint modeling of longitudinal outcomes and time-to-event data

In longitudinal studies, a subject may have different types of outcomes that could be correlated. For example, a response variable of interest would be measured repeatedly over time on the same subject and at the same time, an event time representing a single event or competing-risks event is also observed. In this paper, we propose a joint modeling framework that accounts for the inherent association between such multiple outcomes via frailties (unobserved random effects). Among outcomes, at least one outcome is an event time that has a type of a single event or competing-risks event. For inference we use the hierarchical likelihood (h-likelihood) that provides an unified efficient fitting procedure for the joint models. Numerical studies are provided to show the performance of the proposed method and two data examples are shown.

Analysis of clustered competing risks data using subdistribution hazard models with multivariate frailties

Competing risks data often exist within a center in multi-center randomized clinical trials where the treatment effects or baseline risks may vary among centers. In this paper, we propose a subdistribution hazard regression model with multivariate frailty to investigate heterogeneity in treatment effects among centers from multi-center clinical trials. For inference, we develop a hierarchical likelihood (or h-likelihood) method, which obviates the need for an intractable integration over the frailty terms. We show that the profile likelihood function derived from the h-likelihood is identical to the partial likelihood, and hence it can be extended to the weighted partial likelihood for the subdistribution hazard frailty models. The proposed method is illustrated with a dataset from a multi-center clinical trial on breast cancer as well as with a simulation study. We also demonstrate how to present heterogeneity in treatment effects among centers by using a confidence interval for the frailty for each individual center and how to perform a statistical test for such heterogeneity using a restricted h-likelihood.

Frailty modeling for clustered competing risks data with missing cause of failure

Competing risks data often occur within a center in multi-center clinical trials where the event times within a center may be correlated due to unobserved factors across individuals. In this paper, we consider the cause-specific proportional hazards model with a shared frailty to model the association between the event times within a center in the framework of competing risks. We use a hierarchical likelihood approach, which does not require any intractable integration over the frailty terms. In a clinical trial, cause of death information may not be observed for some patients. In such a case, analyses through exclusion of cases with missing cause of death may lead to biased inferences. We propose a hierarchical likelihood approach for fitting the cause-specific proportional hazards model with a shared frailty in the presence of missing cause of failure. We use multiple imputation methods to address missing cause of death information under the assumption of missing at random. Simulation studies show that the proposed procedures perform well, even if the imputation model is misspecified. The proposed methods are illustrated with data from EORTC trial 30791 conducted by European Organization for Research and Treatment of Cancer (EORTC).

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.

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.

Interval estimation of random effects in proportional hazards models with frailties

Semi-parametric frailty models are widely used to analyze clustered survival data. In this article, we propose the use of the hierarchical likelihood interval for individual frailties. We study the relationship between hierarchical likelihood, empirical Bayesian, and fully Bayesian intervals for frailties. We show that our proposed interval can be interpreted as a frequentist confidence interval and Bayesian credible interval under a uniform prior. We also propose an adjustment of the proposed interval to avoid null intervals. Simulation studies show that the proposed interval preserves the nominal confidence level. The procedure is illustrated using data from a multicenter lung cancer clinical trial.

Frailty modelling for survival data from multi-centre clinical trials

Despite the use of standardized protocols in, multi-centre, randomized clinical trials, outcome may vary between centres. Such heterogeneity may alter the interpretation and reporting of the treatment effect. Below, we propose a general frailty modelling approach for investigating, inter alia, putative treatment-by-centre interactions in time-to-event data in multi-centre clinical trials. A correlated random effects model is used to model the baseline risk and the treatment effect across centres. It may be based on shared, individual or correlated random effects. For inference we develop the hierarchical-likelihood (or h-likelihood) approach which facilitates computation of prediction intervals for the random effects with proper precision. We illustrate our methods using disease-free time-to-event data on bladder cancer patients participating in an European Organization for Research and Treatment of Cancer trial, and a simulation study. We also demonstrate model selection using h-likelihood criteria.

Robust frailty modelling using non-proportional hazards models

Correlated survival times can be modelled by introducing a random effect, or frailty component, into the hazard function. For multivariate survival data, we extend a non-proportional hazards (PH) model, the generalized time-dependent logistic survival model, to include random effects. The hierarchical likelihood procedure, which obviates the need for marginalization over the random effect distribution, is derived for this extended model and its properties are discussed. The extended model leads to a robust estimation result for the regression parameters against the misspecification of the form of the basic hazard function or frailty distribution compared to PH-based alternatives. The proposed method is illustrated by two practical examples and a simulation study which demonstrate the advantages of the new model.

Model selection for multi-component frailty models

Various frailty models have been developed and are now widely used for analysing multivariate survival data. It is therefore important to develop an information criterion for model selection. However, in frailty models there are several alternative ways of forming a criterion and the particular criterion chosen may not be uniformly best. In this paper, we study an Akaike information criterion (AIC) on selecting a frailty structure from a set of (possibly) non-nested frailty models. We propose two new AIC criteria, based on a conditional likelihood and an extended restricted likelihood (ERL) given by Lee and Nelder (J. R. Statist. Soc. B 1996; 58:619-678). We compare their performance using well-known practical examples and demonstrate that the two criteria may yield rather different results. A simulation study shows that the AIC based on the ERL is recommended, when attention is focussed on selecting the frailty structure rather than the fixed effects.

Genetic mixed linear models for twin survival data

Twin studies are useful for assessing the relative importance of genetic or heritable component from the environmental component. In this paper we develop a methodology to study the heritability of age-at-onset or lifespan traits, with application to analysis of twin survival data. Due to limited period of observation, the data can be left truncated and right censored (LTRC). Under the LTRC setting we propose a genetic mixed linear model, which allows general fixed predictors and random components to capture genetic and environmental effects. Inferences are based upon the hierarchical-likelihood (h-likelihood), which provides a statistically efficient and unified framework for various mixed-effect models. We also propose a simple and fast computation method for dealing with large data sets. The method is illustrated by the survival data from the Swedish Twin Registry. Finally, a simulation study is carried out to evaluate its performance.

Dispersion frailty models and HGLMs

In medical research recurrent event times can be analysed using a frailty model in which the frailties for different individuals are independent and identically distributed. However, such a homogeneous assumption about frailties could sometimes be suspect. For modelling heterogeneity in frailties we describe dispersion frailty models arising from a new class of models, namely hierarchical generalized linear models. Using the kidney infection data we illustrate how to detect and model heterogeneity among frailties. Stratification of frailty models is also investigated.

Comparison of the Berman Intubating Airway and the Williams Airway Intubator for fibreoptic orotracheal intubation in anaesthetised patients

Sixty patients with no clinical indicators of a difficult airway were selected to undergo a fibreoptic assessment after induction of general anaesthesia using both the Berman Intubating Airway and the Williams Airway Intubator. The bronchoscopic view and ease of railroading a tracheal tube during fibreoptic orotracheal intubation were studied. The bronchoscopic views obtained by the Berman Airway and the Williams Airway were significantly different (p < 0.008, test of symmetry). The estimated odds ratio of obtaining an obstructed path was 3.06 times higher for the Berman than the Williams Airway. However, if the glottis could be reached with the bronchoscope, there was no significant difference in the degree of ease of intubation between the two airways.

Multilevel mixed linear models for survival data

For the analysis of correlated survival data mixed linear models are useful alternatives to frailty models. By their use the survival times can be directly modelled, so that the interpretation of the fixed and random effects is straightforward. However, because of intractable integration involved with the use of marginal likelihood the class of models in use has been severely restricted. Such a difficulty can be avoided by using hierarchical-likelihood, which provides a statistically efficient and fast fitting algorithm for multilevel models. The proposed method is illustrated using the chronic granulomatous disease data. A simulation study is carried out to evaluate the performance.

Hierarchical-likelihood approach for mixed linear models with censored data

Mixed linear models describe the dependence via random effects in multivariate normal survival data. Recently they have received considerable attention in the biomedical literature. They model the conditional survival times, whereas the alternative frailty model uses the conditional hazard rate. We develop an inferential method for the mixed linear model via Lee and Nelder's (1996) hierarchical-likelihood (h-likelihood). Simulation and a practical example are presented to illustrate the new method.