The data

To demonstrate and compare the applicability of different extensions of the CoxPH model to a real-world data example, injury data were obtained on 35 players from a professional rugby league club competing in the 2008 Australian NRL competition. Injury and participation data were collected from 29 matches (including all trial, fixture and finals matches). Injuries were defined as conditions associated with pain or disability that occurred during match participation, irrespective of the need for first aid, medical attention or time loss.30 In the context of this paper, a recurrent injury was said to have occurred if a person sustained more than one injury over the 29 matches, irrespective of whether or not it was to the same body region or of the same type. For this paper, only data on all contact injuries (defined as those resulting from tackling, being tackled, collision and accidental contact) were extracted. All players received a clear explanation of the study, including the risks and benefits of participation and written consent was obtained. The study was approved by the University of Queensland Human Ethics Committee.

The models

The CoxPH model and four recurrent event generalisations were applied to the sports injury data. The time variable was taken to be the match number (range 1–29). Major statistical challenges with this sort of data are how to address the number of recurrent events and the number of players at risk appropriately. Four components were considered for the recurrent event model:5 (a) risk interval which defines when a player is at risk of having an injury along a given timescale and determines whether a model is either marginal or conditional; (b) risk set or the number of players included in the set at a given point in time; (c) an event-specific or common baseline hazard; and (d) handling of within-subject correlation. Table 1 summarises how these components are defined for each of the four models considered in this paper.

The following three risk intervals were considered: (a) gap (or interoccurrence) time representing the time from the prior injury event and not relative to the actual timeline of observation; (b) calendar time which uses the same timescale for all events, referenced to a fixed point in time, but does not allow an overlap in risk periods across events for a given player; and (c) total time representing the time from the start of the player follow-up. In each case, the interval ends with the current injury. In both the gap time and calendar time representations, the player is at risk for the same length of time. Gap and calendar time models are conditional since a player is at risk of a new injury, conditioned on having sustained a previous injury. For total time, the clock does not reset for each event and the beginning of each event is at the same point in the observation timeline; risk periods for different events for the same player overlap. Total time models are marginal since the player is at risk from the start of play, independent of any previous injury. Irrespective of the definition, the risk interval for the first injury is the same.

Figure 1 describes these risk intervals in more detail, through the specific examples of three players (figure 1A). In the gap time representation, after an injury event, the player resumes play again at time 0 and the time to the next event corresponds to the number of matches that it takes for that player to experience the next event. The occurrence of all events after the first is modelled on a timescale relative to the prior event and not relative to the actual timeline of observation. Thus, the gap time (figure 1B) for our example indicates that player A is at risk of his first injury during 0–2 matches, and of his second, third and subsequent injuries during 0–18 matches, 0–6 matches and 0–1 match, respectively. In the calendar time (figure 1C), player A is at risk for his first injury event during 0–2 matches, and his second, third and subsequent injuries during matches 2–20, 20–26 and 26–27, respectively. The total time (figure 1D) indicates that player A is at risk for his first, second and subsequent injuries during 0–2, 0–20, 0–26 and 0–27 matches, respectively.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Illustrations of the risk interval formulations: (A) three players with recurrent injuries; (B) gap time; (C) calendar time; (D) total time. A circle (•) indicates an injury event and a solid square (▪) indicates censoring. Each time to an event or censoring is a separate risk interval.

The Kaplan-Meier (K-M) method is used to estimate the survival function non-parametrically from observed (censored and uncensored) survival times.31 The CoxPH model with time to first injury event as the outcome variable is a regression model used to estimate the survival probability after adjusting for both baseline hazard and explanatory variables. This model counts the players at risk at the time of this first event, after which they are no longer considered to be at risk. The result is an estimation of the probability of remaining free of injury for a given point in time based on the observed injuries. The steps in the K-M curves show changes in the probabilities of remaining free of injury for various matches across the group of players, when first injuries occur in new players.

The A-G model is a simple extension of the CoxPH model where players contribute to the risk set for an event as long as they are under observation at the time the injury occurs and share the same baseline hazard function. However, the A-G model requires the strongest statistical assumptions including that of an independent increment in which any given injury occurrence is not affected by previous injuries, that is within players, injuries are independent. This restriction means that injury dependence cannot be included and the A-G model inherently assumes that injuries do not change the player and that the player does not learn from previous injuries. Moreover, this model does not allow investigation of effects that might change based on injury-specific covariate effects, but there is the possibility of incorporating injury dependence via time-dependent covariates. Given these limitations, the A-G model is recommended when there is no injury dependence and no covariate/injury effects.

Analysis of recurrent injury data frequently assumes that the player injury histories are all statistically independent (at least conditionally on observed time-fixed covariates) so that the interoccurrence times appear in an independent manner. However, some players are intrinsically more or less susceptible to experiencing an injury than others. The frailty model is characterised by its inclusion of a random effect, or frailty, term can account for the within-player correlation between injuries. If the frailty is less than 1, a player tends to experience the injury at a later time than another player, whereas the opposite occurs if the frailty is greater than 1.

The WLW-TT model is a marginal model and assumes a common baseline hazard for all injuries within a player. Marginal models consider the marginal distribution of each failure time and impose no particular structure of dependence among distinct failure times on each player. Each recurrence is modelled as a different stratum and each stratum is treated as marginal data. This model is marginal with respect to the risk set since each player is at risk from the beginning of the study and can be at risk for several events simultaneously.

The PWP-GT model is a conditional model which allows for event dependence via stratification by event number so that different events can have different baseline hazards. The main difference to the marginal model is that a player cannot be at risk for the later injury until a prior event occurs. This conditional model preserves the order of sequential injuries in the creation of the risk set and therefore incorporates injury dependence. The PWP-GT model is estimated with the data organised in interoccurrence/gap time (ie, gap time risk set or time since the previous injury).

Model estimation and evaluation

The outcome being modelled is the probability of remaining injury-free over the 29 matches. As shown in table 1, different model formulations handle the time variable in relation to injury occurrence differently. All models were fitted using the cph function of the Design package within R (Version 2.12.2).32 ,33 The strata, cluster and frailty functions were used to fit the extended CoxPH models. The proportional hazard assumption test was performed using the cox.zph command. All models were adjusted by age, match experience and body mass of the players as known confounders of injury risk in NRL players. The R code is available from the authors, upon request.

K-M curve representations of the observed probability of remaining free of injury were used to provide a visual comparison of each model fit. The log likelihood (LL), Akaike information criterion (AIC) and Bayesian information criterion (BIC) were used to compare the goodness of fit of the fitted models in terms of fitting the observed data.34 ,35 A lower AIC or BIC indicates a better fit to the observed data and two models can be compared by comparing the differences in the AIC or BIC, with preference being given to the model with the smallest criterion measure.36 A simple rule of thumb is that models are not different if the difference in AIC is less than 2; there is minor evidence of difference when the AIC ranges from 2 to 4, and there is strong evidence for a difference with the AIC difference is more than 10. When comparing BIC, differences ≤2 are considered weak, those >2 but ≤6 are positive, those >6 but ≤10.0 are strong and BIC differences >10 are very strong.37

Model accuracy

The most common criterion for evaluating the performance of a statistical model is its accuracy in terms of data fit. In this sense, the model accuracy is an assessment of the closeness of estimates to the exact (or observed) value and can be computed on a point-by-point basis. The most widely used measures of accuracy are the mean-squared error (MSE), the root MSE, the mean absolute error and the mean absolute percentage error.38 Smaller values of each of measure indicate more accurate and reliable models. Further details about these measures can be found in Hyndman and Koehler.39

Comparing the models

Three test criteria were used to compare the fitted models: likelihood ratio (LRT), F40 and bootstrap tests.41 ,42 All tests compare two models where one model is an extension to the other (ie, the models are nested, with the simpler model being contained as a subset of the more complex one). For example, the A-G model is nested within the frailty model and comparison of the two models can test if there are random effects components for recurrent events that need to be modelled (as considered by the frailty model, but not the A-G model).

As an example, the LRT begins with a direct comparison of the likelihood scores of the two models and tests whether the frailty is necessary for analysing recurrent sports injury events. A significant LRT suggests that a random effect (frailty) accounts for the within-player correlation between injuries. A similar approach is used for the F-ratio test.

In the bootstrapping procedure, a large number of random samples are generated.41 The observed test statistic is then compared with the test statistics calculated from the bootstrap samples. Although there are many ways to use the bootstrap for hypothesis testing, the method of Walters and Campbell42 for computing a bootstrap p-value corresponding to the observed value of a test statistic T was used.

Simulation framework

Finally, a simulation approach for calculating the size and power of the models was undertaken.43 One hundred sets of injury data were simulated from the exponential distribution and the bootstrapping procedure was applied 100 times to each generated dataset to obtain the significance level of the test. Within the context of model selection, power and size estimates are based on the proportion of replications that indicate acceptable fit, with greater numbers of replications resulting in smaller CIs (higher power, more accuracy) around the estimates. Simulations were run on a bi-processed Pentium 4 machine with a 3.20 GHz processor and 2.0 Gb RAM memory. The data-generating process was performed using the SimSurv function of the prodlim package44 from R version 2.12.2,33 operating on a Windows XP professional platform.