Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 https://doi.org/10.1186/s12874-022-01503-1 RESEARCH Open Access Left-censored recurrent event analysis in epidemiological studies: a proposal for when the number of previous episodes is unknown Gilma Hernández‑Herrera1,2, David Moriña3,4,5* and Albert Navarro6,7 Abstract Background: When dealing with recurrent events in observational studies it is common to include subjects who became at risk before follow-up. This phenomenon is known as left censoring, and simply ignoring these prior epi‑ sodes can lead to biased and inefficient estimates. We aimed to propose a statistical method that performs well in this setting. Methods: Our proposal was based on the use of models with specific baseline hazards. In this, the number of prior episodes were imputed when unknown and stratified according to whether the subject had been at risk of present‑ ing the event before t = 0. A frailty term was also used. Two formulations were used for this “Specific Hazard Frailty Model Imputed” based on the “counting process” and “gap time.” Performance was then examined in different sce‑ narios through a comprehensive simulation study. Results: The proposed method performed well even when the percentage of subjects at risk before follow-up was very high. Biases were often below 10% and coverages were around 95%, being somewhat conservative. The gap time approach performed better with constant baseline hazards, whereas the counting process performed better with non-constant baseline hazards. Conclusions: The use of common baseline methods is not advised when knowledge of prior episodes experienced by a participant is lacking. The approach in this study performed acceptably in most scenarios in which it was evalu‑ ated and should be considered an alternative in this context. It has been made freely available to interested research‑ ers as R package miRecSurv. Background how often) at study inception, which is problematic when In epidemiological cohort studies, participants may have the baseline hazard of the event is time dependent. There been at risk of the event of interest before entering follow- are two relevant situations to consider. First, if an event up, especially in observational designs. This can lead to can only occur once and has already happened, the result unawareness of their history, specifically the time at risk, is already determined and will require specific statistical and whether they have experienced the event (including techniques for analysis. Second, if the outcome of interest is a recurrent event that can occur more than once, the number of prior episodes will be unknown despite the *Correspondence: dmorina@ub.edu 5 Facultat d’Economia i Empresa, Universitat de Barcelona (UB), Avinguda potential for new events. This represents left censoring Diagonal, 690‑694, 08034 Barcelona, Spain where the censored variable is of the discrete type and for Full list of author information is available at the end of the article which we can define different baseline hazards. © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​ mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 2 of 9 There is a need to clarify how data are handled when episode is always posterior at the end of the k-1st. The the prior history is unknown in a cohort and the outcome hazard function is shown below (Eq. 1). In the second of interest is a recurrent event with event dependence. format, called “gap time,” time is always considered in Specifically, this concerns situations where we know the relation to the prior episode such that the start of each moment from when individuals become at risk, but we new episode for a given individual is set at zero (Eq. 2) do now know the number of prior episodes. For example, concerning the risk of sick leave in a work force, we will ik (t) = 0k (t)eXi β hik (t) = h0k (t)eXi β (1) likely know the start date for employment; however, espe- cially for people with ample trajectory, we may not know hik (t) = h0k t − tk−1 eXi β (2)   how many sick leaves they had prior to this employment. Another example can be seen in cohorts where the out- where Xi represents the vector of covariables, β repre- come is the incidence of infection with the human papil- sents the coefficients of regression, k is the kth event epi- loma virus in adult women. It would be relatively simple sode for individual i, and h0 k(t) is the function of the to know how long they have been at risk (first sexual baseline hazard dependent on k. intercourse), but we will not know the true number of If the phenomenon under study lacks event depend- infections because they are typically asymptomatic. ence, these models can be simplified to a function of the This study has two objectives. The first is to describe common baseline hazard. These are known as Andersen– an analysis proposal for recurrent phenomena in the Gill models [13] and assign the same baseline hazard presence of event dependence when the prior history independent of the episodes that have already been expe- is unknown for some or all subjects. The second is to rienced (Eq. 3). compare the performance of our proposed model with another that ignores event dependence. hi (t) = h0 (t)eXi β (3) where h0(t) is the common baseline hazard. Methods Event dependence Individual heterogeneity When analyzing a recurrent event, we often observe If a model is perfectly specified, with all possible covari- a phenomenon called event dependency in which the ates accounted for, then the baseline hazard reflects the baseline hazard of an episode depends on the num- randomness of the event time given the value of the ber of episodes that have already occurred. To date, the covariates. In practice, however, it is rarely possible to phenomenon has been estimated for falls [1], sickness account for all relevant covariates [14]. This requires absence [2, 3], hospitalizations in heart failure [4], and another aspect to be considered: the unmeasured effects cardiovascular readmission after percutaneous coronary produced by between-subject variability, presumably due intervention [5]. All of these have shown that the baseline to unobserved exposures. This phenomenon is called hazard increased most significantly with the number of individual heterogeneity and in practice is analyzed by prior episodes. adding a frailty to the model, νi, making an individual Recurrent event analyses in epidemiological studies are random effect to account for the “extra” variability. Given well described [6–9], typically as extensions of Cox’s clas- that νi has a multiplicative effect, we can imagine that it sical model of proportional hazards [10, 11]. Specifically, represents the cumulative effect of one or more omit- baseline hazard models, also known as conditional mod- ted covariates [15]. The most commonly adopted frailty els or Prentice, Williams, and Peterson (PWP) models terms are E [νi] = 1 and V [νi] = θ [16–18]. Thus, the [12], are used to consider the existence of event depend- models specified in Eq. 1, Eq. 2, and Eq. 3 change to those ence. These models assume that the baseline hazard of an specified in Eq. 4 for “Conditional Frailty Model - Count- episode differs as a function of episodes that have already ing Process,” Eq. 5 for “Conditional Frailty Model - Gap occurred, stratifying by how many there have been. This Time,” and Eq. 6 for “Shared Frailty Model,” respectively. allows general or specific effects to be calculated for each episode, with all at-risk individuals included in the first hik (t) = νi h0k (t)eXi β (4) strata, but only those with an episode in the previous strata subsequently considered at risk. hik (t) = νi h0k t − tk−1 eXi β (5)   PWP models can be formulated in two ways depend- ing on the risk interval used (i.e., how time is considered) [6]. In the first, called “counting process,” time is handled hi (t) = νi h0 (t)eXi β (6) typically for the survival analysis, always referenced from the start of follow-up such that the beginning of the kth Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 3 of 9 The problem of being unaware of the previous history episodes, we will mix individuals that have different for some individuals baseline hazards at that instant (e.g., risk starts for Sub- Specific baseline hazard methods, either with or without ject 2 at baseline, when Subject 1 is already at risk of frailties, can be applied when all required information a third). Second, we mix temporal scales: on one side is known, particularly the number of prior episodes for there will be subjects (e.g., Subject 2) whose follow-up each individual. In practice, however, this information time corresponds with their time at risk, but on the is not always available, preventing the inclusion of basic other side there will be subjects (e.g., Subject 1) whose information needed to consider event dependence. follow-up time does not correspond with their time Figure 1 shows representative data for two subjects at risk. Consequently, two subpopulations with differ- with their respective risk intervals. These are shown ent baseline hazards are mixed, which goes against the according to the counting process on the top and accord- assumptions of the Cox model and its extensions for ing to gap time on the bottom. Notice that the difference recurrent phenomena. between these formulations is that the gap time “restarts” When there is missing data for prior history and spe- time to risk at t = 0 after each episode. When starting our cific baseline hazard models are unsuitable, an alter- study, the first subject (id = 1) had two prior episodes native approach is to ignore that history and fit the and had been at risk for a considerable amount of time, models based on common baseline hazards (i.e., eqs. 3 making this a left-censored observation. However, we are and 6). However, these assign the same baseline hazard only aware of the two episodes that occur from the start to all episodes, and as such, they consider neither the of follow-up. By contrast, for the second subject (id = 2), possible effect of the number of prior episodes nor the exposure starts at t = 0, they have an episode at t = 5, and effect of comparing two subjects at the same instant follow-up stops at t = 7. The data tables show the result- who have times to risk that can be radically different. ing analysis: note that the prior history for Subject 1 The use of such models with common baseline hazards “disappears” in both approaches such that the risk only for analyzing recurrent phenomena with event depend- appears at the same instant as it does for Subject 2. ence has been shown to be highly inefficient, generating Thus, two problems exist. First, if there is event high levels of bias in parameter estimates and confi- dependence, and we fail to stratify by number of dence intervals. Results can be unsound, even if when Fig. 1 Example history of two individuals presented according to the counting process and gap time formulation. The counting process is shown in the top image and table, and the gap time approach is shown below for two patients (id 1, 2) at time points (t) 0, 3, 5, and 7 Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 4 of 9 the event dependence is of small intensity [19]. Thus, rules [23] and based on the following steps in a Bayesian alternative analytic methods are necessary for this context: context. 1. Fit the corresponding COMPoisson model   and find Proposal the posterior mean β̂ and variance V β̂ of model Our proposal starts with the assumption that, even parameters β.    though the history of a given subject may be unknown, 2. Draw new parameters β∗ from N β̂,V β̂ . we do know who was at risk before starting follow-up. 3. Draw imputations from the COMPoisson distribu- This is based on three fundamental considerations: 1) tion using the scores obtained in the previous step as imputing k, the number of previous episodes for those parameters. subjects at risk before follow-up starts; 2) treating the subpopulation “previously at risk” separately from that A detailed description of the whole imputation process “not previously at risk” and 3) using a frailty term to cap- can is reported by Hernández et al. (2020) [24]. ture the impact of unobserved effects, including uncer- tainty related to the imputation process. This gives two formulations that we call “Specific Hazard Frailty Model Simulation study Imputed,” as shown in Eq. 7 for the counting process Six populations were simulated through the R package, (SPECIFIC.CP) and Eq. 8 for gap time (SPECIFIC.GT). survsim [25], corresponding to previously described cohorts of workers [19, 26]. The first three simulated hikr (t) = νi h0kr (t)eXi β (7) worker history depending on the occurrence of long- term sick leave with several intensities of event depend- hikr (t) = νi h0kr t − tk−1 eXi β (8) ence. Given that the simulations are carried out through   Weibull distributions with ancillary parameters equal where k is the number of previous episodes in individual to 1 (that is, exponential distributions) the hazard i when they are known or the imputed value when they was deemed constant within each episode, but differ- are not known; r indicates the “previously at risk” or “not ent between episodes (in other words, there was event previously at risk” subpopulation to which the individual dependence). For the fourth, fifth, and sixth populations, belongs. the outcome was sick leave according to diagnosis (i.e., In practice, this proposal means that we stratify by the respiratory system, musculoskeletal system, and mental interaction between risk prior to the start of follow-up and behavioral disorders). In this case, the hazard func- and the number of prior episodes. Therefore, the use of tions were not constant within the episode. the individual random error term υi intends to capture Table 1 shows the characteristics of each episode in the error that will be made when imputing, as well as the each population. The maximum number of episodes that effect of any variable with a non-nil effect not being con- a subject may suffer was not fixed, but the baseline haz- sidered in the analysis. Stratifying by the number of prior ard was considered constant when k ≥ 3. The exposure is episodes intends to safeguard event dependence: doing represented by covariates X1, X2, and X3, with Xj ≈ Ber- it as an interaction with the fact that it is an individual noulli (0.5) and β1 = 0.25, β2 = 0.5, and β3 = 0.75 repre- previously at risk, or not, separates the two subpopula- senting effects of different magnitudes set independently tions to avoid mixing times that are not comparable on of the episode k to which the worker is exposed. the same scale. Four different situations, which were combinations To impute the prior history, a generalized linear model between two possible follow-up times (2 and 5 years) (GLM) based on the Conway–Maxwell Poisson distribu- and two maximum times at risk prior to the beginning tion (COMPoisson) [20] is fitted to the observed number of the cohort (2 and 10 years), were simulated for each of events. Imputed values are randomly sampled from the population. For example, for 10 years’ maximum time corresponding distribution with the parameters obtained at risk and 5 years’ follow-up, we generated dynamic in the previous step, including random noise generated populations with 15 years’ history and selected sub- from a normal distribution. To produce proper estima- jects who were either present at 10 years’ follow-up or tion of uncertainty, the described methodology is used incorporated after that date. Follow-up time was then within a multiple imputation framework. The data in the re-scaled, setting t = 0 at 10 years for subjects already Results section reflect the combination of 5 imputed data present in the population, and setting t = 0 as the begin- sets (following the classic advice to use a low number ning of the follow-up period for those incorporated of imputed data sets and considering that there is little later. This procedure allowed us to obtain a cohort in to gain from using more [21–23]), according to Rubin’s which some subjects had worked up to 10 years and Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 5 of 9 Table 1 Characteristics of the simulated populations ir (t) = νi 0r (t)eXi β (9) Episode Distribution β0 Ancillary HR Additionally, a second simulation study was conducted 1 Weibull 8.109 1 1 with the same technical settings but with m = 1000 simu- 2 Weibull 7.927 1 1.20 lation instances to avoid spurious fluctuations and with ≥3 Weibull 7.745 1 1.44 only one covariate with no effect (β = 0), to assess type I 1 Weibull 8.109 1 1 error rates for rejecting the null hypothesis of no expo- 2 Weibull 7.703 1 1.50 sure effect. Table 2 presents the performance measures ≥3 Weibull 7.298 1 2.25 used for evaluating the different methods. 1 Weibull 8.109 1 1 2 Weibull 7.193 1 2.50 Results ≥3 Weibull 6.276 1 6.25 The presented results refer to n = 1000 and 5 years’ fol- 1 Log-normal 7.195 1.498 1 low-up because the differences observed for n = 250, 2 Log-logistic 6.583 0.924 1.77 n = 500, and 2 years’ follow-up were not considered rele- ≥3 Weibull 6.678 0.923 2.53 vant (see Supplementary Material). Regarding bias, Fig. 2 1 Log-logistic 7.974 0.836 1 highlights that the COMMON model only obtained val- 2 Weibull 7.109 0.758 3.81 ues under 10% in population 1 and in some cases when ≥3 Log-normal 5.853 1.989 7.19 100% of individuals were at prior risk. Models SPECIFIC. 1 Log-normal 8.924 1.545 1 CP and SPECIFIC.GT generally obtain biases < 10% in 2 Log-normal 6.650 2.399 10.13 most situations, except for SPECIFIC.GT in populations ≥3 Log-normal 6.696 2.246 11.19 4 and 5, which is slightly above this threshold, and for −t p SPECIFIC.CP, especially in population 3 when 100% of Weibull distribution: f (t) = pt p−1 e ,  = e−pβ0  1  −1 {log(t)−µ}2 individuals had prior risk at t = 0, which seemed more Lognormal distribution: f (t) = √ e 2σ 2 ,µ = β0 tσ 2π sensitive. For the first three populations, SPECIFIC.GT 1/γ t 1/γ −1 Loglogistic distribution: f (t) =  2 ,  = e −β 0 shows equal or less bias than SPECIFIC.CP, whereas the γ 1+(t)1/γ converse is true for the second three populations, pro- vided that risk starts during follow-up for at least 50% of were treated as unknown with 5 years of effective fol- the cohort. low-up (observed between 10 and 15 years in the origi- The average 95%CIs are largest in the COMMON nal simulated follow-up). In each case, m = 100 samples model, excluding population 1, when up to 30% of indi- of sizes 250, 500, and 1000 were simulated with differ- viduals were at previous risk, always overcoming the ent proportions of subjects at risk prior to entering the accuracy of the SPECIFIC.CP model in populations 1, cohort (0.1, 0.3, 0.5, 0.75 and 1). 2, and 3, with the SPECIFIC.GT model most accurate in Our proposal was compared to a model with frailty populations 4, 5, and 6 (Fig. 3). and common baseline hazard in terms of the number The 95%CI coverage of the COMMON model was of previous events, but with different risk per subpop- clearly under 95% in most cases, decreasing as the event ulation. We named this the “Common Hazard Frailty dependence increases, and for the first four populations, Model with stratification by subpopulation” (COM- the proportion at risk prior to t = 0. The coverages for MON). Importantly, this model does not take event models SPECIFIC.CP and SPECIFIC.GT were gener- dependence into account, as expressed in Eq. 6, but ally acceptable, even when excessively conservative in separates individuals by prior risk: some cases (> 95%). The SPECIFIC.CP model failed in Table 2 Performance evaluation criteria Simulation study 1 Evaluation criteria (covariate j = 1, 2, 3 and sample m = 1, …, 100) Mean relative bias 3 100 |β̂jm −βj | ×100 j=1 m=1 βj 3×100 Average length of 95%CI 3 100   j=1 m=1 2Z1−α/2 SE β̂jm 3×100 Coverage   Percentage of times the 95% confidence interval β̂jm ± Z1−α/2 SE β̂jm includes βj Simulation study 2 Evaluation criterium Type I error Percentage of times the H0:β = 0 is rejected, for m = 1, …, 1000 Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 6 of 9 Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 60 40 20 Mean relative bias (%) 0 Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 60 40 20 0 0.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.00 Proportion of subjects at risk before t=0 Model COMMON SPECIFIC.CP SPECIFIC.GT Fig. 2 Bias according to population and maximum time at risk prior to the beginning of the cohort Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 0.8 Average length of the 95% confidence interval 0.6 0.4 Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 0.8 0.6 0.4 0.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.00 Proportion of subjects at risk before t=0 Model COMMON SPECIFIC.CP SPECIFIC.GT Fig. 3 Average length of the 95% confidence interval according to population and maximum time at risk prior to the beginning of the cohort population 3 when 100% had prior risk, due to the high coverage (92.5–97.5%) criteria. Considering all situa- bias, whereas the SPECIFIC.GT model obtained coverage tions, the COMMON model rarely met these criteria, < 80% even in population 4 when < 100% of individuals while the SPECIFIC.CP model showed the highest level were at prior risk (Fig. 4). of compliance. Tables S1 and S2 (Supplementary Material) summa- Figure 5 shows the type I error distribution for each method, rize the results using arbitrary relative bias (< 10%) and with acceptable rates of around 5% observed in all cases. Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 7 of 9 Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 100 80 60 40 Coverage (%) 20 Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 100 80 60 40 20 0.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.000.10 0.30 0.50 0.75 1.00 Proportion of subjects at risk before t=0 Model COMMON SPECIFIC.CP SPECIFIC.GT Fig. 4 Coverage according to population and maximum time at risk prior to the beginning of the cohort Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Max(t) before t0 = 2 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 6 4 Type I error rate (%) Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Max(t) before t0 = 10 years Population 1 Population 2 Population 3 Population 4 Population 5 Population 6 6 4 0.25 0.50 0.75 1.00 0.25 0.50 0.75 1.00 0.25 0.50 0.75 1.00 0.25 0.50 0.75 1.00 0.25 0.50 0.75 1.00 0.25 0.50 0.75 1.00 Proportion of subjects at risk before t0 Model COMMON SPECIFIC.CP SPECIFIC.GT Fig. 5 Type I error rate according to population and maximum time at risk prior to the beginning of the cohort Discussion should be to describe the baseline hazard form of the Our results indicate that not knowing the number of phenomenon under study in each episode. The SPE- prior episodes an individual has experienced should not CIFIC.GT model should be selected if this is a function justify the use of models with a common baseline hazard, of constant risk, while the SPECIFIC.CP model would which in most situations show higher bias and less cover- also be acceptable when information about the number age than specific baseline hazard models. The first step, of prior episodes is lacking in up to 30–50%; however, the which is indispensable for choosing the correct method, SPECIFIC.CP model should be selected for phenomenon Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 8 of 9 ruled non-constant when there is a lack of prior informa- process; GLM: Generalized linear model; GT: Gap time; HR: Hazard ratio; PWP: Prentice, Williams, and Peterson model; SPECIFIC.CP: Specific hazard frailty tion in up to 50% of cases. Neither model is sufficiently model imputed (Counting process formulation); SPECIFIC.GT: Specific hazard robust if the number of previously experienced episodes frailty model imputed (Gap Time formulation). is unknown for all individuals, even though the SPE- CIFIC.GT model could be adequate in some cases. Supplementary Information Some readers may be surprised by the finding that The online version contains supplementary material available at https://​doi.​ bias for the COMMON models in populations 5 and 6 org/​10.​1186/​s12874-​022-​01503-1. decreased when increasing the proportion of workers at risk before t = 0. This may be because these populations Additional file 1. present high levels of recurrence: thus, more subjects are Additional file 2. at risk prior to t = 0 meaning that more subjects will have already had ≥2 episodes (and will be imputed with val- Acknowledgements ues ≥2). Given that our simulation set had the same base- Not applicable. line hazard from the second episode onwards (a realistic Authors’ contributions convergence phenomenon), many subjects will have had G.H-H. and A.N. conceived the theoretical framework, and D.M. and G.H-H. the same baseline hazard and the COMMON model will designed the simulation study and participated in the analysis of the data. All have performed better. authors have read and agreed to the published version of the manuscript. Other unexpected results are seen with the partially Funding non-monotonous bias and CI coverage functions associ- Not applicable. ated with increasing the proportion of missing data (e.g., Availability of data and materials the SPECIFIC.CP and SPECIFIC.GT models in popula- R codes used in the described analyses are available as supplementary mate‑ tion 1). These situations were in the minority, and the rial and in the GitHub repository https://​github.​com/​dmori​nya/​miRec​Surv. reader should keep in mind that the actual differences were very low despite appearing relevant (the Y-axis Declarations scales in some figures are very accurate). The results for Ethics approval and consent to participate n = 250 and n = 500 (Supplementary Material) did not Not applicable. follow the same pattern, which indicates that these differ- ences are probably spurious. Consent for publication Not applicable. Finally, it is important to note that the effectiveness of imputing prior missing episodes to account for left-cen- Competing interests sored data depends strongly on the assumption of ‘miss- The authors declare no competing interests. ing at random’ (i.e., that measured covariates actually Author details predict the episode history of an individual). Therefore, 1 Instituto de Investigaciones Médicas, Facultad de Medicina, Universidad a correctly specified imputation (and similarly, an appro- de Antioquia, Medellín, Colombia. 2 Methodology of Biomedical Research and Public Health, Autonomous University of Barcelona, Cerdanyola del priate final survival model specification) is essential to Vallès, Spain. 3 Department of Econometrics, Statistics and Applied Econom‑ obtain unbiased and efficient estimates. ics, Riskcenter‑IREA, University of Barcelona (UB), Barcelona, Spain. 4 Centre de Recerca Matemàtica (CRM), Cerdanyola del Vallès, Spain. 5 Facultat d’Economia i Empresa, Universitat de Barcelona (UB), Avinguda Diagonal, 690‑694, Conclusions 08034 Barcelona, Spain. 6 Psychosocial Risks, Organization of Work and Health (POWAH), Autonomous University of Barcelona (UAB), Cerdanyola del Vallès, This simulation study reveals that simply ignoring prior Spain. 7 Biostatistics Unit, Faculty of Medicine, Autonomous University of Barce‑ episodes of the outcome of interest and using common lona (UAB), Cerdanyola del Vallès, Spain. baseline hazard models (the most common approach in Received: 16 July 2021 Accepted: 3 January 2022 epidemiological cohort studies) can lead to extremely biased and inefficient estimates. In turn, this will gener- ate inaccurate conclusions, especially when the amount of missing information is high. However, the perfor- References mance of these estimates can be improved by using the 1. Navarro A, Ancizu I. Analyzing the occurrence of falls and its risk factors: proposed methodology, which is already available on some considerations. Prev Med (Baltim). 2009;48:298–302. 2. Navarro A, Reis RJ, Martin M. Some alternatives in the statistical analysis of standard software and ready to be used by any interested sickness absence. Am J Ind Med. 2009;52(10):811–6. researcher as R package miRecSurv [27]. 3. Reis RJ, Utzet M, La Rocca PF, et al. Previous sick leaves as predictor of subsequent ones. Int Arch Occup Environ Health. 2011;84:491–9. 4. Braga JR, Tu JV, Austin PC, et al. Recurrent events analysis for examination Abbreviations of hospitalizations in heart failure: insights from the enhanced feedback COMMON: Common hazard frailty model with stratification by subpopulation; for effective cardiac treatment (EFFECT) trial. Eur Hear J Qual Care Clin CI: Confidence interval; COMPoisson: Conway-Maxwell Poisson; CP: Counting Outcomes. 2018;4:18–26. Hernández‑Herrera et al. BMC Medical Research Methodology (2022) 22:20 Page 9 of 9 5. Vasudevan A, Choi JW, Feghali GA, et al. Event dependence in the analysis of cardiovascular readmissions postpercutaneous coronary intervention. J Investig Med. 2019;67:943–9. 6. Kelly PJ, Lim LL. Survival analysis for recurrent event data: an application to childhood infectious diseases. Stat Med. 2000;19:13–33. 7. Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. New York: Springer; 2000. 8. Amorim LDAF, Cai J. Modelling recurrent events: a tutorial for analysis in epidemiology. Int J Epidemiol. 2015;44:324–33. 9. Box-Steffensmeier JM, De Boef S. Repeated events survival models: the conditional frailty model. Stat Med. 2006;25:3518–33. 10. Cox DR. Partial likelihood. Biometrika. 1975;62:269–76. 11. Cox DR. Regression models and life table. J R Stat Soc. 1972;B34:187–220. 12. Prentice R, Williams B, Peterson A. On the regression analysis of multivari‑ ate failure time data. Biometrika. 1981;68:373–9. 13. Andersen PK, Gill RD. Cox’s regression model counting process: a large sample study. Ann Stat. 1982;10:1100–20. 14. Balan TA, Putter H. A tutorial on frailty models. Stat Methods Med Res. 2020;29:3424–54. 15. O’Quigley J, Stare J. Proportional hazards models with frailties and ran‑ dom effects. Stat Med. 2002;21:3219–33. 16. Rondeau V, Commenges D, Joly P. Maximum penalized likelihood estima‑ tion in a gamma-frailty model. Lifetime Data Anal. 2003;9:139–53. 17. Balakrishnan N, Peng Y. Generalized gamma frailty model. Stat Med. 2006;25:2797–816. 18. Govindarajulu US, Lin H, Lunetta KL, et al. Frailty models: applications to biomedical and genetic studies. Stat Med. 2011;30:2754–64. 19. Navarro A, Casanovas G, Alvarado S, et al. Analyzing recurrent events when the history of previous episodes is unknown or not taken into account: proceed with caution. Gac Sanit. 2017;31:227–34. 20. Shmueli G, Minka TP, Kadane JB, et al. A useful distribution for fitting discrete data: revival of the COM-Poisson. Appl Stat. 2005;54:127–42. 21. Schafer J, Olsen M. Multiple imputation for multivariate missing- data problems: a data Analyst’s perspective. Multivariate Behav Res. 1998;33:545–71. 22. Schafer J. Analysis of incomplete multivariate data. London: Chapman & Hall; 1997. 23. Rubin D. Multiple imputation for nonresponse in surveys. New York: John Wiley & Sons; 1987. 24. Hernández-Herrera G, Navarro A, Moriña D. Regression-based imputation of explanatory discrete missing data. arXiv preprint; 2020. Available at: https://​arxiv.​org/​abs/​2007.​15031 25. Moriña D, Navarro A. The R package survsim for the simulation of simple and complex survival data. J Stat Softw. 2014;59:1–20. 26. Navarro A, Moriña D, Reis R, et al. Hazard functions to describe patterns of new and recurrent sick leave episodes for different diagnoses. Scand J Work Environ Health. 2012;38:447–55. 27. Moriña D, Hernández-herrera G, Navarro A. miRecSurv Package: Prentice- Williams-Peterson models with multiple imputation of unknown number of previous episodes. R J. 2021;13(2). https://​journ​al.r-​proje​ct.​org/​archi​ve/​ 2021/​RJ-​2021-​082/​RJ-​2021-​082.​pdf. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub‑ lished maps and institutional affiliations. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions