Joint modeling of interval counts of recurrent events and death

When a recurrent event process is ended by death, this may imply dependent censoring if the two processes are associated. Such dependent censoring would have to be modeled to obtain a valid inference. Moreover, the dependence between the recurrence process and the terminal event may be the primary topic of interest. Joint frailty models for recurrent events and death, which include a separate dependence parameter, have been proposed for exactly observed recurrence times. However, in many situations, only the number of events experienced during consecutive time intervals are available. We propose a method for estimating a joint frailty model based on such interval counts and observed or independently censored terminal events. The baseline rates of the two processes are modeled by piecewise constant functions, and Gaussian quadrature is used to approximate the marginal likelihood. Covariates can be included in a proportional rates setting. The observation intervals for the recurrent event counts can differ between individuals. Furthermore, we adapt a score test for the association between recurrent events and death to the setting in which only individual interval counts are observed. We study the performance of both approaches via simulation studies, and exemplify the methodology in a biodemographic study of the dependence between budding rates and mortality in the species Eleutheria dichotoma.

In Section 4, we assess the performance of the estimation method and the test in simulation studies. In Section 5, we apply the proposed methods to the data on E. dichotoma, followed by a discussion in Section 6.

JOINT FRAILTY MODEL AND INTERVAL COUNTS
Before presenting the estimation and test procedure in the next section, we will introduce in the following the joint frailty model, which allows us to model the dependence of the recurrent event process and the terminal event. We then derive the likelihood function for data that only contain interval counts of the recurrent events.
We consider a sample of independent individuals denoted by , = 1, … , . Each individual is observed from time 0 = 0 until the end of its follow-up . The time is either a censoring time , which is independent of the recurrent event process and the terminal event, such as end of study; or it is the time of the terminal event, whichever comes first: = min ( , ). For simplicity, we assume that the terminal event is death, and denote by = The joint frailty model for recurrent events and death, as proposed by Liu et al. (2004), is then specified as with baseline rates 0 ( ) and ℎ 0 ( ) that are common to all individuals. The shared frailty enters both the recurrent event rate and the hazard rate of the terminal event, thereby introducing both dependence between the recurrences within one individual, as well as the association between the recurrences and the terminal event. The parameter determines the direction and the strength of the association between the two processes. If > 0, a higher rate of recurrence implies a higher mortality risk; if < 0, a higher rate of recurrence implies a lower mortality risk. If = 0, the rate of recurrence does not affect the mortality risk.
The frailties are often assumed to follow a gamma distribution with mean one and variance . In the following, we will more generally assume that the stem from a distribution with density ( ) with parameter .
The covariates enter in (2) in a proportional rates or a proportional hazards formulation, but can have different effects and on the terminal event and the recurrence process, respectively. In the simplest setting, the exact times of event occurrence are observed. In many cases, however, only the number of events that occurred in a sequence of consecutive time intervals is available. More precisely, we observe individual interval counts as realizations of = ( ) − ( −1 ). shows the number of recurrent events experienced by individual in the interval = ( −1 , ], = 1, … , . correspond to the observation times of individual , for instance, times of hospital visits in medical studies or, as in our example, generated by the lab logistics. Thus, both the positions of the and the total number of intervals can vary across individuals. The follow-up times are still exactly observed so that = . As the frailties are unobservable, the inference is based on the marginal likelihood that is obtained by integrating the conditional likelihood given the frailties over the frailty distribution ( ). The conditional likelihood of the joint frailty model (2) based on exactly observed recurrence times was developed in Liu et al. (2004). For the current setting, the likelihood factor for the recurrence times (formula (7) in Liu et al., 2004) has to be replaced by the contribution of the interval counts of the recurrent events. From (1) and (2), it follows that ( | , , ≥ ) has a Poisson distribution with mean , where = ∫ ′ 0 ( )d . Therefore, the likelihood contribution ( ) ( ) of individual conditional on its frailty value is given by This leads to the marginal likelihood contributions In general, this integral does not have a closed-form expression.

METHODS
In the first part of this section, we elaborate the estimation procedure of the joint frailty model based on interval counts with different individual observation intervals. Then, in Section 3.2, we demonstrate how the score test for dependence between the two processes can be adapted to the case of interval-count data.

Estimation of the joint frailty model based on interval counts
For observed recurrence times, Liu and Huang (2008) suggested using Gaussian quadrature to approximate the marginal likelihood of the joint frailty model. The approximated likelihood can then be maximized directly because the integral is replaced by a weighted sum of function values. More specifically, Liu and Huang (2008) applied Gauss-Hermite quadrature, which for a function ( ) uses the approximation The quadrature points are the roots of the th -order Hermite polynomial, and are the corresponding weights. This approach is applicable to marginal likelihoods that are integrated over normal random effects, for which with the standard normal density (⋅), and modified quadrature points̃= √ 2 together with weights̃= √ 2 2 . For nonnormal random effects, Liu and Huang (2008) used the probability integral transformation proposed by Nelson et al. (2006). If the random effect density is ( ) with corresponding distribution function ( ), then the integral over the density ( ) is transformed into an integral over standard normal random effects. This is achieved by noting that = Φ −1 ( ( )) follows a standard normal distribution if the ( ), which follow a standard uniform, are transformed by the inverse of the standard normal distribution function Φ(⋅).
We apply this quadrature approach to the marginal likelihood of the joint frailty model based on interval counts of recurrent events; see Equation (3). Substituting = −1 (Φ( )) in the marginal likelihood contributions, we obtain These can directly be approximated using Gauss-Hermite quadrature as with̃and̃as defined in (4). The approximate marginal likelihood of the joint frailty model is then given by To actually maximize the approximate likelihood (5), we have to specify the baseline rates 0 ( ) and ℎ 0 ( ), as well as the frailty distribution ( ). Similar to Liu and Huang (2008), we model the baseline rates as piecewise constant functions This choice is particularly suitable if no prior knowledge of the shapes of the two rates 0 ( ) and ℎ 0 ( ) is available. The specifications of the intervals = ( −1 , ], = 1, … , , and = ( −1 , ], = 1, … , , can differ between the recurrent event process and the death process in terms of both their lengths Δ = − −1 and Δ = − −1 and their total numbers and . (The intervals for the piecewise constant baseline rates should not, however, be confused with the intervals in which the numbers of recurrent events are observed; see Section 2.) The baseline rate 0 ( ) of the recurrence process enters the likelihood (5) through the conditional means of the interval counts given the frailty value . Under the piecewise constant rate model, is computed as Piecewise constant baseline rates offer more flexibility than parametric models, such as the Weibull model, while at the same time remaining more tractable than purely nonparametric models. Previous studies have suggested that a moderate number of intervals, that is, between 8 and 10 intervals, yields satisfactory estimation results (Lawless & Zhan, 1998;Liu & Huang, 2008). The performance of the piecewise constant model is usually improved when the interval cut-points are based on quantiles of the recurrence and the survival times, respectively. In the current setting in which only interval counts of recurrent events are observed, the exact recurrence times are unknown. If, however, the individual observation intervals are relatively short compared with the total follow-up, we can approximate quantiles by creating a set from the observation times , with each repeated times, = 1, … , , = 1, … , ; and then determining the cut-points as quantiles of this set of times.
Parameter estimation in the joint frailty model is then done by maximizing the approximate marginal log-likelihood; that is, the logarithm of (5). The standard errors for the parameter estimates can be obtained from the inverse of the negative Hessian of the approximate marginal log-likelihood. Further computational details can be found in Section A.1.

Score test for the association between recurrences and death
The joint frailty model for recurrent events and death is rather complex, with estimation procedures that are more involved than those needed for the fitting of two separate models, one for the recurrent event process and one for the survival process. Investigating whether the two processes are associated-and, consequently, whether joint modeling is required-is a useful first step. Moreover, the question of whether there is an association-and, if so, whether it is positive or negativecan be a stand-alone issue, not necessarily followed by fitting a joint model. Balan et al. (2016) proposed a correlation score test for the association between recurrences and terminal events in the joint frailty model for settings with observed recurrence times. Their test can be performed by fitting separate models for the recurrence times and the survival data, and, thus, without fitting the joint frailty model. In addition, the sign of the test statistic is an indicator of the direction of the association. In the following, we show that this score test can be adapted to the setting with interval counts of recurrent events, provided the survival times are exactly observed.
A test for association in the joint frailty model (2) corresponds to a test of 0 ∶ = 0 against 1 ∶ ≠ 0. All other parameters, including those for the baseline rate models, are treated as nuisance parameters, and denoted by . The score test of Balan et al. (2016) is based on the score function for under the null hypothesis; that is, where is the marginal log-likelihood of the joint frailty model. The authors showed that the score, evaluated at the maximum likelihood estimate under 0 , that is, (0,̂0), is proportional to the covariance of the estimated martingale residuals of the terminal event and the "posterior" estimates of the log-frailties for the recurrent events given the observed data. More formally, defining with the cumulative baseline rates Λ 0 ( ) = ∫ 0 0 ( ) d and 0 ( ) = ∫ 0 ℎ 0 ( ) d , Balan et al. (2016) derived that Theˆare estimates of the martingale residuals of the terminal event model log ( ) are the posterior estimates of the log-frailty values given the observed data from the recurrent event process: For gamma distributed frailties with mean one and variance , one obtainŝ where (⋅) is the digamma function.
Because of the zero-mean constraint of theˆ, the second line of (8) is proportional to the correlation = cor(ˆ,log( )) between the martingale residuals and the estimated log-frailties. Consequently, Balan et al. (2016) based the correlation score test on the test statistic which, under the null hypothesis, asymptotically follows a -distribution with − 2 degrees of freedom.
It turns out that in the setting in which only interval counts of the recurrent events are available, but exact survival times (or censoring times) are observed, Equation (8) still holds. We show this result in Section A.2. Therefore, the test statistic in (10) is still valid. Furthermore, for gamma distributed frailties, the estimatesl og ( ) can still be determined using formula (9). The latter formula involves estimateŝof the frailty variance, the covariate effect̂on the recurrence rate, and the cumulative baseline rateΛ 0 ; all determined under 0 . These can be obtained by fitting a mixed Poisson model to the interval counts of the recurrent events (see, for instance, Lawless & Zhan, 1998, who assume a piecewise constant baseline rate function 0 ). We generally recommend to estimate the shared frailty model for the interval counts using a flexible specification such as (6) for the baseline rate. A simple parametric model like, for example, the Weibull model, although appealing due to its parsimony, always bears the risk of misspecification and consequently misleading test results. The estimates of the martingale residualsˆcan be derived from a Cox proportional hazards model fitted to the survival data { , , ; = 1, … , }.

Performance of the estimation method
To evaluate the performance of the proposed method for estimating the parameters of the joint frailty model based on interval counts of recurrent events and survival times, we conducted a simulation study. Several different aspects will affect the estimation results, both on the part of the model specification but also on the part of the observable data. The latter include the number and lengths of the intervals for the recurrent event counts and whether they are the same for all individuals in the sample or not. The amount of independent censoring will also have an impact. Among the former aspects, the size of the frailty variance and the sign of the dependence parameter are expected to influence the estimation, but also the number of intervals used in the piecewise constant specification of the baseline rates (and consequently the total number of parameters to be estimated) will matter.
We generated data for = 200 individuals from the joint frailty model (2). For the baseline rates 0 ( ) and ℎ 0 ( ), we chose the form of Weibull hazards, ( ∕ ) ( ∕ ) −1 , with shape parameter equal to 1.5 or 3 and scale parameter equal to 1∕3 and 1.35, respectively. A single binary covariate that takes values 0 or 1 with probability 0.5 was included, and had the same effect on the hazard of death and the rate of recurrence ( = = 1). Frailties were simulated from a gamma distribution with mean one and different variances ∈ {0.25, 0.5, 0.75}. We considered both cases of positive ( = 1) and negative association ( = −1) between the recurrence process and the terminal event. Additionally, we looked at one setting with a relatively small value for the frailty variance, = 0.05, that was inspired by the results of the data set on E. dichotoma, see Section 5. In this setting two values for the dependence parameter | | = 1 or | | = 5 were studied.
Independent censoring was considered in two versions: either an end-of-study censoring at time = 2 for all individuals still alive then or individual censoring times, which occurred uniformly over the total follow-up window [0,2].
Regarding the observation times , which determine the intervals during which the recurrent events are counted, we examined two scenarios. In Scenario I the observation times were the same for all individuals. We set = 0.2 , = 0, … , 10, such that we observed individual interval counts in up to 10 intervals of equal length 0.2. We also studied one scenario with = 0.1 , = 0, … , 20, leading to up to 20 individual interval counts. In Scenario II the observation times varied across individuals. This scenario mimicked a study in which participants have scheduled visit times, but actually could be early or late for their visits. This was implemented as follows: the scheduled times were 0 = 0.2 , = 0, … , 10, but the actual observation times for each individual were obtained by adding a random noise to the scheduled times so that = 0 + . The were drawn from a uniform distribution on [−0.1, 0.1] (except for 0 = 0 and 10 as uniform on [−0.1, 0]) and hence the varied across individuals.
For the estimation of the joint frailty model, we assumed that the frailties were gamma distributed, and that the baseline rates were specified as piecewise constant functions on 10 intervals of equal length, = = 0.2 . Thus, under Scenario I and = 0.2 the intervals for the rate pieces coincided with the observation intervals for the counts. We also studied the impact of an increase of the number of intervals to 20. Approximation of the likelihood was based on = 30 quadrature points. All computations were run in R (R Core Team, 2019), see Section A.1 for further details. In all settings we ran 200 replications.
The results of the simulation study under Scenario I with frailty variance = 0.5 are shown in Figures 1 and 2, both From the box plots of the estimates of the covariate effects and , the dependence parameter , and the frailty variance in Figure 1, we can see that the method performs satisfactorily. Censoring scheme A is considerably milder than scheme B (for an illustration, see Table S2), and the corresponding loss in information is reflected in a moderately increased variability in the estimates and, as could be expected, also in the estimated standard errors (Figure 1, bottom). For the frailty variance , variability in estimates and estimated standard errors is higher for negative dependence ( = −1) than for positive dependence. This is also true for the dependence parameter itself. The estimated standard errors for the covariate effects and , the dependence parameter , and the frailty variance are compared with the empirical standard deviations of the respective estimates across the replications in the bottom panels of Figure 1. The general magnitude of the standard errors is well captured. F I G U R E 2 Estimates (gray, solid) of the cumulative rate of recurrence (left) and of death (right) based on samples of size = 200, which were generated from a joint frailty model with positive dependence for two schemes of independent censoring. Top: end-of-study censoring at = 2 (CensA), bottom: uniform censoring on [0,2] (CensB). Red, solid line gives the true cumulative baseline rate; black, dashed line is the mean of the 200 estimates Figure 2 displays the estimates of the cumulative baseline rates for positive dependence for both censoring schemes A and B. The averages of the estimates are very close to the true underlying rates. The stronger loss of observations in scheme B leads to markedly increased variability of the estimates toward the end of the follow-up window, which is an obvious consequence of the decreasing number of individuals under study over time. This low number of observations results in part from individuals experiencing a terminal event. Additionally, even fewer individuals (in this setting, less than 5%) are observed in the interval [1.8,2), which corresponds to the last piece of the baseline rates, due to independent censoring under scheme B. In contrast, censoring scheme A is benign.
For the other settings of the simulation study we restricted our attention to the independent censoring scheme B. With its uniformly distributed censoring times it produces a relatively challenging loss of observations. So the results we report here are conservative in the sense that they hold despite this appreciable amount of censoring. The Supporting Information illustrates the results of the other simulation settings analogous to Figures 1 and 2. Here we summarize the main results of the various scenarios.
If we change the frailty variance to = 0.25 or = 0.75, we note that also in these settings the biases of the parameter estimates, if present at all, are small. In general, varying the true underlying frailty variance mainly affects the variability of the estimates. The estimates of the covariate effects and the frailty variance show increased variability for larger frailty variance. In contrast, the estimates of the dependence parameter show less variability for larger frailty variance because the increased heterogeneity among the individual patterns makes it easier to assess the dependence between the recurrence process and survival.
If the true frailty variance is close to zero ( = 0.05) the estimates of the dependence parameter are still only modestly biased but the variability of thêincreases strongly. This had to be expected because the identification of the dependence hinges on the variation in the individual frailty. If the true value of | | = 1, then the estimated values may turn out with the wrong sign of the dependence parameter, particularly if the true dependence is negative. Thus, weak dependence is difficult to identify in case of low variability in frailty. However, if the dependence is strong, | | = 5, then, despite the variability of the estimates, the sign of the dependence is correctly estimated.
If we modify the width of the observation intervals from 0.2 to 0.1 we first note the following: if we keep the specification of the baseline rate of recurrence as piecewise constant on the 10 intervals with = 0.2 as before, then the additional information provided by the finer observation intervals would not be used by the estimation method (see Section A.3 for a proof). Therefore, we also specify the baseline rates as piecewise constant functions on 20 intervals, that is, = = 0.1 . Apart from increased variability in the rate estimates near the end of the follow-up window, there is little change to the estimation results with width 0.2.
If we allow the observation intervals to vary across individuals (Scenario II), then we find that the method performs equally well as for fixed observation intervals, with the exception of the estimates for the last baseline rate pieces. In this simulation, we used fixed intervals for the baseline rates though, to unify the presentation of results. In applications we would recommend to choose the cut-points and based on (approximate) quantiles of the event times, as described earlier, to increase precision of the baseline rate estimates.
If we replace the baseline rates by a parametric Weibull specification (which here is correct), then there is only very modest change in the other parameter estimates of the model as well as their standard errors, so the main advantage is in the less variable estimation of the two baseline rates. This is, however, counterbalanced by the risk of model misspecification so we advise the use of piecewise constant rates unless one has a good understanding of the underlying processes that justifies a parametric choice.
Overall, the proposed method for estimating the joint frailty model based on interval counts yielded reliable results in this simulation study and different scenarios behaved in the way we had anticipated beforehand.

Performance of the score test
In a second set of simulations, we investigated the performance of the score test in the setting with interval counts of recurrent events. For this purpose, we again generated data for = 200 individuals from the joint frailty model (2). The covariate effects, values for the variance of the gamma frailties, and baseline rates were the same as in Section 4.1. Here, however, we considered not only different directions but also different strengths of the association between the recurrence process and the terminal event; namely, ∈ {−1, −0.5, 0, 0.5, 1}. Counts were again observed in 10 intervals of equal length 0.2 (Scenario I) or of varying length (Scenario II). We studied two schemes of independent censoring, scheme A with censoring time = 2, and scheme B with uniformly distributed censoring times on [0,2]. We ran 1000 replications for each setting to determine the size or power of the test.
The score test involves fitting separate models to the recurrence data and the survival data. First, we fitted a shared gamma frailty model to the individual interval counts of the recurrent events, assuming that ( | ) follows a Poisson distribution with meañ, wherẽ= ∫̃′̃0( ) d . The baseline recurrence ratẽ0(⋅) was modeled as piecewise constant, as in (6), with pieces defined by the cut-points = 0.2 , = 0, … , 10. The estimatesl og ( ) were then determined according to formula (9). Second, we estimated a Cox proportional hazards model from the survival data { , , ; = 1, … , } using function coxph() from package survival (Therneau & Grambsch, 2000) to obtain the martingale resid-ualsˆof the terminal event. Finally, we calculated the test statistic based on the correlation between theˆand thê log ( ). Table 1 reports the size and the power of the score test, performed at a level of 5%, depending on the true underlying dependence parameter and frailty variance . The proportion of a type I error, that is, falsely rejecting the hypothesis of no dependence ( = 0), was affected most strongly by the censoring scheme. If = 2, which implies a modest proportion of independently censored cases (see Table S4), then the level of the test is met or exceeded only slightly. For strong independent censoring (scheme B), which implies that roughly half of the observations are censored and quite some of TA B L E 1 Power and size of the score test, performed at the 5% level, in the joint frailty model with varying dependence parameter ∈ {−1, −0.5, 0, 0.5, 1}, frailty variance ∈ {0.25, 0.5, 0.75}, and independent censoring ∼  [0, 2] or = 2, across 1,000 replications each them early during the follow-up, then this loss of information increases the proportion of a type I error, particularly for the frailty variance = 0.75. This is a noteworthy result, and hence the test should be regarded with caution, if there is strong (and early) censoring in the data.
Regarding the power of the score test, as expected, we found that the power increased with the strength of the dependence, that is, | |, and with the magnitude of the frailty variance. In the settings with the larger frailty variances = 0.5 and = 0.75 and stronger dependence | | = 1, the score test detected the association in all cases. We note that for the samples for which the association was detected by the score test, the direction of the dependence was always identified correctly. An extension of the simulation settings suggested that the score test can detect associations even for small values of the frailty variance as long as the association is sufficiently strong, that is, | | is sufficiently large (see Table S5).
Furthermore, we assessed the performance of the score test for Scenario II in which the observation intervals vary across individuals. Results for the setting with 10 scheduled visit times and 10 pieces for the baseline rate are displayed in Table S8. The proportion of false rejections of the hypothesis of no dependence are again a bit higher than the nominal level due to the high percentage of censoring. The power of the test is comparable to the values obtained in Scenario I, although the power in the settings with negative dependence is a bit lower here.
Additional results on the performance of the score test for modifications of Scenario I and II are given in the Supporting Information.
In conclusion, the results of the simulation study provide further evidence that the score test is a powerful method for assessing the association between recurrent events and the terminal event, also in a setting in which only interval counts of the recurrent events are observed.

FERTILITY AND MORTALITY IN ELEUTHERIA DICHOTOMA
To illustrate the proposed methods, we use data from a biodemographic study on the fertility and mortality of E. dichotoma (Dańko et al., 2020), which was briefly introduced in Section 1. Eleutheria dichotoma is a marine organism that passes through several life-cycle stages: that is, planula larva, polyp, and medusa stages. The colonial polyps (hydroids) asexually produce medusae. The medusae can then reproduce both sexually (by producing larvae) and asexually (by producing medusa buds). In this study, asexually budded medusae were collected directly from the hydroid colony and reared for three generations-each of which was, in turn, obtained through the asexual budding of medusae. In our analysis, we focus on the age trajectory of the budding rate (asexual reproduction) and the mortality of one of the medusa generations, as well as on the association between the two processes.
Age 0 = 0 of a medusa corresponds to the point in time when it detaches from its hydroid colony or from its ancestor medusa. The medusae were followed individually, and were observed until death or censoring. This occurred when either the study ended, a laboratory accident (e.g., water evaporation) led to the loss of the medusa, or the medusa was absorbed by a large bud of the same individual. The animals were checked for newly released larvae and buds roughly three times per week. The resulting observation times differed across individuals, with interval lengths between 1 and 11 days.
Salinity is an important factor that affects the physiological responses of species like E. dichotoma, both at the level of the hydroid colony that produced the medusae under study, and at the level of the medusae themselves. Four combinations of salinity levels were studied here (low(hydroid)-low(medusa), medium-medium, low-medium, and medium-low); for more details, see Dańko et al. (2020). The data set contains = 141 individuals, with the following group sizes in the four experimental conditions: 36 lowlow, 40 medium-medium, 32 low-medium, and 33 medium-low. The individuals produced between 0 and 27 buds over their life course, with a mean of 8.99. Follow-up times varied between 9 and 217 days, with a median of 98 days; and 34 individuals (24%) were censored. Figure 3 exemplarily shows the data for 14 individuals.
To assess whether the recurrent budding process and survival are associated in E. dichotoma, we first conducted the score test that was presented in Section 3.2. For fitting the shared frailty model to the individual bud counts, we assumed that the frailty variable was gamma distributed. Moreover, for the piecewise constant baseline budding rate, we defined cut-points at 0, 16,22,27,32,40,47,63,85,111, and 220 days. These cut-points were obtained as (approximative) deciles of the recurrent event times, as described at the end of Section 3.1. We included two binary covariates for the salinity levels at the polyp stage and the medusa stage, respectively.
The martingale residuals of the terminal event were obtained from a Cox proportional hazards model fitted to the survival data, including the two covariates on the salinity levels.
The parameter estimates of the two separately fitted models, the Cox model for the survival data, and the shared frailty model for the budding rate based on interval counts, are shown in Table 2.
The correlation between the martingale residuals from the Cox regression and the estimated log-frailties from the budding model is = −0.434, yielding a test statistic of = −5.676 with a -value of 7.729 ⋅ 10 −8 . Thus, the result of the score test clearly suggests that reproduction and mortality are negatively associated, that is, that a higher budding rate is associated with lower mortality, which implies that a joint model should be used to analyze these data on E. dichotoma.
Therefore, we estimated a joint frailty model for reproduction and mortality for these data, including the two covariates on salinity. Frailties were again assumed to stem from a gamma distribution. For the baseline budding rate, the same specification was used as in the separate model (see above). For the baseline hazard of death, we used a piecewise constant F I G U R E 4 Estimated baseline rates of budding (left) and death (right) for the E. dichotoma data set function with cut-points at 0, 60,67,72,85,98,103,114,132,153, and 220 days; again, the cut-points were taken from approximate deciles of the survival times. For the Gaussian quadrature in the marginal likelihood (see Equation (5)), = 30 quadrature points were used. The parameter estimates of the joint model are also displayed in Table 2. Interestingly, the dependence parameter was estimated as −4.941 along with a frailty variance of 0.051. The negative value of̂indicated that individuals with higher rates of asexual reproduction tended to have a lower mortality risk than individuals with lower rates of asexual reproduction. This finding that higher rates of reproduction were accompanied by longer survival stands in contrast to the idea of a trade-off between reproduction and survival. Regarding the salinity levels, we found that the salinity experienced by the polyps did not have a noticeable effect on the reproduction or survival of the medusae. In contrast, individuals who were exposed to low salinity at the medusa stage were found to have both higher fertility and mortality rates than those exposed to medium salinity. Finally, Figure 4 shows the estimates of the age-specific budding and death rates. The budding rate peaked at about 20 days of age before gradually declining to a nonzero level at later ages. The death rates increased markedly after about 60 days of age.
To conclude, although the results presented here are in agreement with the findings of Dańko et al. (2020), they also provide additional insight into the association between the asexual reproduction and survival of E. dichotoma. In addition to the biological implications of the dependence itself, due to the association between the two processes, the recurrent event process is censored by the terminal event (death) nonindependently, which warrants the joint modeling of the two processes.

DISCUSSION
We have presented a method for estimating the joint frailty model for recurrent events and death in situations in which only individual interval counts of the recurrence process are observed. When modeling the baseline rates as piecewise constant, the marginal likelihood can be approximated using Gaussian quadrature, and can then be maximized directly. In addition, we have shown that the score test for the association between recurrences and death (Balan et al., 2016) is also applicable in the setting with interval counts. The test is based on the correlation between the martingale residuals of the terminal event and the estimates of the log-frailties, which are obtained by separately fitting a Cox proportional hazards model to the survival data and a shared frailty model to the interval counts of recurrent events. Our simulation studies demonstrated that both the estimation method and the score test perform well. We also found that when applying the proposed methods to data on fertility and mortality in E. dichotoma, the rate of asexual reproduction and the mortality risk are negatively associated. This finding is interesting from a biological point of view, and it also demonstrates the necessity of allowing positive and negative dependence in a model such as the relatively complex joint frailty model (2). Moreover, the E. dichotoma example illustrates the advantages of using the piecewise constant rate model. On the one hand, the shape of the budding rate, which is characterized by a sharp peak at younger ages and a gradual leveling off to a nonzero level at older ages, is hard to capture using a simple parametric model. On the other hand, the data structure of the interval counts, in which the observation times vary across individuals, makes it difficult to construct a purely nonparametric estimate of the baseline rate.
Implementing the estimation method in the joint frailty model involves making several choices, such as decisions regarding the distribution of the frailties, the number of quadrature points, and the number of pieces included in the baseline rate models. For the frailty distribution, we assumed that frailties follow a gamma distribution in both the simulation study and the application. The use of a gamma distribution is a popular choice for the distribution of frailties, and, with respect to the score test, it has the benefit of yielding closed-form expressions for the estimated log-frailties. However, the quadrature approach is equally able to accommodate the log-normal distribution or any frailty distribution that has a closed-form inverse distribution function. Yet, the performance of the quadrature approach relies on the quality of the approximation of the marginal likelihood, which, in turn, depends on the number of quadrature points. Liu and Huang (2008) suggested using = 30 quadrature points for gamma frailty models, and in our experience, this number yields reliable results. In practice, we recommend fitting the model for several increasing values of until the estimates stabilize. For the piecewise constant rate models, using a moderate number of up to 10 intervals for the rates seems to produce good results. Using a larger number of intervals for the baseline rates generates a larger number of parameters to be estimated, which, in turn, increases computational costs, and might affect the numeric stability of the method.
One of the limitations of the approach presented here is that, in the joint frailty model, the association between recurrences and death is modeled via a dependence parameter acting on a shared frailty. If, however, the individuals are not sufficiently heterogeneous, that is, if the frailty variance is not sufficiently large, the dependence parameter is not meaningful, and the association cannot be assessed.
Another restriction is imposed by the piecewise constant rate models. Although these models can capture a variety of different shapes of the baseline rates, using more flexible rates-and, in particular, smooth rates-might be desirable in some applications. Further work is needed on how to incorporate smooth rate models with automatic smoothing parameter selection.
Finally, the observation times in our application are fixed by the experimental setup, even though in real-world applications, particularly in medicine, the observation times may depend on the recurrence process. For instance, patients may visit the doctor more often when they are in worse condition. This is taken into account in a model developed by Zhao et al. (2013), which considers interval counts of recurrent events in the presence of death and a dependent observation process. However, the authors adopted a marginal approach that left the association between the recurrences and death unspecified.

C O N F L I C T O F I N T E R E S T
The authors have declared that there is no conflict of interest.

O P E N R E S E A R C H B A D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available in the Supporting Information section. This article has earned an open data badge "Reproducible Research" for making publicly available the code necessary to reproduce the reported results. The results reported in this article were reproduced partially for data confidentiality reasons and due to their computational complexity.

S U P P O R T I N G I N F O R M AT I O N
Additional supporting information may be found online in the Supporting Information section at the end of the article.

APPENDIX A A.1 Computational details
We implemented our estimation approach in R (R Core Team, 2019) using function gauss.quad() from package statmod (Smyth, 2005) to determine the quadrature points and weights, and function nlm() for numerical optimization of the approximate marginal log-likelihood.
The Hessian of the log-likelihood can be obtained directly from the output of the function nlm(). However, in some cases, computation of the Hessian using function hessian() from package numDeriv (Gilbert & Varadhan, 2019) yields more stable results.
To ensure that the Hessian of the approximate log-likelihood with piecewise constant baseline rates is invertible, it can be necessary to fit the joint frailty model with small, fixed ridge penalties on the logarithm of the baseline rates. The standard errors are then calculated based on the Hessian of the penalized log-likelihood.

A.2 Derivation of the score test
In this section, we show that the score ( , ) of the joint frailty model (2) has the same form independently of whether exact recurrence times or only interval counts of the recurrent events are available, as long as the recurrence process is observed up to exactly known follow-up times.
Let us start by rewriting the individual contributions (3) to the marginal likelihood for interval counts and ( , ) as defined in (7) which is the first line of (8).

A.3 Likelihood with fixed observation times
In this section, we study the likelihood of the joint frailty model (2) based on individual interval counts of recurrent events in one particular setting. Specifically, the observation times are assumed to be the same for all individuals and the baseline rate of recurrence 0 ( ) is modeled as a piecewise constant function. We will show that if the observation intervals are finer than the intervals for the rate pieces, the score depends on the individual's interval counts only through the sums of these counts over each baseline rate piece. For this purpose, suppose the observation times are given by = for = 0, … , − 1, and the last observation time = is equal to the follow-up time . The interval counts , observed over the intervals = ( −1 , ], enter the likelihood contribution ( ) ( ) of individual given its frailty value (see (3) (6), we have where | | denotes the length of the interval . Now let be the index of the interval for which ∈ . We can write the sum in the third term in (A.1) as which depends only on the individual's follow-up time , but not on the other observation times , = 0, … , − 1. Using (A.2), the last term in (A.1) equals Suppose now we choose the intervals for the baseline rate of recurrence to be rougher than the observation intervals, but such that the are a subset of the . Then, each , = 1, … , , is a subset of exactly one , thus, Consequently, if the observation intervals are finer than the rate intervals the likelihood is proportional to an expression that depends only on recurrent event interval counts corresponding to the larger rate intervals.