The Poisson-exponential model for recurrent event data: an application to bowel motility data

This paper presents a new parametric model for recurrent events, in which the time of each recurrence is associated to one or multiple latent causes and no information is provided about the responsible cause for the event. This model is characterized by a rate function and it is based on the Poisson-exponential distribution, namely the distribution of the maximum among a random number (truncated Poisson distributed) of exponential times. The time of each recurrence is then given by the maximum lifetime value among all latent causes. Inference is based on a maximum likelihood approach. A simulation study is performed in order to observe the frequentist properties of the estimation procedure for small and moderate sample sizes. We also investigated likelihood-based tests procedures. A real example from a gastroenterology study concerning small bowel motility during fasting state is used to illustrate the methodology. Finally, we apply the proposed model to a real data set and compare it with the classical Homogeneous Poisson model, which is a particular case.


Introduction
Recurrent event data are usually observed in longitudinal studies involving multiple subjects. This kind of data set arises in several areas, such as biomedical studies with reappearance of cancerous tumors, criminology studies with relapse of an offender in a crime, industrial studies with recurrent failures of an equipment, demographic studies with repeated migrations, among others.
In the literature, two types of time scale are often used to analyze recurrent event data, namely the event time of recurrence [13,15]  The remainder of the paper is organized as follows. In Section 2, we present the proposed model and its properties. The likelihood function for statistical inference, model assessment and hypothesis testing are detailed in Section 3. Section 4 contains the results of a simulation study on the behavior of the maximum likelihood estimates when a possibly small or moderate number of recurrent events per subject is observed. Likelihood-based tests procedures for testing the adequacy of the proposed model and their particular case are also discussed in this section. A real data analysis on a gastroenterology study concerning small bowel motility during fasting state for 19 individuals are presented in Section 5. Finally, Section 6 concludes the paper with final comments.

Model formulation
The proposed model in this paper is described in the context of recurrent events and it is founded on a gap time model formulation. For explicit representation of our model, we denote t as a calendar time (referred to as the event time) and w as another time of interest, such as a gap (interevent) time. Motivated by the idea of Zhao and Zhou [41], we consider a model in which the general form of the rate function is given by where λ(w|t) is the rate function for the recurrence process up to time t + w and λ 0 (t) is a deterministic function describing the general behavior of a subject over time.
Hereafter, we suppose that n independent subjects are submitted to experience an initial event and the recurrences of the same event are registered along the study period. Subjects are indexed by i, i = 1, 2, . . . , n, and each subject's recurrences are indexed by j, j ≥ 1. For subject i, let (W ij ) j≥1 be the sequence of gap times between successive occurrences of the event. Let 0 < T i1 < T i2 < · · · represent the times at which event occurs, where T ij = W i1 + . . . + W ij is the occurrence time of the jth event of subject i.
Specifically, conditional on T i,j−1 = t i,j−1 , we consider that the baseline rate function has a Poisson-exponential form and then the rate function of recurrence process N i (t i,j−1 + w) for subject i is defined by where θ, α > 0 and w > 0. The rate function increases with recurrence times but stabilizes at α. Besides, when θ approaches zero the rate function will be close to α, which is the upper bound of rate function, and hence the model reduces to the classical homogeneous Poisson process (HPP) with constant rate (i.e. the gap times W ij between successive events are independent and identically distributed exponential random variables with mean α −1 ). Figure 1 illustrates rate function shapes for some selected values of θ and fixed α. The PEre parameters have a direct interpretation in terms of complementary risks. The θ/(1 − e −θ ) represents the mean of the number of complementary risks, while α denotes the failure rate of the distribution of time-to-event due to individual complementary risks. It can be noted that the expected number of complementary risks tends to 1 as the parameter θ tends to zero. We now consider the corresponding cumulative rate function over the interval (t i,j−1 , t i,j−1 + w] given by of subject i is assumed to be a nonhomogeneous Poisson process with rate function (2), and together with (3), we can write So for a process that is independent of the history of events prior to t i,j−1 , the conditional distributions of gap times W ij are given by Cook and Lawless [8] Pr Let S(w|t i,j−1 ) denote the survival function for subject i. Considering that the rate function expressed by Equation (2) is deterministic and integrable over (t i,j−1 , t i,j−1 + w] and (t i,j−1 , w) is continuous, the survival function of W ij , given T i,j−1 = t i,j−1 , for each subject is given by and thus the each subject's gap times are not in general statistically independent. Moreover, it follows from Equation (5) that conditional on T i,j−1 = t i,j−1 , the jth gap time W ij of subject i has rate function (2) and density function given by The proof that the function in Equation (7) is a probability density function is trivial and a known result, then it is omitted.
Hereafter, we present a general expression for the raw moments of the jth gap time W ij of subject i, given by Equation (8), which follow directly from the moments of Cancho et al. [4] considering Equation (7). Then for r ∈ N, the raw moments are given by is the generalized hypergeometric function (see [25, p. 405] ). Hence, considering Equation (8), the mean and variance of the jth gap time W ij of subject i are given, respectively, by Finally, a method for data generation is presented below. The numerical study of the inferential method properties requires the simulation of the gap times, driven by the corresponding rate function (2), and of the times of events.
The gap times (w ij ) j≥1 are simulated using the iterative inverse-transform algorithm [30,Chapter 2]. The idea of the inverse-transform method is that for a uniform random variable U on the interval (0, 1), and for any continuous distribution function F, the random variable X, defined as X = F −1 (U), has distribution function F.
Consider a random variable W ij , the gap time between the (j − 1)st and the jth recurrent event times of a subject i. The conditional distribution function F of W ij may be used for inversion. Specifically, Then the gap times and times of events are generated through the following steps: (1) Draw u ij uniform(0, 1); (2) Use the Equations (9) and (2) and obtain the distribution of gap time W ij (3) Solve Equation (10) for w, and obtain w ij , a realization of the random variable W ij . Then the general expression for generated gap times is easily obtained: ( 1 1 ) (4) Get the times of events by doing t ij = t i,j−1 + w ij , for i = 1, . . . , n and j ≥ 1, with t i0 = 0.

Statistical inference
The parameters (α, θ) of the model are estimated using the maximum likelihood method, described in the current section. We assume that the recurrences of subjects are driven by a rate function as defined in Section 2. Let i be a subject observed over the interval time (0, τ i ], i = 1, . . . , n, where t = 0 corresponds to the start of the recurrence process. More precisely, for subject i, if m i recurrences are observed at times . . , m i and t i0 = 0. We assume that the stopping time τ i coincides with the occurrence time of the m i failure, that is, Assuming that each subject's gap times, W i1 , . . . , W im i , are associated with a censoring variable, and the censoring mechanism is noninformative, the contribution to the likelihood from each individual with respect to all occurrences and all intervals is where c ij is a censoring indicator, which is equal to 0 if the gap time is censored or 1 if the gap time is completely observed. Then, from Equation (6) we obtain the log-likelihood function as The maximum likelihood estimates (MLEs) of the parameters can be obtained by direct maximization of the log-likelihood function (13) using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization procedure of the R system [29].
Considering the usual large-sample approximation, inference for the parameters can be based on the MLEs and their estimated variances, which are available from Fisher information matrix. Then, considering appropriate conditions on the model, to obtain interval estimates for ϑ = (α, θ) we can act as though (θ − ϑ) N(0, I −1 (θ)), where I(ϑ) denotes the observed information matrix and I −1 (θ) = I −1 (ϑ)| ϑ=(α,θ) T . The components of the observed information are given by I kl = −∂ 2 /∂ϑ k ∂ϑ T l , k, l = 1, 2, and they are provided in Appendix A. Assessment of fit and model selection is also a important issue. The PEre model (2) includes as particular case the HPP. We may be interested in verify if such simpler model could be considered. Thus, we may test the hypotheses H 0 : θ = 0 versus H 1 : θ > 0, which leads to the particular case of Equation (2). Note that, for the H 0 , θ = 0 is on the boundary of the parameter space of θ.
We consider the likelihood ratio and the score tests. Letθ = arg max (α,θ) L(α, θ) the MLEs obtained from fitting full model andθ 0 = arg max (α,θ =0) L(α, θ) the corresponding MLE obtained under the restricted hypothesis H 0 , both under a sample of size n. The likelihood ratio statistic (LRS) is defined by where (ϑ) is the log-likelihood function. Following Maller and Zhou [23], since the test is performed in the boundary of the parameter space, the LRS, D 2 , is assumed to be asymptotically distributed as a symmetric mixture of a chi-squared random variable with one degree of freedom (χ 2 1 ) and a mass at zero. This result can be written as lim n→∞ P(D 2 ≤ x) The score test is defined by Peng and Xu [27] as where A = −∂ 2 (ϑ)/∂θ∂α and B = −∂ 2 (ϑ)/∂α 2 . The null distribution of the score test can be approximated by the chi-square distribution with one degree of freedom. Again, large values of this statistic, Z 2 , provide evidence against H 0 . In addition, checking the model fit is also critical to ensure that assumptions underlying the model are plausible to the available data. However, informal graphical methods can provide information, and in this sense we consider the Cox-Snell residuals to the assessment of fit. These residuals are defined asÊ ij =ˆ (t i,j−1 , w ij ), i = 1, . . . , n and j = 1, . . . , m i , whereˆ (·) is the cumulative rate function obtained from the fitted model [8]. The Cox-Snell residuals for the PEre model are defined, with the estimates being the MLEs, aŝ The Cox-Snell goodness-of-fit plot is a visual diagnostic tool where the Nelson-Aalen estimate of the cumulative hazard function is generated based on the Cox-Snell residuals, which are used as the lifetimes failure and the original censoring used as failures [6].

Simulation study
In this section we describe the results of a simulation study performed in order to assess the applicability of asymptotic results of the estimation procedure previously described. In addition, we also conducted a simulation study to investigate the empirical power of the hypothesis tests and the behavior of Cox-Snell residuals for recurrent events. The study was based on the generation of data set as described in Section 2, assuming initially that each one of n subjects experienced the same number of event occurrences m i = m = 2, 5, 7, 15 with n = 30, 50 and 100. In the simulation we fixed α = 3, but the assessment of the inferential procedure does not depend on the specific value of α selected for data generation. We considered different values of θ , θ ∈ {0.75, 2, 4}, which correspond to different degrees of inclination of the rate function. Then we simulated 36 different cases, each one with 1000 samples. Besides, to address the impact of censoring, the study was repeated with right-censored samples with m = 5, 7 and 15. For m = 5 we consider one censored observation, for m = 7 we consider two censored observations and for m = 15 we consider four censored observations, which correspond approximately to 25%, 30% and 30% of censoring.
Then, we first evaluate the accuracy of the approximation of the variance and covariance of the MLEs obtained from the Fisher information matrix. Table 1 shows the simulated values of Var(θ), Var(α) and Cov(θ ,α) as well as the approximate values from the average the corresponding values determined through the observed information matrix. It is noted that the approximated values obtained from the observed information matrix are close to the simulated ones for a small and moderate numbers of recurrence. Moreover, it is observed that the approximation becomes quite accurate when the sample size and the number of recurrences per subject increase.
In addition, simulation studies have been performed in order to analyze the frequentist properties of the estimation procedure. To examine the frequentist properties, we constructed the 95% confidence intervals for the model parameters and calculated their coverage probabilities (CP). The averages of the 1000 MLEs as well as their standard errors and the empirical CP for different sets of parameters, different sample sizes and different numbers of recurrence are condensate in  Table 2. It can be observed that the empirical CP are closer to the nominal coverage level for a larger sample size and moderate numbers of recurrence. Moreover, as the number of observations increase the estimates approaching the real value of the parameter and the standard errors of MLEs decrease, which allows us to conclude that the established estimators of parameters are consistent. Thus, we can guarantee the form of data generation given by Equation (11). The results are similar for all sets of parameters. We also analyze the likelihood ratio test and the score test under the null hypotheses, H 0 : θ = 0, and their power to detect the alternative hypothesis. The results are summarized in Tables 3  and 4, which report, respectively, the observed proportions of type I error and the empirical powers from the likelihood ratio test and the score test at a 5% nominal significance level. The rejection regions are obtained from a symmetric mixture of a chi-squared random variable with one degree of freedom and a point-mass at zero for the likelihood ratio test D 2 and from the chi-square distribution with one degree of freedom for the score test Z 2 .
It can be seen that empirical significance levels of the likelihood ratio test are very close to the nominal level while those from the score test are slightly greater than the nominal level under 30 and 50 sample sizes, particularly for recurrence numbers lower than 15. As the sample size increases, the empirical levels of the tests all approaching the nominal level. The larger sample sizes also correspond to larger power for the tests, but small and moderate numbers of recurrence do not harm the empirical powers. Furthermore, it can be concluded that the empirical power of the likelihood ratio test is greater than those from the score test for values of θ closer to zero, and this difference decreases when the parameter θ is farther to zero. When θ = 1.5, a larger power than 90% is achieved, for both tests, when the sample size is 100. However, for θ = 3 an even larger power than 90% is achieved even for a small sample size. Finally, we also assess the fit of proposed model, based on the simulated data, using the Cox-Snell residuals. Figure 2 presents the Cox-Snell goodness-of-fit plots for the proposed model, for different numbers of occurrences m = 2, 7 and 15 and sample sizes n = 30, 50 and 100. Similar results hold for k = 5 recurrences. To obtain these residuals, we have calculated the cumulative Cox-Snell residual. This is necessary since the recurrent event data require a cumulative residual for each occurrence, given that in the last recurrence of each subject is recorded a cumulative residual total. A good fit is observed when the points have a closely fitting to the diagonal.
The results presented in Figure 2 have no departure from linearity, even with a small or moderate number of recurrences and the presence of censored observations, reflecting a satisfactory fit of the proposed model. Again, we can also ensure the form of data generation given by Equation (11).

Bowel motility data
In this section, the methodology is illustrated in a real data set which concerns the occurrence of certain cyclic movements in the small bowel during the fasting state. We consider the data discussed by Aalen and Husebye [2], which describe a study of small bowel motility (muscular activity) involving 19 healthy human subjects. Initially, catheters were positioned nearby the small bowel to monitor intraluminal pressure. Afterward, the individuals were examined continuously from 5:45 p.m. to 7:25 a.m. on the next day, giving a total of 13 h and 40 min observation. Then, a standardized mixed meal was given to each individual at 6:00 p.m. to induce a fed state. After a time period in which irregular bowel contractions occur, a fasting (interdigestive) state begins, with a cyclical bowel activity pattern. The time between two consecutive fasting cycles is also called the migrating motor complex (MMC) period. The lengths of successive motility (or MMC) cycles make up the data set, so that the motility cycle patterns are not associated with the duration of the fed state. In this study, events are associated with the beginning of a motility cycle and the length of the last cycle for each individual is censored by the end of the monitoring period.
Intestinal dysmotility is an alteration in bowel functions and occurs due abnormalities in the muscles and nervous of the bowels. As pointed out in Section 1, even though the complementary risks are unobserved, one may speculate on some possible complementary risks for the bowel motility data. For example, the bowel motility can be affected by diet, drugs, diseases, such as neuropathy, autonomic neuropathy, amyloidosis, myopathies, and other causes that even unobserved lead to the intestinal dysmotility [39]. In this context, the PEre model allows to analyze the effect of these factors as latent complementary risks for a better interpretation and adaptation to the situation. Then the PEre model (2) was fitted to the bowel motility data, as well as its HPP particular case. Furthermore, by fitting the submodel of (2) we can evaluate whether the gap times for an individual are possibly independent. Table 5 presents the MLEs and the corresponding 95% confidence intervals (in parenthesis) of the parameters, which were based on the observed information matrix, for both PEre and HPP models. The last column of this table provides the log-likelihood for both models.
The results of Table 5 show that both PEre and HPP models are well supported by the data. However, the HPP model assumes that the gap times of an individual are independent. Thus, for the considered data, where the assumption of independence between events may be inappropriate, it is expected that the PEre model is more advantageous with respect to its particular case (HPP model), as can be seen in Table 5.
Therefore, the advantage of the proposed model is that it can be used both in situations where the events are independent as those in which there is no certainty of independence between  recurrent events of an individual, providing a wider applicability relative to the classical HPP model and also to other renewal processes. Testing of independence between gap times can be done based on likelihood ratio and score tests by considering the null hypothesis is H 0 : θ = 0. The LRS, of the hypothesis that θ = 0 gives an observed value for the D 2 of 12.713. The 95th percentile, x 0.95 , of the asymptotic distribution of the LRS can be calculated from 1 2 + 1/2P(χ 2 1 ≤ x 0.95 ) = 0.95, so P(χ 2 1 ≤ x 0.95 ) = 0.9, giving x 0.95 = 2.71 [23]. Since 12.713 > 2.71, we reject H 0 , with a p − value < 0.001, and considered there is evidence in favor of the proposed model. Similarly, the score test provides a value for the Z 2 equal to 13.710, which is compared with a chi-square distribution with one degree of freedom, χ 2 1,0.05 . Since χ 2 1,0.05 = 3.841 and 13.710 > 3.841, we reject H 0 , with a pvalue < 0.001, and again we considered there is evidence in favor of the PEre model. Therefore, the value ofθ = 4.119 is significantly greater than zero, with p-value < 0.001 for both tests.
In order to verify the overall goodness of fit of the PEre model we calculated the Cox-Snell residuals. Figure 3 shows theCox-Snell goodness-of-fit plots for the PEre and HPP models based on the bowel motility data. Again, the advantage of the PEre model is demonstrated by closely fitting to the diagonal.
In addition, the proposed model is a satisfactory option for data analysis, pointing out the importance of marginal features on the modeling of recurrent event data and the effect of unobservable causes that also contributes to the occurrence of the event.

Final comments
In this paper, we consider situations where individuals in a longitudinal study are subject to recurrences of an event of interest. The proposed model is an application of the Poisson-exponential model, proposed by Cancho et al. [4], for recurrent event data, where we have investigated the gap times between events. We obtained the conditional distributions of gap times from the rate function, which is an attractive formulation for recurrent event data with direct interpretations frequently preferred by professionals on medical and biostatistical areas [41]. The parameters of PEre model have a direct interpretation in terms of complementary risks. Moreover, its specific parametric form is analytically convenient and easily interpretable. We discussed MLEs, obtained by direct maximization of the log-likelihood function. The inferential procedure based on the maximum likelihood approach is implemented straightforwardly. The results of simulation study showed the effectiveness of the parameter estimation approach, allowing us to conclude that small and moderate numbers of recurrences do not affect the estimates of the parameters, even in presence of censoring. The PEre model allowed a straightforwardly nested hypothesis testing procedures for comparison with its HPP particular case. The Cox-Snell residual plots allowed to assess the model fit and ensure that assumptions underlying the model are plausible to the available data. The applicability of the PEre model was demonstrated in a real data set.
The analysis considered in this paper is a preliminary step for the development of a more complete model, and our modeling may be generalized in some directions. An immediate extension of this model is to consider the effects of covariates for each subject and also a possible heterogeneity between them. As an alternative to the Poisson-exponential distribution other distributions, such as Weibull-geometric [3], beta-Weibull geometric [9], Geometric Birnbaum-Saunders [5], can be considered for modeling. These issues are of interest and will be investigated in future.