Model selection and diagnostics for joint modeling of survival and longitudinal data with crossing hazard rate functions

Comparison of two hazard rate functions is important for evaluating treatment effect in studies concerning times to some important events. In practice, it may happen that the two hazard rate functions cross each other at one or more unknown time points, representing temporal changes of the treatment effect. Also, besides survival data, there could be longitudinal data available regarding some time‐dependent covariates. When jointly modeling the survival and longitudinal data in such cases, model selection and model diagnostics are especially important to provide reliable statistical analysis of the data, which are lacking in the literature. In this paper, we discuss several criteria for assessing model fit that have been used for model selection and apply them to the joint modeling of survival and longitudinal data for comparing two crossing hazard rate functions. We also propose hypothesis testing and graphical methods for model diagnostics of the proposed joint modeling approach. Our proposed methods are illustrated by a simulation study and by a real‐data example concerning two early breast cancer treatments. Copyright © 2014 John Wiley & Sons, Ltd.     

An evaluation of resampling methods for assessment of survival risk prediction in high-dimensional settings

Resampling techniques are often used to provide an initial assessment of accuracy for prognostic prediction models developed using high-dimensional genomic data with binary outcomes. Risk prediction is most important, however, in medical applications and frequently the outcome measure is a right-censored time-to-event variable such as survival. Although several methods have been developed for survival risk prediction with high-dimensional genomic data, there has been little evaluation of the use of resampling techniques for the assessment of such models. Using real and simulated datasets, we compared several resampling techniques for their ability to estimate the accuracy of risk prediction models. Our study showed that accuracy estimates for popular resampling methods, such as sample splitting and leave-one-out cross validation (Loo CV), have a higher mean square error than for other methods. Moreover, the large variability of the split-sample and Loo CV may make the point estimates of accuracy obtained using these methods unreliable and hence should be interpreted carefully. A k-fold cross-validation with k = 5 or 10 was seen to provide a good balance between bias and variability for a wide range of data settings and should be more widely adopted in practice.     

A new measure of prognostic separation in survival data

Multivariable prognostic models are widely used in cancer and other disease areas, and have a range of applications in clinical medicine, clinical trials and allocation of health services resources. A well-founded and reliable measure of the prognostic ability of a model would be valuable to help define the separation between patients or prognostic groups that the model could provide, and to act as a benchmark of model performance in a validation setting. We propose such a measure for models of survival data. Its motivation derives originally from the idea of separation between Kaplan-Meier curves. We define the criteria for a successful measure and discuss them with respect to our approach. Adjustments for 'optimism', the tendency for a model to predict better on the data on which it was derived than on new data, are suggested. We study the properties of the measure by simulation and by example in three substantial data sets. We believe that our new measure will prove useful as a tool to evaluate the separation available-with a prognostic model.     

A non-parametric procedure for evaluating treatment effect in the meta-analysis of survival data

This paper addresses the problem of combining information from independent clinical trials which compare survival distributions of two treatment groups. Current meta-analytic methods which take censoring into account are often not feasible for meta-analyses which synthesize summarized results in published (or unpublished) references, as these methods require information usually not reported. The paper presents methodology which uses the log(-log) survival function difference, (i.e. log(-logS2(t))-log(-logS1(t)), as the contrast index to represent the multiplicative treatment effect on survival in independent trials. This article shows by the second mean value theorem for integrals that this contrast index, denoted as theta, is interpretable as a weighted average on a natural logarithmic scale of hazard ratios within the interval [0,t] in a trial. When the within-trial proportional hazards assumption is true, theta is the logarithm of the proportionality constant for the common hazard ratio for the interval considered within the trial. In this situation, an important advantage of using theta as a contrast index in the proposed methodology is that the estimation of theta is not affected by length of follow-up time. Other commonly used indices such as the odds ratio, risk ratio and risk differences do not have this invariance property under the proportional hazard model, since their estimation may be affected by length of follow-up time as a technical artefact. Thus, the proposed methodology obviates problems which often occur in survival meta-analysis because trials do not report survival at the same length of follow-up time. Even when the within-trial proportional hazards assumption is not realistic, the proposed methodology has the capability of testing a global null hypothesis of no multiplicative treatment effect on the survival distributions of two groups for all studies. A discussion of weighting schemes for meta-analysis is provided, in particular, a weighting scheme based on effective sample sizes is suggested for the meta-analysis of time-to-event data which involves censoring. A medical example illustrating the methodology is given. A simulation investigation suggested that the methodology performs well in the presence of moderate censoring.     

A testing procedure for survival data with few responders.

In the course of designing a clinical trial, investigators are often faced with the possibility that only a fraction of the patients will benefit from the experimental treatment. A proper clinical trial design requires prospective specification of the testing procedures to be used in the analysis. In the absence of reliable prognostic factors capable of identifying the appropriate subset of patients, there is a need for a test procedure that will be sensitive to a range of possible fractions of responders. Focusing on survival data, we propose guidelines for selecting a proper test procedure based on the anticipated proportion of responding patients. These guidelines suggest that the logrank test should be used when the fraction of responders is expected to be greater than 0.5, otherwise procedures based on weighted linear rank tests are preferable. Overall this approach provides good power properties when the treatment affects only a small proportion of patients while protecting against substantial loss of power when all patients are affected. Use of the procedure is illustrated with data from two published randomized studies.     

Assessing potentially time‐dependent treatment effect from clinical trials and observational studies for survival data,with applications to the Women's Health Initiative combined hormone therapy trial

For risk and benefit assessment in clinical trials and observational studies with time‐to‐event data, the Cox model has usually been the model of choice. When the hazards are possibly non‐proportional, a piece‐wise Cox model over a partition of the time axis may be considered. Here, we propose to analyze clinical trials or observational studies with time‐to‐event data using a certain semiparametric model. The model allows for a time‐dependent treatment effect. It includes the important proportional hazards model as a sub‐model and can accommodate various patterns of time‐dependence of the hazard ratio. After estimation of the model parameters using a pseudo‐likelihood approach, simultaneous confidence intervals for the hazard ratio function are established using a Monte Carlo method to assess the time‐varying pattern of the treatment effect. To assess the overall treatment effect, estimated average hazard ratio and its confidence intervals are also obtained. The proposed methods are applied to data from the Women's Health Initiative. To compare the Women's Health Initiative clinical trial and observational study, we use the propensity score in building the regression model. Compared with the piece‐wise Cox model, the proposed model yields a better model fit and does not require partitioning of the time axis. Copyright © 2015 John Wiley & Sons, Ltd.     

Doubly robust inference procedure for relative survival ratio in population-based cancer registry data

Cancer registry system has been playing important roles in research and policy making in cancer control. In general, information on cause of death is not available in cancer registry data. To make inference on survival of cancer patients in the absence of cause of death information, the relative survival ratio is widely used in the population-based cancer research utilizing external life tables for the general population. Another difficulty arising in analyzing cancer registry data is informative censoring. In this article, we propose a doubly robust inference procedure for the relative survival ratio under a certain type of informative censoring, called the covariate-dependent censoring. The proposed estimator is doubly robust in the sense that it is consistent if at least one of the regression models for the time-to-death and for the censoring time is correctly specified. Furthermore, we introduced a doubly robust test assessing underlying conditional independence assumption between the time-to-death and the censoring time. This test is model based, but is doubly robust in the sense that at least one of the models for the time to event and for the censoring time is correctly specified, it maintains its nominal significance level. This notable feature entails us to make inference on cancer registry data relying on assumptions, which are much weaker than the existing methods and are verifiable empirically. We examine the theoretical and empirical properties of our proposed methods by asymptotic theory and simulation studies. We illustrate the proposed method with cancer registry data in Osaka, Japan.     

Case-cohort methods for survival data on families from routine registers

In the Nordic countries, there exist several registers containing information on diseases and risk factors for millions of individuals. This information can be linked to families by the use of personal identification numbers, and represents a great opportunity for studying diseases that show familial aggregation. Due to the size of the registers, it is difficult to analyze the data by using traditional methods for multivariate survival analysis, such as frailty or copula models. Since the size of the cohort is known, case-cohort methods based on pseudo-likelihoods are suitable for analyzing the data. We present methods for sampling control families both with and without replacement, and with or without stratification. The data are stratified according to family size and covariate values. Depending on the sampling method, results from simulations indicate that one only needs to sample 1-5 per cent of the control families in order to obtain good efficiency compared with a traditional cohort analysis. We also provide an application to survival data from the Medical Birth Registry of Norway.     

A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non‐linear and non‐proportional effects of covariates

The excess hazard regression model is an approach developed for the analysis of cancer registry data to estimate net survival, that is, the survival of cancer patients that would be observed if cancer was the only cause of death. Cancer registry data typically possess a hierarchical structure: individuals from the same geographical unit share common characteristics such as proximity to a large hospital that may influence access to and quality of health care, so that their survival times might be correlated. As a consequence, correct statistical inference regarding the estimation of net survival and the effect of covariates should take this hierarchical structure into account. It becomes particularly important as many studies in cancer epidemiology aim at studying the effect on the excess mortality hazard of variables, such as deprivation indexes, often available only at the ecological level rather than at the individual level. We developed here an approach to fit a flexible excess hazard model including a random effect to describe the unobserved heterogeneity existing between different clusters of individuals, and with the possibility to estimate non‐linear and time‐dependent effects of covariates. We demonstrated the overall good performance of the proposed approach in a simulation study that assessed the impact on parameter estimates of the number of clusters, their size and their level of unbalance. We then used this multilevel model to describe the effect of a deprivation index defined at the geographical level on the excess mortality hazard of patients diagnosed with cancer of the oral cavity. Copyright © 2016 John Wiley & Sons, Ltd.     

Accuracy of predictive ability measures for survival models

One aspect of an analysis of survival data based on the proportional hazards model that has been receiving increasing attention is that of the predictive ability or explained variation of the model. A number of contending measures have been suggested, including one measure, R²(β), which has been proposed given its several desirable properties, including its capacity to accommodate time‐dependent covariates, a major feature of the model and one that gives rise to great generality. A thorough study of the properties of available measures, including the aforementioned measure, has been carried out recently. In that work, the authors used bootstrap techniques, particularly complex in the setting of censored data, in order to obtain estimates of precision. The motivation of this work is to provide analytical expressions of precision, in particular confidence interval estimates for R²(β). We use Taylor series approximations with and without local linearizing transforms. We also consider a very simple expression based on the Fisher's transformation. This latter approach has two great advantages. It is very easy and quick to calculate, and secondly, it can be obtained for any of the methods given in the recent review. A large simulation study is carried out to investigate the properties of the different methods. Finally, three well‐known datasets in breast cancer, lymphoma and lung cancer research are given as illustrations. Copyright © 2017 John Wiley & Sons, Ltd.     

Application of functional ANOVA models for hazard regression to the Hisayama data

A methodology for modeling covariate effects on the time-to-event data is developed. The covariates are allowed to be time dependent and their effects are modeled using polynomial splines in order to account for possibly non-linear effects. The methodology is applied to examine the effects on the incidence brain infarction based on a cohort study in Hisayama, Japan. The results indicate that at least two non-linear effects are significant (body mass index and systolic blood pressure) and there is a time-varying drug effect. The resulting significant risk factors are assessed by the proposed method that is more flexible and hence less biased than the traditional procedures where linear effects are imposed. These results are extremely important to the local medical investigation. In particular, more insight has been gained by examining the non-linear effects.     

A stochastic multiple imputation algorithm for missing covariate data in tree-structured survival analysis

Missing covariate data present a challenge to tree-structured methodology due to the fact that a single tree model, as opposed to an estimated parameter value, may be desired for use in a clinical setting. To address this problem, we suggest a multiple imputation algorithm that adds draws of stochastic error to a tree-based single imputation method presented by Conversano and Siciliano (Technical Report, University of Naples, 2003). Unlike previously proposed techniques for accommodating missing covariate data in tree-structured analyses, our methodology allows the modeling of complex and nonlinear covariate structures while still resulting in a single tree model. We perform a simulation study to evaluate our stochastic multiple imputation algorithm when covariate data are missing at random and compare it to other currently used methods. Our algorithm is advantageous for identifying the true underlying covariate structure when complex data and larger percentages of missing covariate observations are present. It is competitive with other current methods with respect to prediction accuracy. To illustrate our algorithm, we create a tree-structured survival model for predicting time to treatment response in older, depressed adults.     

Comparison of methods for survival analysis of dependent data

Analysis of dependent survival data by conventional partial likelihood methods produces unbiased estimates of the regression coefficients but incorrectly estimates their variance. Here we compared the conventional partial likelihood methods with two alternative methods for analyzing dependent survival data. The first alternative method estimated the regression coefficient by the partial likelihood approach but adjusted the variance to account for clustering. The second alternative method used marginal likelihoods to estimate both the regression coefficient and its variance. We evaluated the performance of the three methods using simulated and actual data. Simulated data were used to examine bias, efficiency, type I errors, and power. An Old Order Amish genealogy was analyzed under these models to illustrate their performance on real data. The simulation study showed that all three methods provided unbiased estimates of the regression coefficient, but the efficiency of the estimated regression coefficient varied according to the simulation conditions. The standard partial likelihood method showed increasing type I error as the dependence increased within clusters. Both alternative methods had acceptable levels of type I errors at all dependence levels. In the analysis of genealogic data, the regression coefficient was similar in the three methods showing stable estimates of the regression coefficients. The variance estimates from the alternative methods were slightly different from the conventional method, suggesting a flow level of dependence. This study displays the effect of violating the independence assumption and provides guidelines for using alternative statistical methods. © 1996 Wiley-Liss, Inc.     

Interval estimation of risk difference for data sampled from clusters

Risk difference (RD) is an important measure in epidemiological studies where the probability of developing a disease for individuals in an exposed group, for example, is compared with that in a control group. There are varying cluster sizes in each group and the binary responses within each cluster cannot be assumed independent. Under the cluster sampling scenario, Lui (Statistical Estimation of Epidemiological Risk. Wiley: CA, 2004; 7-27) discusses four methods for the construction of a confidence interval for the RD. In this paper we introduce two very simple methods. One method is based on an estimator of the variance of a ratio estimator (Sampling Techniques (3rd edn). Wiley: New York, 1977; 30-67) and the other method is based on a sandwich estimator of the variance of the regression estimator using the generalized estimating equations approach of Zeger and Liang (Biometrics 1986; 42:121-130). These two methods are then compared, by simulation, in terms of maintaining nominal coverage probability and average coverage length, with the four methods discussed by Lui (Statistical Estimation of Epidemiological Risk. Wiley: CA, 2004; 7-27). Simulations show at least as good properties of these two methods as those of the others. The method based on an estimate of the variance of a ratio estimator performs best overall. It involves a very simple variance expression and can be implemented with a very few computer codes. Therefore, it can be considered as an easily implementable alternative.     

A covariate-specific time-dependent receiver operating characteristic curve for correlated survival data

Several studies for the clinical validity of circulating tumor cells (CTCs) in metastatic breast cancer were conducted showing that it is a prognostic biomarker of overall survival. In this work, we consider an individual patient data meta-analysis for nonmetastatic breast cancer to assess the discrimination of CTCs regarding the risk of death. Data are collected in several centers and present correlated failure times for subjects of the same center. However, although the covariate-specific time-dependent receiver operating characteristic (ROC) curve has been widely used for assessing the performance of a biomarker, there is no methodology yet that can handle this specific setting with clustered censored failure times. We propose an estimator for the covariate-specific time-dependent ROC curves and area under the ROC curve when clustered failure times are detected. We discuss the assumptions under which the estimators are consistent and their interpretations. We assume a shared frailty model for modeling the effect of the covariates and the biomarker on the outcome in order to account for the cluster effect. A simulation study was conducted and it shows negligible bias for the proposed estimator and a nonparametric one based on inverse probability censoring weighting, while a semiparametric estimator, ignoring the clustering, is markedly biased. Finally, in our application to breast cancer data, the estimation of the covariate-specific area under the curves illustrates that the CTCs discriminate better patients with inflammatory tumor than patients with noninflammatory tumor, with respect to their risk of death.     

A review of the use of time-varying covariates in the Fine-Gray subdistribution hazard competing risk regression model

In survival analysis, time-varying covariates are covariates whose value can change during follow-up. Outcomes in medical research are frequently subject to competing risks (events precluding the occurrence of the primary outcome). We review the types of time-varying covariates and highlight the effect of their inclusion in the subdistribution hazard model. External time-dependent covariates are external to the subject, can effect the failure process, but are not otherwise involved in the failure mechanism. Internal time-varying covariates are measured on the subject, can effect the failure process directly, and may also be impacted by the failure mechanism. In the absence of competing risks, a consequence of including internal time-dependent covariates in the Cox model is that one cannot estimate the survival function or the effect of covariates on the survival function. In the presence of competing risks, the inclusion of internal time-varying covariates in a subdistribution hazard model results in the loss of the ability to estimate the cumulative incidence function (CIF) or the effect of covariates on the CIF. Furthermore, the definition of the risk set for the subdistribution hazard function can make defining internal time-varying covariates difficult or impossible. We conducted a review of the use of time-varying covariates in subdistribution hazard models in articles published in the medical literature in 2015 and in the first 5 months of 2019. Seven percent of articles published included a time-varying covariate. Several inappropriately described a time-varying covariate as having an association with the risk of the outcome.     

A method for determining groups in multiple survival curves

Survival analysis includes a wide variety of methods for analyzing time-to-event data. One basic but important goal in survival analysis is the comparison of survival curves between groups. Several nonparametric methods have been proposed in the literature to test for the equality of survival curves for censored data. When the null hypothesis of equality of curves is rejected, leading to the clear conclusion that at least one curve is different, it can be interesting to ascertain whether curves can be grouped or if all these curves are different from each other. A method is proposed that allows determining groups with an automatic selection of their number. The validity and behavior of the proposed method was evaluated through simulation studies. The applicability of the proposed method is illustrated using real data. Software in the form of an R package has been developed implementing the proposed method.     

The median hazard ratio: a useful measure of variance and general contextual effects in multilevel survival analysis

Multilevel data occurs frequently in many research areas like health services research and epidemiology. A suitable way to analyze such data is through the use of multilevel regression models (MLRM). MLRM incorporate cluster‐specific random effects which allow one to partition the total individual variance into between‐cluster variation and between‐individual variation. Statistically, MLRM account for the dependency of the data within clusters and provide correct estimates of uncertainty around regression coefficients. Substantively, the magnitude of the effect of clustering provides a measure of the General Contextual Effect (GCE). When outcomes are binary, the GCE can also be quantified by measures of heterogeneity like the Median Odds Ratio (MOR) calculated from a multilevel logistic regression model. Time‐to‐event outcomes within a multilevel structure occur commonly in epidemiological and medical research. However, the Median Hazard Ratio (MHR) that corresponds to the MOR in multilevel (i.e., ‘frailty’) Cox proportional hazards regression is rarely used. Analogously to the MOR, the MHR is the median relative change in the hazard of the occurrence of the outcome when comparing identical subjects from two randomly selected different clusters that are ordered by risk. We illustrate the application and interpretation of the MHR in a case study analyzing the hazard of mortality in patients hospitalized for acute myocardial infarction at hospitals in Ontario, Canada. We provide R code for computing the MHR. The MHR is a useful and intuitive measure for expressing cluster heterogeneity in the outcome and, thereby, estimating general contextual effects in multilevel survival analysis. © 2016 The Authors. Statistics in Medicine published by John Wiley & Sons Ltd.