Advertisement
Research Article| Volume 7, ISSUE 2, P253-260, June 2019

Download started.

Ok

Survival analysis in longitudinal studies for recurrent events: Applications and challenges

Open AccessPublished:February 23, 2019DOI:https://doi.org/10.1016/j.cegh.2019.01.013

Abstract

Background and objective

Recurrent event data analysis is most commonly used in biomedical research. However, the researchers dealing with recurrent events in survival analysis have ignored the assumption that the recurrent events are correlated. There are methods available that takes into account dependency between recurrent events. The main objective of this study was to demonstrate the recurrent event models using upper respiratory infection (URI) among Indian infants.

Methods

The frequency of URI among a birth cohort of 210 babies was evaluated monthly with nasopharyngeal swabbing. Data on 11 potential risk factors were noted. The extended Cox models, such as Andersen-Gill counting process (CP), Prentice-Williams-Peterson (PWP-CP), PWP–Gap time model, Marginal Model and Cox frailty model were applied. The better model was assessed based on Loglikehood test statistics.

Results

Of the four models, PWP-CP model had lower log likelihood value. The URI incidence rate was estimated to be 2.27 (95%CI: 1.70–3.03)times significantly higher in the month of July–October and 1.43 (95%CI: 1.19–1.71) times higher in the month of November–February as compared to March–June (p < 0.001). Nasopharyngeal colonization with S. pneumoniae was also another important risk factor (HR = 1.18, 95%CI:1.01–1.39, p = 0.03).

Conclusion

In the current study, PWP-CP model was found to be better model as compared to other models. Also biologically appropriate as subsequent events depend on the first event of URI. Hence, the choice of an appropriate method for analyzing the recurrent event data should not be decided only on statistical basis but also based on the research question.

Keywords

1. Introduction

Cox Proportional Hazard Model has been widely used by most researchers in the recent past due to its versatility and simplicity in nature for interpreting the results. This model is used in recurrent events such as repeated asthma attacks, episodes of upper respiratory infections, repeated myocardial infarctions, recurrent urinary tract infection among the renal transplant patients, etc. are very common in medical research. However, the researchers dealing with recurrent events in survival analysis have ignored the assumption that the recurrent events are correlated.
  • Cox D.R.
Regression models and life-tables.
  • Anker S.D.
  • McMurray J.J.
Time to move on from ‘time-to-first’: should all events be included in the analysis of clinical trials?.
  • Twisk J.
  • Smidt N.
  • de Vente W.
Applied analysis of recurrent events: a practical overview.
In such situation either they have used the latest event and the time related to that event as outcome or, they have assumed the recurrent events are independent and analysed data using survival analysis.
  • Purroy F.
  • Jiménez Caballero P.E.
  • Gorospe A.
  • et al.
Recurrent transient ischaemic attack and early risk of stroke: data from the PROMAPA study.
,
  • Gill D.P.
  • Zou G.Y.
  • Jones G.R.
  • Speechley M.
Comparison of regression models for the analysis of fall risk factors in older veterans.
If the correlations between the recurrent events are ignored, then the null hypothesis is mostly rejected, because the Cox model does not incorporate within subject correlation. However, methods have been developed that make use of all available data, while accounting for the lack of independence of recurrent events within subjects. The two popular approaches are namely, “Variance-corrected Cox based models” and “Frailty/random effects” models.
  • Box-Steffensmeier J.M.
  • De Boef S.
Repeated events survival models: the conditional frailty model.
  • Guo Z.
  • Gill T.M.
  • Allore H.G.
Modeling repeated time-to-event health conditions with discontinuous risk intervals: an example of a longitudinal study of functional disability among older persons.
  • Ullah S.
  • Gabbett T.J.
  • Finch C.F.
Statistical modelling for recurrent events: an application to sports injuries.
Variance-corrected models were developed to account for correlation by using robust (sandwich) standard errors. However, the theory behind frailty models is that some subjects are intrinsically more or less prone to experience events of interest than others; frailty can be considered as a random covariate in the model that corrects dependence among recurrent event times. Limitations of applying these variance-corrected models and frailty models include their complexity and difficulty in implementation. Generalised Estimating Equation (GEE) analysis is a variance corrected model that requires the specification of a working correlation matrix but are still inefficient as the information on time to event is ignored. Most statistical models have been developed for analysing the recurrent event data which largely focused on binary outcome having discontinuous risk intervals using GEE.
  • Guo Z.
  • Gill T.M.
  • Allore H.G.
Modeling repeated time-to-event health conditions with discontinuous risk intervals: an example of a longitudinal study of functional disability among older persons.
,
  • Amorim leila D.A.F
  • Cai J.
Modelling recurrent events: a tutorial for analysis in epidemiology.
Secondly, longitudinal data are analyzed using binary logistic regression ignoring the time, which is inappropriate. We intend to formulate a model which will consider the actual time of recurrence of each outcome. Though these concepts have been disseminated in statistical journals this is seldom practised in developing countries. Hence the aim of this paper is to summarize various methods for modelling recurrent event data. We would also show the differences in estimation and interpretation of recurrent event approaches, as well as to sensitize appropriate models, based on research objectives for the longitudinal study.

2. Materials and methods

2.1 Data

This study was conducted in K.V Kuppam rural development block, which belongs to the service area of RUHSA (Rural Unit for Health and Social Affairs) Christian Medical College, Vellore, India between February 2009 to August 2010. After taking an informed consent, a detailed socio-demographic history was obtained. Patient information was obtained from their parents who were interviewed at each visit regarding recurrent colds, allergic symptoms, overcrowding, family size, breast feeding, smoke exposure and day care attendance. At birth and at monthly scheduled visits, nasopharyngeal swabbing was performed with a calcium alginate swab stick. Then, the presence of upper respiratory infection was noted.
  • Rupa V.
  • Isaac R.
  • Manoharan A.
  • Jalagandeeswaran R.
  • Thenmozhi M.
Risk factors for upper respiratory infection in the first year of life in a birth cohort.

2.2 Standard Cox PH model

The standard Cox proportional hazard model for the survival data specifies the hazard of ith individual as,
λi(t) = λ0(t) exp{βxi}
(1)


Where (t) λ0(t) is an unspecified baseline hazard function and β is the vector of regression coefficients, xi is the vector of covariates of the ith subject.

2.3 Extended Cox model

The extended Cox models were used to model recurrent time-to-event outcomes within a subject comprehensively than the Cox model. The extended Cox models were: 1) the Andersen-Gill counting process (CP), 2) the Prentice-Williams-Peterson (PWP-CP/Total time),3) PWP – Gap time (PWP-GT) model, 4) Marginal (Wei, Lin and Weissfeld)Model and 5) Cox frailty model.

2.4 Andersen Gill model

Andersen Gill model assumes that the correlation between event times for a subject can be explained by the past events. AG model is suitable model when correlations among events for each individual are induced by measured covariates. The counting process style of data input is seen in AG model where each subject is represented as series of observation with recurrence time given as (t0, t1], (t1, t2] … (tm, last follow-up time] where, each recurrent event for the ith subject is assumed to follow a proportional hazard model is given as
λi(t) = λ0(t)exp{βk xi(t)}
(2)


Under this model, the risk of recurrent event for a subject follows the usual Cox PH assumption but the number of recurrence is not taken into account. Every subject risk intervals contribute to the risk set for every event, irrespective of the number of events for each individual (Fig. 1b).
  • Andersen P.K.
  • Gill R.D.
Cox's regression model for counting processes: a large sample study.
Fig. 1
Fig. 1Upper respiratory infection recurrent time to event data of birth cohort in the first year of life.

2.5 Prentice, William and Peterson model (PWP)

Another model for analyzing recurrent events is PWP model.
  • Prentice R.L.
  • Williams B.J.
  • Peterson A.V.
On the regression analysis of multivariate failure time data.
PWP CP model (total time) and PWP Gap time model. PWP CP model is similar to the AG-CP model but stratified by events. The baseline hazards vary from event to event, the hazard function for the kth event for the ith subject with the PH form is written as
λik(t) = λ0(t)exp{βk xi(t)}
(3)


The PWP - GT model describes an intensity process from the occurrence of an immediately preceding event, with the gap time defined as (ttk1). Both PWP approaches are conditional models as an individual is not considered in the riskset for the kth event until experiencing the (k−1)th event.
λik(t) = λ0(t-tk-1)exp{βk xi(t)}
(4)


λ0k(t) represents the event-specific baseline hazard for the kth event over time. AG model and both PWP Models are adjusted by estimating the sandwich type estimators and hence they are known as variance corrected models
  • Kleinbaum D.G.
  • Klein M.
Survival Analysis: A Self-Learning Text.
,
  • Clayton D.
Some approaches to the analysis of recurrent event data.
(Fig. 1d).

2.6 Marginal (Wei, Lin and Weissfeld) model

Wei, Lin and Weissfeld (1989) proposed a Cox-type model to analyse repeated events data.
  • Wei L.J.
  • Lin D.Y.
  • Weissfeld L.
Regression analysis of multivariate incomplete failure time data by modeling marginal distributions.
In most applications the analysis has been the “time from the study entry” scale, since all the time intervals start at zero
  • Hosmer Jr., David W.
  • Lemeshow S.
  • May S.
Applied Survival Analysis: Regression Modeling of Time to Event Data.
(Fig. 1c). The hazard function for the kth event for the ith individual is,
λik(t) = λ0(t)exp{βk xi(t)}
(5)


Unlike the AG model, this model allows a separate underlying hazard for each event. When an event is zero, it means that subject is no longer at risk after the last given event.
  • Cook R.J.
  • Lawless J.
The Statistical Analysis of Recurrent Events.
,
  • Therneau T.M.
  • Grambsch P.M.
Modeling Survival Data: Extending the Cox Model.

2.7 Cox frailty model

The frailty model is an extension of the Cox PH model, in which, the hazard function depends on an unmeasured random variable.
  • Therneau T.M.
  • Grambsch P.M.
Modeling Survival Data: Extending the Cox Model.
,
  • Hougaard P.
Analysis of Multivariate Survival Data.
The term ‘frailty’ means that each subject has his/her own disposition to failure, in additional to any effects that will be quantified using regression. Hazard function λij(t) for the recurrent time of the kth event in the ith subject (j = 1,2,…ki; i = 1,2,…n) conditional on the frailty Zi, follows the PH form and its given by:
λik(t)=λ0k(t)Zi{exp{xi(t)βk},t>0
(6)


Where, λ0k(t) is the common baseline hazard function, Xi is a vector of observable covariates and β is a vector of unknown regression coefficients. Frailty Zi is the unobserved (random) common risk factors shared by all subjects in cluster ‘i’ and is assumed to be i.i.d random variable with unit mean and unknown variance θ.
  • Hougaard P.
Analysis of Multivariate Survival Data.
,
  • Duchateau L.
  • Janssen P.
The Frailty Model.
The Frailty effects occur when the observed sources of variation in the observed or unobserved explanatory variables fail to account for the true difference in risk. That is, when there are other important but omitted variables presented, the effect of omitted variable can be captured by frailty.

3. Results

Number of recurrence experienced by infants ranged between zero to ten during the follow-up period. Seventeen infants (8.1%) out of 210 infants did not return to the study area after birth. The upper respiratory infection recurred at least once in 193 subjects and highest recurrence events (9 and 10 times) were observed in 7 patients.
Table 1 shows a summary of follow-up times and number of patients with and without URI event for the consecutive recurrent events. The median follow-up time to the first URI event were 98 days and starts increasing for the higher consecutive recurrent event.
Table 1Summary of time between consecutive URI recurrent Events.
Follow-up time (in days)No of patient with URI
MinMaxMedianEventCensoredTotal
1st recurrence23399818518203
2nd recurrence4938219716213175
3rd recurrence7739324313512147
4th recurrence10540527610219121
5th recurrence152424311691887
6th recurrence17539832144953
7th recurrence20338633832234
8th recurrence23138433611516
9th recurrence287363347617
10th recurrence315377346112
A total of 163 infants (77.6%) had 6–13 visits where as 30 infants (14.3%) made <5 visits. The median number of visits for these 193 infants was 9 visits. In thousand days of life, 845 records from 747 upper respiratory patients were followed–up during the study period and three infants died during the period of the study. The socio-demographic data were obtained at birth and child characteristics were presented in Table 2a and Table 2b. More parents resided in tiled/pucca houses than thatched houses (66.7%), and were labourers/unemployed. The majority of parents (61.9%) had at least high school education (84.8% fathers and 86.7% mothers). Majority of household had (98.1%) < 3 under-five children, 76.9% of the households had more than 4 family members.
Table 2aSocio demographic and baseline characteristics.
VariablesBaseline (n = 210)
n%
Sex
 Male121 (57.6)
 Female89 (42.4)
Type of House
 Thatched70 (33.3)
 Tiled/Terraced/Group House140 (66.7)
Parental Occupation
 Nil/Laborer130 (61.9)
 Petty Business/Professional/Others80 (38.1)
Father's Education
 Illiterate/Primary32 (15.2)
 High/Higher Secondary and above178 (84.8)
Mother's Education
 Illiterate/Primary28 (13.3)
 High/Higher Secondary and above182 (86.7)
Birth weight (Grams)
 ≤ 250076 (36.2)
 > 2500134 (63.8)
Smoke
 Yes15 (7.1)
 No195 (92.9)
No. of members in the house
 ≤448 (23.1)
 > 4160 (76.9)
Firev
 Yes85 (40.5)
 No125 (59.5)
Water
 Bore well124 (59.0)
 River/Open Well86 (41.0)
Nasophryngeal Swab Report
 Positive8 (3.8)
 Negative201 (96.2)
Table 2bSocio Demographic and baseline characteristics by URI recurrent events.
VariableEvent 1Event 2Event 3Event 4Event 5
n (%)n (%)n (%)n (%)n (%)
Sex
Female85 (41.9)73 (41.7)62 (42.2)50 (41.3)93 (46.7)
Male118 (58.1)102 (58.3)85 (57.8)71 (58.7)106 (53.3)
Type of house
Tiled/Terraced/Grouped Houses136 (67.0)120 (68.6)104 (70.7)82 (67.8)134 (67.3)
Thatched67 (33.0)55 (31.4)43 (29.3)39 (32.2)65 (32.7)
Occupation
Farmer/Bigbusiness/Petty business75 (36.9)63 (36.0)52 (35.4)39 (32.2)54 (27.1)
Nil/Labourer128 (63.1)112 (64.0)95 (64.6)82 (67.8)145 (72.9)
Father Education
High school/Secondary and above173 (85.2)150 (85.7)126 (85.7)102 (84.3)164 (82.4)
Illiterate/Primary30 (14.8)25 (14.3)21 (14.3)19 (15.7)35 (17.6)
Mother Education
High school/Secondary and above186 (91.6)161 (92.0)135 (91.8)111 (91.7)180 (90.5)
Illiterate/Primary17 (8.4)14 (8.0)12 (8.2)10 (8.3)19 (9.5)
Birth weight
>2.5 kg130 (64.0)109 (62.3)90 (61.2)71 (58.7)116 (58.3)
≤2.5 kg73 (36.0)66 (37.7)57 (38.8)50 (41.3)83 (41.7)
Smoking
No189 (93.1)163 (93.1)138 (93.9)112 (92.6)180 (90.5)
Yes14 (6.9)12 (6.9)9 (6.1)9 (7.4)19 (9.5)
Mem5
≤4188 (93.5)173 (98.9)145 (98.6)120 (99.2)196 (98.5)
>413 (6.5)2 (1.1)2 (1.4)1 (0.8)3 (1.5)
Fire
No123 (60.6)107 (61.1)91 (61.9)75 (62.0)101 (50.8)
Yes80 (39.4)68 (38.9)56 (38.1)46 (38.0)98 (49.2)
Water
Bore well121 (59.6)106 (60.6)95 (64.6)78 (64.5)123 (61.8)
Open well/River82 (40.4)69 (39.4)52 (35.4)43 (35.5)76 (38.2)
Swab
Negative139 (68.8)102 (58.3)80 (54.4)83 (68.6)166 (83.4)
Positive63 (31.2)73 (41.7)67 (45.6)38 (31.4)33 (16.6)
Season(Months)
March to June24 (11.8)28 (16.0)42 (28.6)56 (46.3)128 (64.3)
July to October118 (58.1)56 (32.0)23 (15.6)8 (6.6)21 (10.6)
November to February61 (30.0)91 (52.0)82 (55.8)57 (47.1)50 (25.1)
Fig. 1a: We have considered 10 children to demonstrate the risk intervals using URI data. Among those 10 children, seven had at least three events. Remaining three of the children was censored at 365 days; Subject 4 and 6 had largest number of events (6 events). Subject 3 had only two event at times 154 and 336 days.
Variance corrected models and Frailty model results are presented in Table 3, the ‘seasonal’ variable was the only consistently significant risk factor for an URI. Cox PH model had lower log-likelihood values. However, these lower likelihood values do not represent ‘good fit’ because it does not consider the subsequent events within each child.
Table 3Risk factors for upper respiratory infection (URI) recurrent event data using Variance-Corrected model and frailty model in the first year of life.
VariableModel 1 (AG Model)Model 2 (PWP Total Time Model)Model 3 (PWP Gap time Model)Model 4 (Marginal Model)Model 5 (Cox Frailty Model)
HR95% CIP ValueHR95% CIP ValueHR95% CIP ValueHR95% CIP ValueHR95% CIP Value
Season
March–June1.001.001.001.001.00
July to October2.271.70–3.03<0.0012.601.14–5.940.0232.221.66–3.00<0.0011.581.05–2.370.0272.301.84–2.86<0.001
November to February1.431.19–1.71<0.0011.501.12–2.020.0071.371.11–1.690.0032.501.11–1.69<0.0011.441.17–1.76<0.001
Sex
Male0.950.83–1.090.4920.930.86–1.010.0980.930.81–1.070.3081.090.81–1.070.5040.950.82–1.110.530
Female1.001.001.001.001.00
Swap
Positive1.231.07–1.420.0031.181.08–1.390.0391.130.97–1.310.1271.491.14–1.950.0031.221.03–1.440.014
Negative1.001.001.001.001.00
Mem5
≤41.010.48–2.160.9661.020.49–2.110.9540.910.42–1.970.8110.970.42–1.970.8491.660.84–3.310.960
>41.001.001.001.001.00
Smoking
Yes1.040.86–1.250.6931.040.83–1.300.7371.050.84–1.310.6580.990.84–1.310.9631.040.79–1.370.790
No1.001.001.001.001.00
Water
Open well/River0.970.84–1.140.7480.960.80–1.160.6911.010.87–1.160.9471.080.87–1.160.5620.970.83–1.120.750
Borewell1.001.001.001.001.00
Fire
Yes1.140.98–1.330.0871.151.04–1.270.0041.181.02–1.350.0221.221.02–1.350.1311.140.98–1.340.920
No1.001.001.001.001.00
Father Education
Illterate/Primary1.070.89–1.280.4831.070.92–1.250.3721.020.86–1.230.7901.010.86–1.230.9601.070.87–1.320.540
High school & above1.001.001.001.001.00
Mother Education
Illterate/Primary0.970.74–1.280.8500.890.71–1.120.3300.990.76–1.310.9770.780.96–1.270.2230.980.75–1.270.850
High school & above1.001.001.001.001.00
Birth_weight
≤2.5 kg1.110.96–1.280.1511.141.04–1.270.0091.110.96–1.270.1451.301.08–1.270.0351.110.95–1.300.190
>2.5 kg1.001.001.001.001.00
Parent's Occupation
Nil/Labourer1.110.95–1.300.1771.110.97–1.270.1221.080.91–1.270.3670.920.86–1.160.5761.100.94–1.290.190
Professional & Others1.001.001.001.001.00
Type of House
Pucca/Kacha0.990.85–1.160.9631.010.82–1.230.9430.990.86–1.160.9750.980.86–1.160.9081.020.87–1.200.960
Tiled/Terraced/Grouped house1.001.001.001.001.00
Frailty Variance0.000.870
Log likelihood−3712.08−2578.19−2918.27−5627.17−3706.00
R Square0.0830.0920.0720.1290.083

3.1 Variance corrected model comparison

Children were most predisposed to upper respiratory infection in the July–October months and November–February months, which is statistically significant with slight differences in their parameter estimates in all the models. Using AG model, the upper respiratory infection was estimated to be 2.27 (95%CI: 1.70–3.03) significantly higher in July–October months and 1.43 (95%CI: 1.19–1.71) in November–February months as compared to March–June Months (p < 0.001). The PWP gap time model showed HR of 2.22 (95%CI: 1.66–3.00) for the July–October months and 1.37 (95%CI: 1.11–1.69) for November–February months respectively compared to March–June Months, which is statistically significant (p < 0.01). The WLW model for total time up to the 10th URI recurrence since the study entry yielded a HR of 1.58 (95%CI: 1.05–2.37, P value = 0.027) in July–October months and 2.50 (95%CI: 1.87–3.32) in November–February months as compared to March–June Months (p < 0.001). Nasopharyngeal colonization with S. pneumoniae was another important risk factor which was significant in all recurrent event models except PWP Gap time model. (AG model: HR = 1.23, 95% CI = 1.07–1.42, p value = 0.003; PWP total time model: HR = 1.18, 95% CI: 1.01–1.39, p value = 0.03; Marginal model: HR = 1.49, 95% CI: 1.14–1.95, p value = 0.003). Birth weight of an infant <2.5 kgs had a risk of 1.14 (95% CI: 1.04–1.27) times of URI infection as compared to normal birth weight infant (Table 3).

3.2 AG model and frailty model comparison

The parameter estimates obtained from the frailty model with counting process time scale and AG models were almost same without frailty term. In other words, when the frailty model, with a variance almost close to zero (θ = 0) would indicate that the frailty component does not contribute to the model. Based on the AG model, children with nasopharyngeal colonization with S. pneumoniae positive had high risk of recurring URI, which was 23% higher as compared to children without nasopharyngeal colonization by S. pneumoniae.Children were most susceptible to URI in July–October months (HR: 1.43, 95%CI: 1.70–3.03) and November–February months (HR: 2.27, 95%CI: 1.19–1.71) as compared to March–June months, which is statistically significant (p < 0.001). The cumulative hazard plot in Fig. 2a showed that both AG model and frailty model have estimated same cumulative hazard in the study and it clearly showed that if frailty variance is not significant. The frailty variance θ was estimated to be 0 and 0.153 for counting process time scale and gap respectively. Fig. 2b shows that both models have the different estimated cumulative hazard over a time. Recurrent event data structure and how to organize the data for each recurrent event models, R code is given in the Appendix.
Fig. 2
Fig. 2Cumulative hazard plot for upper respiratory infection recurrence over a time of follow-up for AG model and frailty models.

4. Discussion

Upper respiratory infections are the most common cause of morbidity in the first year of life among the Indian infants. We described the relevant methods, its importance, how to fit and interpret the results for different methods. Cox PH model does not examine the effect of the risk factors on the number of recurrence over the follow up time.
  • Pandeya N.
  • Purdie D.M.
  • Green A.
  • Williams G.
Repeated occurrence of basal cell carcinoma of the skin and multifailure survival analysis: follow-up data from the Nambour Skin Cancer Prevention Trial.
,
  • Dancourt V.
  • Quantin C.
  • Abrahamowicz M.
  • Binquet C.
  • Alioum A.
  • Faivre J.
Modeling recurrence in colorectal cancer.
Many researchers continuously used logistic regression, GEE, Poisson and negative binomial approaches for estimating the risk factors for recurrent events.
  • Guo Z.
  • Gill T.M.
  • Allore H.G.
Modeling repeated time-to-event health conditions with discontinuous risk intervals: an example of a longitudinal study of functional disability among older persons.
,
  • Pandeya N.
  • Purdie D.M.
  • Green A.
  • Williams G.
Repeated occurrence of basal cell carcinoma of the skin and multifailure survival analysis: follow-up data from the Nambour Skin Cancer Prevention Trial.
However, they considered the total number of events per fixed period of time interval, ignoring the actual time to event concept between repeated occurrences.
  • Rupa V.
  • Isaac R.
  • Manoharan A.
  • Jalagandeeswaran R.
  • Thenmozhi M.
Risk factors for upper respiratory infection in the first year of life in a birth cohort.
,
  • Pandeya N.
  • Purdie D.M.
  • Green A.
  • Williams G.
Repeated occurrence of basal cell carcinoma of the skin and multifailure survival analysis: follow-up data from the Nambour Skin Cancer Prevention Trial.
Observations from the same child are expected to be correlated and hence Cox PH model is not suitable method to account of the extra variability of the recurrent event data. So, variance corrected models and frailty model were used.
Several methods have been proposed to account for intra-individual correlation that rises from recurrent events setting in survival analysis. The biological reason for the infection/disease is essential when selecting a model for the recurrent events for example if it is possible that after experiencing the first URI infection, the risk to the next infections may increase. If AG model is reasonable to assume that the risk of the repeated infections remains constant, irrespective of the number of previous infections, then the AG model is recommended.
  • Therneau T.M.
  • Grambsch P.M.
Modeling Survival Data: Extending the Cox Model.
AG Model provides more powerful inference for a covariate effect than the Cox model. A robust sandwich estimate is used to evaluate the standard errors.
  • Andersen P.K.
  • Gill R.D.
Cox's regression model for counting processes: a large sample study.
The PWP models are reasonable to assume that the child can only be at risk for a given event after he/she experiences the previous event. The PWP model means that the underlying hazard function is assumed to be the same from event to event. Also that, when infections increase with subsequent recurrence, the PWP model may be more appropriate than the AG model.
The WLW model overestimates the covariate effects due to the fact that every child has as many records as the maximum number of the event occurred in our data.
  • Rupa V.
  • Isaac R.
  • Manoharan A.
  • Jalagandeeswaran R.
  • Thenmozhi M.
Risk factors for upper respiratory infection in the first year of life in a birth cohort.
The variance corrected models handle correlation between the recurrent events occurring in the same patient by only correcting for the variance.
  • Kelly P.J.
A review of software packages for analyzing correlated survival data.
,
  • Lee E.
  • Wei L.
  • Amato D.
Cox-type regression analysis for large numbers of small groups of correlated failure time observations.
Random effect/Frailty models leads to a person specific interpretation of the estimates which is similar to that of mixed models in longitudinal data in order to account for the dependency between the recurrent event and unobserved heterogeneity among patients properly as this cannot be explained by covariates alone.
  • Guo G.
Use of sibling data to estimate family mortality effects in Guatemala.
Frailty model gives consistent estimates based on the distribution of the number of events and sample size. A small frailty variance implies very low correlation between the event times.
  • Chen L.M.
  • Ibrahim J.G.
  • Chu H.
Sample size determination in shared frailty models for multivariate time-to-event data.
Frailty model is mainly applied for a moderate to large number of events but even for a small number of events, it is quite adequate.
  • Therneau T.M.
  • Grambsch P.M.
Modeling Survival Data: Extending the Cox Model.
,
  • Wienke A.
Frailty Models in Survival Analysis.
In this study the within subject correlation is very low (0.07). However, this need not necessarily be the case always, when dealing with data. Therefore use of the entire model is based on the concept of the problems.
There is strong evidence in the literature that if frailty is present but ignored, then the covariate effects will be underestimated.
  • Finkelstein D.M.
  • Schoenfeld D.A.
  • Stamenovic E.
Analysis of multiple failure time data from an aids clinical trial.
,
  • Pickles A.
  • Crouchley R.
A comparison of frailty models for multivariate survival data.
If the common baseline hazard between each event time is not appropriate for repeated event data and also, when a robust variance to any of these models does not adequately account for within-subject correlation, then it has been suggested to apply the frailty model which is also a similar finding from the present study.
  • Kelly P.J.
A review of software packages for analyzing correlated survival data.
However, if the primary interest of investigation is a measurement of dependence of within subject, then the frailty model is more adequate.
  • Hougaard P.
Frailty models for survival data.
The difference between WLW model and frailty models is driven from the fact that, the frailty model is more naturally related to the fundamental performance of recurrences, while the unconditional WLW model does not provide understanding of the interrelationship among recurrences. However, it has been suggested that, if the frailty distribution is correctly defined, then the frailty model is expected to be more efficient than the WLW model.
  • Lin D.Y.
Cox regression analysis of multivariate failure time data: the marginal approach.
In summary, the choice of appropriate model for analyzing recurrent event data will be influenced by many factors, such as number of events, relationship among subsequence events, within subject correlation and varying covariates and the sample size. AG model is appropriate, when the assumption of a common underlying hazard over recurrent event observation is reasonable and when only interested in the overall rate of recurrences. When the dependence from the past event is strong and consistent, then the PWP model is appropriate. However, when the distribution of event per subject is small or prediction of time to the next event is of interest, the PWP gap time model is the appropriate method. However, for the present study, each episode of URI within a child is biologically related though the estimated correlation was very small. This was because the risk of infection to the same sero group/sero type was less in subsequent events within a child; thereby the estimated correlation coefficient was close to 0. However, if the researcher was interested to study the measurement of the dependence between recurrent event times within the subject, the frailty model would be appropriate.

5. Conclusion

The present study finding suggests that the choice of an appropriate method for analyzing the recurrent event data should not mainly depend on statistical basis such as model with low likelihood values; rather the selection should also be based on the research question, a thorough clinical knowledge on the events of interest followed by organization of the data. Thus the PWP-CP model fit the data appropriately while the biological process also suggested the same model.

Conflicts of interest

None.

Appendix. Data structures for modelling recurrent event data

Organisation of the dataset is a more complication than the usual discontinuous risk intervals. Each subject is represented by several rows of data dependent on number of events child had, with time organized as intervals that represent (entry time, 1st event], (1st event, 2nd event],……, (kth event, end time]. A key difference for fitting recurrent event models is the creation of appropriate datasets. To show up the important features of data structure we present some information about the datasets that we used in the present paper.
Let's consider the example of first five children details from the Upper Respiratory Infection data are presented below:
Tabled 1
Study IDStartStopURI StatusGapSexSwabMonthsURI count
1022612261031
12262821561032
12823101281133
13103380281014
20841841011
2841271431022
21271471201023
21471681211024
216832211541035
3013211322111
31322021702122
32022301282123
32303001702034
33003281282035
33283561282135
4015411541011
415433611821032
50351351011
53527612411032
52763021261033
53023441421134
A pair of variable (start, stop) is used to define the time interval of the URI. The start time is generally equal to 0 for the 1st URI and equals to the last recurrence time for further URI. The stop time is a recurrent URI time (URI status = 1) or censored time (URI-status = 0). The study ID variable identifies the child's. 1st child study ID = 1 from 226 to 282, 310 and 338 – with start time equal to 0 and stop time equal to follow up time, while child have four rows (study ID = 1 and 5). Child with no censoring in the end of the follow-up but the child five have 3 event and end of the visit he/she became censoring. Child study ID 4 had an event time at 154 days and second event time at 336 days. For five child data corresponding covariates are presented in the following column in the above table. This structure of the data can used to fit AG model and Frailty Model. In the PWP total time model, Gap time model addition information we added as URI counts based on the number of URI occurrence in the study duration, which is going to be used for stratification. The PWP Gap time model and frailty gap time model the time is defined as stop minus start time and the data structure as same.
Marginal approach focuses on total survival time from study entry until the occurrence of a specific (e.g., Kth) event. Suggested when recurrent events are viewed to be of different types also. Each subject is considered to be at risk for all failures that might occur, regardless of no: of events a subject actually experienced. For example in our study, every child to be at risk as the maximum number of recurrent events occurred in the study (k = 10) event if a child has one recurrence. i.e, every child has 10 observations, one in each stratum. In this data the URI event indicator, which is going to be used for stratification. Strata will correspond to the number of URI. Risk set determined from time since study entry. Marginal model is stratified model. The below data structure can be used to fit the marginal models.
Tabled 1
Study IDStartStopURI StatusSexSwabMonthsURI count
1022611031
1028211032
1031011133
1033801014
1033811015
1033811016
1033811017
1033811018
1033811019
10338110110
R Code for the entire Model:
The library survival in R allows all recurrent event models, which is discussed in this paper.
  • library (survival)
  • library (foreign)
  • uri < -read.spss ("file location", use.value.labels = TRUE, to.data.frame = FALSE)
AG Model:
  • AG_Model < -coxph (Surv (Start, Stop,URI_status)∼Mon_R + Sex_r + Swap_r.
  • +smk_r + water_r + fire_r + bwt_r + Pocc_r2+Toh_r + cluster (StudyID), data = uri)
  • summary (AG_Model)
Stratification Models: For the below models are stratified Models, the argument strata (URI_Count) identifies stratification variable to obtain their estimates. Estimates are obtained for event-specific effects for each covariates.
1. PWP-Total time Model:
  • PWP_TT < -coxph (Surv (Start, Stop,URI_status)∼Mon_R + Sex_r + Swap_r.
  • +smk_r + water_r + fire_r + bwt_r + Pocc_r2+Toh_r + cluster (StudyID)+Strata (URI_Count),data = uri)
  • summary (PWP_TT)
2. PWP-Gap time Model:
  • PWP_GT < -coxph (Surv (Stop-Start,URI_status)∼Mon_R + Sex_r+
  • Swap_r + smk_r + water_r + fire_r + bwt_r + Pocc_r2+Toh_r + cluster.
  • (StudyID)+Strata (URI_Count),data = uri)
  • summary (PWP_GT)
3. Marginal Model:
  • Marginal < -coxph (Surv (Start, Stop,URI_status)∼Mon_R + Sex_r+
  • Swap_r + smk_r + water_r + fire_r + bwt_r + Pocc_r2+Toh_r + cluster (StudyID)+Strata (URI_Count),data = uri)
  • Summary (Marginal)
Frailty Model: By default gamma distribution is associated to the random effect for the frailty model in R software. However, we can specify the distributions such as gamma and Gaussian. Other way frailty.gamma (Study_ID) and frailty.gaussian (Study_ID) instead of frailty (id,dist = ”gamma”)
  • Frailty < -coxph (Surv (Start, Stop,URI_status)∼Mon_R + Sex_r + Swap_r.
  • +mem5_r + smk_r + water_r + fire_r + Fathedu_r + MothEdu_r + bwt_r + Pocc_r2+Toh_r + frailty (StudyID, dist = "gamma"), data = uri)
  • summary (Frailty)

References

    • Cox D.R.
    Regression models and life-tables.
    J R Stat Soc Ser B Methodol. 1972; 34: 187-220
    • Anker S.D.
    • McMurray J.J.
    Time to move on from ‘time-to-first’: should all events be included in the analysis of clinical trials?.
    Eur Heart J. 2012; 33: 2764-2765
    • Twisk J.
    • Smidt N.
    • de Vente W.
    Applied analysis of recurrent events: a practical overview.
    J Epidemiol Community Health. 2005; 59: 706-710
    • Purroy F.
    • Jiménez Caballero P.E.
    • Gorospe A.
    • et al.
    Recurrent transient ischaemic attack and early risk of stroke: data from the PROMAPA study.
    J Neurol Neurosurg Psychiatry. 2013; 84: 596-603
    • Gill D.P.
    • Zou G.Y.
    • Jones G.R.
    • Speechley M.
    Comparison of regression models for the analysis of fall risk factors in older veterans.
    Ann Epidemiol. 2009; 19: 523-530
    • Box-Steffensmeier J.M.
    • De Boef S.
    Repeated events survival models: the conditional frailty model.
    Stat Med. 2006; 25: 3518-3533
    • Guo Z.
    • Gill T.M.
    • Allore H.G.
    Modeling repeated time-to-event health conditions with discontinuous risk intervals: an example of a longitudinal study of functional disability among older persons.
    Methods Inf Med. 2008; 47: 107-116
    • Ullah S.
    • Gabbett T.J.
    • Finch C.F.
    Statistical modelling for recurrent events: an application to sports injuries.
    Br J Sports Med. 2014; 48: 1287-1293
    • Amorim leila D.A.F
    • Cai J.
    Modelling recurrent events: a tutorial for analysis in epidemiology.
    Int J Epidemiol. 2015; 44: 324-333
    • Rupa V.
    • Isaac R.
    • Manoharan A.
    • Jalagandeeswaran R.
    • Thenmozhi M.
    Risk factors for upper respiratory infection in the first year of life in a birth cohort.
    Int J Pediatr Otorhinolaryngol. 2012; 76: 1835-1839
    • Andersen P.K.
    • Gill R.D.
    Cox's regression model for counting processes: a large sample study.
    Ann Stat. 1982; 10: 1100-1120
    • Prentice R.L.
    • Williams B.J.
    • Peterson A.V.
    On the regression analysis of multivariate failure time data.
    Biometrika. 1981; 68: 373-379
    • Kleinbaum D.G.
    • Klein M.
    Survival Analysis: A Self-Learning Text.
    Springer Science & Business Media, 2005: 616
    • Clayton D.
    Some approaches to the analysis of recurrent event data.
    Stat Methods Med Res. 1994; 3: 244-262
    • Wei L.J.
    • Lin D.Y.
    • Weissfeld L.
    Regression analysis of multivariate incomplete failure time data by modeling marginal distributions.
    J Am Stat Assoc. 1989; 84: 1065-1073
    • Hosmer Jr., David W.
    • Lemeshow S.
    • May S.
    Applied Survival Analysis: Regression Modeling of Time to Event Data.
    2 edition. Wiley-Interscience, Hoboken, N.J2008 (416 p)
    • Cook R.J.
    • Lawless J.
    The Statistical Analysis of Recurrent Events.
    Springer Science & Business Media, 2007 (415 p)
    • Therneau T.M.
    • Grambsch P.M.
    Modeling Survival Data: Extending the Cox Model.
    Springer Science & Business Media, 2000: 372
    • Hougaard P.
    Analysis of Multivariate Survival Data.
    Springer, New York2001 (542 p)
    • Duchateau L.
    • Janssen P.
    The Frailty Model.
    2008 edition. Springer, New York2008 (316 p)
    • Pandeya N.
    • Purdie D.M.
    • Green A.
    • Williams G.
    Repeated occurrence of basal cell carcinoma of the skin and multifailure survival analysis: follow-up data from the Nambour Skin Cancer Prevention Trial.
    Am J Epidemiol. 2005; 161: 748-754
    • Dancourt V.
    • Quantin C.
    • Abrahamowicz M.
    • Binquet C.
    • Alioum A.
    • Faivre J.
    Modeling recurrence in colorectal cancer.
    J Clin Epidemiol. 2004; 57: 243-251
    • Kelly P.J.
    A review of software packages for analyzing correlated survival data.
    Am Stat. 2004; 58: 337-342
    • Lee E.
    • Wei L.
    • Amato D.
    Cox-type regression analysis for large numbers of small groups of correlated failure time observations.
    in: Klein J. Goel P. Survival Analysis: State of the Art. Kluwer Academic Publishers, Dordrecht1992: 237-247
    • Guo G.
    Use of sibling data to estimate family mortality effects in Guatemala.
    Demography. 1993; 30: 15-32
    • Chen L.M.
    • Ibrahim J.G.
    • Chu H.
    Sample size determination in shared frailty models for multivariate time-to-event data.
    J Biopharm Stat. 2014; 24: 908-923
    • Wienke A.
    Frailty Models in Survival Analysis.
    CRC Press, 2010 (322 p.)
    • Finkelstein D.M.
    • Schoenfeld D.A.
    • Stamenovic E.
    Analysis of multiple failure time data from an aids clinical trial.
    Stat Med. 1997; 16: 951-961
    • Pickles A.
    • Crouchley R.
    A comparison of frailty models for multivariate survival data.
    Stat Med. 1995; 14: 1447-1461
    • Hougaard P.
    Frailty models for survival data.
    Lifetime Data Anal. 1995; 1: 255-273
    • Lin D.Y.
    Cox regression analysis of multivariate failure time data: the marginal approach.
    Stat Med. 1994; 13: 2233-2247