1 Introduction

Older adults are at increased risk of experiencing poor health outcomes, such as frailty, disability, and dementia. Identifying risk factors for these outcomes is a major goal of prospectively conducted epidemiological studies of aging. Many aging-related phenotypes can only be assessed at in-person study visits. For these phenotypes, a participant’s time of onset cannot be known precisely. It is known up to the time interval between the study visit in which the outcome was first observed and the previous visit. This type of data in which the true time of onset is only known to be in a certain group is called grouped event time data, and is the focus of our paper. Grouped data arises whenever the outcome of interest is assessed at discrete study visits, so it is common not only in aging research but also in many other fields. Some examples of epidemiological studies that have collected grouped event data include the Cardiovascular Health Study, the Multiethnic Study of Atherosclerosis, and the Women’s Health and Aging Study (Fried et al. 1991; Bild et al. 2002; Onder et al. 2005).

The structure of event times for grouped data differs from that of interval censored event times. In the grouped data setting, the intervals are fixed by the number and timing of visits given by the study design, while in the interval-censored data setting the intervals are random. In the grouped data setting, all the observed event times will cluster into a small number of groups, which would correspond to the study visits, whereas in the case of interval censored data the number of intervals can be quite large.

We are interested in examining how the grouped nature of the data affects risk factor analysis in prospective cohort studies. One modeling approach is the Prentice-Gloeckler (PG) model, which accounts for the grouped nature of the data (Prentice and Gloeckler 1978). A straightforward, but less principled approach is to use the Cox proportional hazards (PH) model, the workhorse of survival analysis, substituting the unknown precise onset of the outcome with the visit time at which the outcome is first observed. The Cox model estimates the impact of covariates on the hazard function, which is based on continuous event time, whereas the observed event times in the grouped data setting are discrete. Due to the grouped nature of the event time, there will be a substantial number of ties in event times, and this needs to be addressed while using the Cox PH model. There are three well-known tie-handling techniques: the exact method, the method by Breslow, and the method by Efron. A challenge of the exact method is that it can be computationally intractable when there is a large number of ties in the data. We propose two new approximations to the exact method that are computationally efficient, even for a very large number of tied failure times. The first approximation is the Laplace approximation to an integral (Erdelyi 1956, p. 36–37), and the second is an analytic approximation to the integral.

Through simulation, we compare our two approximations to the PG model, and to the Cox PH model using Breslow’s and Efron’s tie-handling schemes. In general, the proposed methods yield estimators with low bias, and confidence intervals with the desired coverage probability. Also, we apply the proposed methods to the Cardiovascular Health Study (Fried et al. 1991), a cohort study on the development of clinical frailty in older adults. The methods give comparable results to the Prentice-Gloeckler model, in our data application. The R code for the proposed methods is provided online. The link can be found in the Supplementary Materials section at the end of this paper. Issues of competing events, although relevant given the type of events considered in our motivating example, are not addressed in this paper.

Currently, Breslow’s approximation is the default method for handling tied data in most software packages, such as Stata and SAS. Hertz-Picciotto and Rockhill (1997) and Scheike and Sun (2007) found that Breslow’s approximation tends to underestimate the coefficients in the Cox model and that Efron’s approximation performs better than Breslow’s approximation. The simulation results of this paper are consistent with those findings.

This paper is organized as follows. The assumptions are stated in Section 2. The PG model and the Cox PH Model, including the exact tie-handling method and the tie-handling methods by Breslow and Efron, are introduced in Section 3. We present our proposed approximations to the exact method in Section 4. The data application to the Cardiovascular Health Study and the simulation studies are described in Sections 5 and 6, respectively. Concluding remarks are given in Section 7.

2 Context and assumptions

We consider a population whose members have not experienced the event of interest at baseline (t = 0), which is entry into the study. For each member, we assume her/his hazard function λβ(t|Z) satisfies:

which is the Cox proportional hazards assumption. In (1), Z=(Z1,Z2,…,Zp) is a vector of covariates collected at baseline, and β=(β1,…,βp) is a vector of unknown coefficients corresponding to Z. The function of time, λ0(t), is the hazard function corresponding to Z=(0,…,0). The Cox model assumes that the event time T is absolutely continuous. Our objective is to estimate the coefficient vector β, using data from a prospective cohort study where the participants are assessed at in-person visits. During the baseline visit (t = 0), the covariates Z are collected for each participant. Let Zi denote the vector of baseline covariates (Zi1,Zi2,…,Zip) for participant i. After the baseline visit, follow-up visits are used to determine whether or not participants have experienced the event of interest. We assume that the spacing between visits is uniform across participants, though the visits themselves need not occur at regular intervals. Let m denote the number of follow-up visits. We label the baseline visit as Visit 0, and the follow-up visits as Visits 1 through m. Let tk denote the length of time between Visit 0 and Visit k. Missing data can result due to early drop-out. However, we assume no intermittent missing data, i. e. no cases where a participant missed Visit jbut attended Visit j', where j' > j.

For any given participant i, the observed data vector is (Zi,δi,Vi), where Zi is the vector of covariates, δi is an indicator of censoring, and Vi is an integer between 0 and m. A participant is censored if the study ended or the participant dropped out before she/he experienced the event. The censoring indicator δi takes the value 0 if participant i is censored, and 1 otherwise. If δi = 0, Vi is the last visit that participant i attended. If δi = 1, Vi is the first visit where participant i was observed to have experienced the event, and we refer to Vi as the “observed event time.” From the observed data, the true (i. e. precise) event time is unknown, but it is known to occur between Visits Vi−1 and Vi. We make two assumptions about the observed data. First, the data vectors (Zi,δi,Vi), i=1,…,n, are independent and identically distributed from a true, unknown joint distribution P on (Z,δ,V). Second, for participants who are censored due to dropout, their dropout time is independent of their true event time, conditional on the observed covariates.

A common attribute of grouped data is heavy ties in the observed event times. This occurs because adjacent visits can be months to years apart. For example, in the Cardiovascular Health Study, the first follow-up assessment of the clinical frailty outcome occurred three years after the baseline visit. There were 307 participants whose observed event time was the first follow-up visit, and 187 participants whose observed event time was the second follow-up visit. Define Ik as the set of participants who had the event at Visit k, i. e. Ik={i:δi=1,Vi=k}. Let dk denote the number of participants tied at Visit k, i. e. the number of elements in the set Ik.

We compare the performance of various approaches to estimate β. Existing approaches include the Prentice-Gloeckler model and the Cox model with tie-handling using the exact method, Breslow’s method, or Efron’s method. We propose two new, computationally efficient approximations to the exact method.

3 Existing approaches

4 Proposed approximations to the exact tie-handling method

As discussed in Section 3, computing the exact likelihood is difficult when the tie sizes are not small. For any given participant i and Visit k, where k∈{1,…,m}, define

For conciseness, we refer to this as ai,k throughout the rest of the paper. DeLong et al. (1994) showed that the contribution to the exact likelihood at Visit k (i. e. the average of the dk! possibilities for Visit k) can be represented as dk! multiplied by

Our proposed method is to approximate (14), using two new approaches: the Laplace approximation and a simple analytic approximation. We also considered a Gaussian quadrature approach, the Gauss-Laguerre quadrature (Davis and Rabinowitz 1984, p. 222–223). However, it did not work as well as the other approaches discussed here, mainly due to the difficulty in choosing the appropriate number of quadrature points (results not shown).

5 Application to the cardiovascular health study

The Cardiovascular Health Study (CHS) was a study of 5,201 older adults enrolled between 1989 and 1990, and followed up at annual clinic visits until 1999 (Fried et al. 1991). One of the outcomes collected in CHS was clinical frailty, which is a syndrome defined as the presence of three or more of the following: unintentional weight loss, exhaustion, low grip strength, slow walking speed, and low physical activity (Fried et al. 2001). Frailty status was assessed at the baseline visit and at two follow-up visits, occurring three and seven years after the baseline visit, respectively. We fitted the Cox and PG models to time-to-frailty data collected in CHS. It should be noted that, in addition to the cohort of 5,201 participants, a second cohort of 687 participants was recruited in 1992. We focus on the first cohort in our analysis. This is because the follow-up schedule for frailty differs between the two cohorts. The second cohort had only a single follow-up assessment of frailty.

In our analysis of time-to-frailty, five covariates were used: age, gender (0 = female, 1 = male), whether the person had ever smoked (0 = no, 1 = yes), whether the person had received education beyond high school (0 = no, 1 = yes), and C-reactive protein (CRP) level. These covariates were collected at the baseline visit. We standardized both age and CRP level. Age, in units of years, was standardized by subtracting 65 and dividing by 10. CRP level, in units of mg/L, was standardized by dividing by the interquartile range (2.15 mg/L). The fitted models included main terms for the covariates and no interaction terms. We let δ be an indicator of whether the participant was assessed as frail during her/his follow-up (0 if no, 1 if yes). If δ = 0, we let V be the last follow-up visit at which the participant’s frailty status was assessed. If δ = 1, we let V be the first follow-up visit at which the participant was assessed as frail.

We included only participants who satisfied the following conditions: (i) were non-frail at the baseline visit, (ii) frailty status was assessed at the first follow-up only, or at both follow-ups, (iii) covariates were non-missing. The resulting sample size was n = 3,605. The purpose for (i) was, since we sought to study the time-to-frailty development after the baseline visit, we focused on those who were not yet frail at the baseline visit. Criterion (ii) ensured that there was some follow-up on the frailty status after the baseline visit because those without any such follow-up would not affect the model fits. Also, criterion (ii) excluded those whose frailty status was observed at the second follow-up but not the first follow-up. If the status at the second follow-up was frail, then δ = 1 but it could not be deduced whether V should be 2 or 1 (i. e. the person was actually frail by the first follow-up). If the status at the second follow-up was non-frail, there was also some ambiguity since the person might have been non-frail through the second follow-up, or have been frail at the first follow-up and transitioned to being non-frail by the second follow-up. A potential impact of excluding participants who did not attend the first follow-up visit is that the sample used in our analysis may under-represent a specific subtype of participants. For example, the participants who missed the first follow-up visit may be more likely to have severe health problems, e. g. some of them might have missed their study visit due to their health problems. Criterion (iii) was required since participants without the full set of covariates would not impact the model fits.

Web Figure 1 in the Supplementary Materials presents the empirical marginal distributions of age, gender, smoking status, education level, and CRP level, after excluding participants as described above. The mean age is 72 (range 64–95) and the proportion of female participants is 0.58. Table 1 presents the empirical joint distribution on (δ,V), after excluding participants as described above. Each cell gives the number of participants, followed by the corresponding sample proportion in parentheses. We evaluated the validity of the proportional hazards assumption by testing the null hypothesis that the correlation between the scaled Schoenfeld residuals and time is zero, separately for each covariate. The results supported the proportional hazards assumption; the p-value was 0.10 for age, 0.77 for education level, 0.48 for gender, 0.33 for smoking status, and 0.80 for CRP level.

Table 2 presents the regression coefficient estimates in the top panel, with the corresponding standard error estimates in the middle panel. The values are shown up to four decimal places. The results in this table and those throughout the paper were computed using R. When fitting the Cox model, we applied our proposed tie-handling methods Laplace and Analytic, as well as the existing methods by Efron and Breslow. When assessing each of these methods, we use PG as a reference since it accounts for the discreteness of the data directly. All of the methods were computationally efficient; in this data application, the running time was 0.01 seconds for Efron’s method, 0.02 seconds for Breslow’s method, 0.3 seconds for PG, 1.2 seconds for Analytic, and 1.4 seconds for Laplace.

In the upper two panels of Table 2, the results for Laplace perfectly match those for PG. For all covariates, the absolute difference between the Laplace estimate and the corresponding PG estimate was on the order of 10−6 or less. This applied to both the regression coefficient estimates and standard error estimates. For each covariate, the magnitude of the percent difference of Laplace compared to PG, i. e. |(Laplace - PG)/PG× 100|, was at most 0.01 %. The results for Analytic and Efron’s method were also close to those of PG. For both Efron’s method and Analytic, the absolute difference from PG was at most 6 ×10−3 in magnitude. For each covariate, the magnitude of the percent difference from PG was at most 0.8 % for Efron’s method and 0.6 % for Analytic. Compared to PG and Laplace, the regression coefficients for Efron’s method and Analytic were closer to the null. Those for Breslow’s method were further in the direction of the null. Breslow’s method yielded the most different results from PG. For the regression coefficient estimates, the magnitude of the percent difference of Breslow’s method compared to PG was as much as 8 %.

The bottom panel of Table 2 shows the 95 % confidence intervals for the exponentiated regression coefficients computed using each approach. For interpretability, we present the confidence intervals for the exponentiated regression coefficients rather than for the regression coefficients themselves. Each confidence interval is computed by taking the corresponding regression coefficient estimate, adding and subtracting the quantity 1.96 multiplied by the standard error estimate, and exponentiating both values. Beside each 95 % confidence interval, we also present the corresponding regression coefficient estimate (after exponentiation). Although there are some differences between the confidence intervals computed by the different approaches, the qualitative conclusions are the same. The various approaches all find that a person’s hazard for frailty increases with higher age, holding the other factors constant. Also, the hazard for frailty is lower for those who have education beyond high school than for those who have a high school education or less. Females have a higher hazard than males, holding other factors constant. All of the approaches concluded that the hazard for frailty is not statistically different between those who have ever smoked versus those who have never smoked, but the point estimate suggests that smoking raises the hazard. The hazard for frailty increases with higher CRP level, holding other factors constant.

6 Simulation studies

Through simulation, we evaluate the Cox and PG models in the grouped data setting. For the Cox model, the Laplace, Analytic, Efron’s, and Breslow’s tie-handling methods are separately applied. We compare the regression coefficient estimators with respect to bias, standard error, and root-mean-square error. Also, we evaluate the actual coverage probability of the nominal 95 % confidence interval for each method.

Each simulation study includes a large number of simulations, nsim=10,000. For any given simulation s, a study with sample size n is simulated. The contents of the study data set are described in the following paragraph. Using the data set, the parameter β is estimated using the Cox model and the PG model, respectively. For a given estimator βˆ=(βˆ1,…,βˆp), we approximate the bias, standard error, and root-mean-square error of the i-th element βˆi as:

where βˆs,i is the estimate of βi in the s-th simulation. Also, we estimate the coverage probability of the 95 % confidence interval for βi by

where SEˆs,i is the estimate of the standard error of βˆi from the s-th simulation.

Each simulation follows n participants according to a schedule which includes a baseline visit (Visit 0) and m follow-up visits (Visits 1,…,m), occurring at fixed times after the baseline visit. The collected data set includes covariate, censoring, and time information of the n participants, designated as (Zi,δi,Vi), i=1,…,n. These random variables were previously defined in Section 2. For any given i, participant i does not miss any visits preceding visit Vi. The simulation studies were organized into two sets, based on the structure of the data generating distribution. We refer to the two sets as the Single Covariate Simulations and the Data-Based Simulations.

7 Discussion

The exact method of tie-handling is computationally intractable when there is a large number of ties in the observed event times. Hence, we propose two, new and computationally efficient approximations, Laplace and Analytic. Both approximations performed well in the simulation studies and a real data application to the Cardiovascular Health Study. The Laplace approximation yielded results that were almost identical to that of the PG model for all scenarios. Laplace and PG had low bias and generated confidence intervals with approximately the nominal coverage probability. Our results provide strong evidence that the Cox model, using the Laplace approximation for tie-handling, can serve as a substitute for the PG model in the grouped data setting. The Analytic approximation, given its remarkable simplicity, performed reasonably well in most scenarios, although it tended to have substantial bias for large β, as shown in the Single Covariate Simulations when β = 2. Hence, Analytic may be an attractive option in grouped event time data where the relative risks are not expected to be large. In fact, in the CHS application, the Analytic approximation provided results very similar to PG and Laplace. Our new approximations use a different approach than Breslow’s and Efron’s approximations because the former approximate the exact likelihood with ties, while the latter approximate the actual likelihood. A potential advantage of our approach compared to using the PG model is that fitting the PG model for high-dimensional data might be computationally intensive, relative to our simple approximations.

Both PG and Cox models make the proportional hazards assumption. Although the Cox model estimates the baseline hazard nonparametrically, and the PG model estimates it as the intercept terms, these two are essentially equivalent (with a complementary log-log link in the PG model). Thus, the essential comparison is in terms of different ways of handling the regression part, i. e. the covariate part of the model.

A potential area of future research is to examine the accuracy of our methods under cases with more covariates. Another potential area of future research is to determine the impact of a large number of visits on the results. In this paper, we considered the case of a small number of follow-up visits, e. g. up to 3 follow-up visits in the simulation studies and 2 in the CHS data example.

We would not recommend the use of Breslow’s tie-handling method in the grouped data setting. In simulation, it was observed that Breslow’s method had larger bias than the other methods and can have poor coverage probability, as low as 0.5. In the CHS application, the results of Breslow’s method were also quite different from the other methods.

For those who use the R programming language, we note that the exact method implemented in the “survival” package of R is not the exact method for tie-handling discussed in this paper. The exact method in the “survival” package corresponds to the discrete logistic model proposed by Cox (1972). The discrete logistic model has a different parameter of interest than the models considered in this paper. This is because, rather than making the proportional hazards assumption given in (1), the discrete logistic model assumes that the odds ratios of the hazards are proportional, i. e. logitλβ(t|Z)=logitλ0(t)exp(Z⋅β). Thus, the term exp(Z⋅β) represents the increase in the odds ratio, rather than the increase in the hazard. This model simplifies to the proportional hazards model in (1), if the event times are continuously distributed (Cox 1972).

The strong performance of the Cox PH model, using Laplace, Analytic, and Efron’s tie-handling methods, was an impressive and highly surprising finding, given the extremely discrete nature of the event times that we considered. There were 307 and 187 tied event times, respectively, at the first and second follow-up visits of frailty in the CHS study. One would not think that the Cox PH model would be appropriate for estimating relative risks in such a situation. To our knowledge, this type of result has not been reported in the literature before.