Additive-Multiplicative Rates Model for Recurrent Event Data with Intermittently Observed Time-Dependent Covariates

Regression methods, including the proportional rates model and additive rates model, have been proposed to evaluate the eﬀect of covariates on the risk of recurrent events. These two models have diﬀerent assumptions on the form of the covariate eﬀects. A more ﬂexible model, the additive-multiplicative rates model, is considered to allow the covariates to have both additive and multiplicative eﬀects on the marginal rate of recurrent event process. However, its use is lim-ited to the cases where the time-dependent covariates are monitored continuously throughout the follow-up time. In practice, time-dependent covariates are often only measured intermittently, which renders the current estimation method for the additive-multiplicative rates model inappli-cable. In this paper, we propose a semiparametric estimator for the regression coeﬃcients of the additive-multiplicative rates model to allow intermittently observed time-dependent covariates. We present the simulation results for the comparison between the proposed method and the simple methods, including last covariate carried forward and linear interpolation, and apply the proposed method to an epidemiologic study aiming to evaluate the eﬀect of time-varying streptococcal infections on the risk of pharyngitis among school children. The R package implementing the proposed method is available at www.github.com/TianmengL/rectime.


Introduction
In various clinical and biomedical studies, the event of interest may happen multiple times, which is referred to as a recurrent event. Examples include recurrent bleedings in patients with hematologic malignancies (Stanworth et al., 2015) and recurrent cardiovascular events in subjects with diabetes (Van Der Heijden et al., 2013). During the follow-up of recurrent events, it is common to have repeated measurements of time-dependent covariates, and it is often of interest to investigate the effect of such covariates on the occurrence of recurrent events. Our motivation is from an observational study about pharyngitis among school children. Pharyngitis is often caused by viruses, but some bacteria, including streptococci, can cause pharyngitis as well. The goal of this study was to explore the effect of streptococci on the risk of pharyngitis. In this study, weekly visits were scheduled to monitor the recurrent occurrence of pharyngitis, and the status of streptococci infection was determined for those diagnosed with pharyngitis. In the meantime, monthly visits were scheduled for each participant to monitor the streptococci infection status regularly.
For the analysis of recurrent events, Prentice et al. (1981) and Andersen and Gill (1982) proposed the multiplicative model on the intensity function, which is interpreted as the instantaneous risk of event conditioning on the event history. To achieve a better interpretation and to allow flexible dependence structure among the recurrent events, various authors (Pepe and Cai, 1993;Lawless et al., 1997;Lin et al., 2000) discussed the regression models on mean or rate function of recurrent event process and assumed that the covariate effects were in a multiplicative fashion. Alternatively, the covariate effects may add to the baseline function. Liu and Wu (2011) considered the additive intensity model, and Schaubel et al. (2006) proposed the additive rates model where the covariate effects were in absolute forms. The multiplicative and additive models have different assumptions on the relationship between the covariate effects and the event process; thus it is desirable to estimate both types of effects under a general model setting. For univariate survival data where the event of interest only happens once, Lin and Ying (1995) proposed the additive-multiplicative model for the hazard function and Scheike and Zhang (2002) considered the Cox-Aalen model to allow the covariate effects to be time-varying. Recently, Cai et al. (2017a) studied modeling additive and multiplicative effects simultaneously in the mean residual life function. For recurrent event data, Han et al. (2016) proposed the additive-multiplicative models focusing on the markers contingent on recurrent events with an informative terminal event. Liu et al. (2010) considered the additive-multiplicative rates models for recurrent events, which allows the covariates to have both multiplicative and additive effects on the rate function of the recurrent event process.
Although the additive-multiplicative rates model allows the covariates to be time-dependent, it requires the observation of the complete history of the covariates, which is nearly impossible in practice. Typically, the time-dependent covariates are measured intermittently during the follow-up. One way to deal with the infrequently updated covariates is to predict the missing covariate values by smoothing the observed values, as was discussed in Raboud et al. (1993); Tsiatis et al. (1995); Boscardin et al. (1998); Bycott and Taylor (1998); Dafni and Tsiatis (1998) and summarized in Andersen and Liestøl (2003). Another commonly used approach is to jointly model the longitudinal covariate process and the event process. The joint models of repeated measured longitudinal data and time-to-event data have been studied extensively. Many authors considered the joint models of the two processes through latent random effects, including Wulfsohn and Tsiatis (1997); Xu and Zeger (2001);Vonesh et al. (2006), among others. In the setting of recurrent event data, Henderson et al. (2000) proposed to model the relationship between the time-dependent covariates and recurrent events by a latent Gaussian process, and Li (2016) considered the joint model of the recurrent event process and the binary covariate process. More complex models which considered the covariate process, recurrent event process and the terminal event simultaneously also have been studied, including Kim et al. (2012) and Cai et al. (2017b) among others. To our knowledge, little research has been done to explore the additive-multiplicative rates model with intermittently observed time-dependent covariates.
Recently, Li et al. (2016) and Cao et al. (2015) have proposed kernel smoothed estimators for the proportional rates model and hazard model, respectively, with intermittently observed time-dependent covariates. Specifically, Cao et al. (2015) considered the scenario where the timedependent covariates were not measured at event times but during regular visits only and Li et al. (2016) focused on the case where the covariates were measured at both recurrent times and regular visits. Lyu et al. (2021) and Sun et al. (2021a) considered kernel weighted estimation procedures for additive rates model and additive hazard regression model, respectively. In addition, Cao and Fine (2021) proposed a weighted last covariate carried forward approach for proportional hazard model with time-dependent covariates not observed at failure times and Sun et al. (2021b) proposed an estimation procedure based on inverse-rate-weighting and kernel-smoothing to estimate the proportional rates model with intermittently observed timedependent covariates measured at informative clinical visits. In this paper, we propose to extend the kernel smoothing method to the parameter estimation of the additive-multiplicative rates model. The proposed estimator is expected to be accurate; it relies on fewer assumptions of the underlying covariate and recurrent event processes.
The rest of the paper is organized as follows. The additive-multiplicative rates model and the proposed estimator are introduced in Section 2. Simulation studies to evaluate the performance of the proposed estimator are presented in Section 3. Section 4 includes the analysis of the pharyngitis data as introduced in the motivating example. Finally, a concluding remark is included in Section 5.

Additive-Multiplicative Rates Model
Suppose that n subjects are recruited in a study. Let i = 1, . . . , n index the subjects. Let N * i (t) denote the number of events that subject i has experienced at or prior to time t when there is no censoring. Let W i (t) = (Z i (t) , X i (t) ) be a p × 1 vector of possibly time-dependent covariates and let θ 0 = (γ 0 , β 0 ) be the corresponding true regression parameters. Following Liu et al. (2010), we assume the rate function of the counting process N * i (t) has the following form, where λ 0 (t) is an unspecified baseline rate function, and the link functions g and h are assumed to be known. Specifically, if we let g(x) = x and h(x) = exp(x), then model (1) becomes Therefore, the model can be regarded as a generalization of the semiparametric additive rates model and proportional rates model for recurrent event process. Let C i denote the censoring time of subject i and we assume that C i is independent of the counting process N * and denote by N i (t) = N * i (min(t, C i )) the observed number of events. Let [0, τ ] be a pre-specified time interval of interest, and the recurrent event process could potentially be observed beyond τ . For model estimation, we define the following process, where μ 0 (t) = t 0 λ 0 (u)du is the baseline mean function. Following Lin and Ying (1995), Liu et al. (2010) proposed the following estimating function for model (1): where D i (θ, t) is a p-dimensional smooth process involving W i (t) and θ, and According to Lin and Ying (1995) and Liu et al. (2010), a possible choice for D i (θ , t) is For a given θ, by solving n i=1 t 0 dM i (u, θ) = 0, the baseline mean function μ 0 (t) in (3) can be estimated by thus the estimating function in (3) becomes Since we are mostly interested in simultaneously modeling the relative and absolute difference in rate function due to the covariates, we focus on the model in (2) in the remainder of this paper. Under model (2), we have and thusD The estimating function in (5) can be written as Then the estimator θ can be obtained by solving U (θ) = 0.

Estimation When Covariates Are Intermittently Observed
When time-dependent covariates are present, the evaluation of the estimating function requires the time-dependent covariates to be observed continuously throughout the entire follow-up time. However, as is illustrated in the motivating example, time-dependent covariates were measured at intermittent visits, including the event times and monthly regular visits. Specifically, denote by O i (t) the number of measurements of the covariates at regular visits at or prior to time t, and we have dO i (t) = 1 when the ith subject has a regular visit at t. Moreover, we assume that the observation process . . , n}, which are assumed to be independent and identically distributed. Hence the values of covariates between measurement times are unknown, which renders that the estimating function in (6) cannot be evaluated based on the observed data only. A simple approach to deal with the intermittently observed timedependent covariates is to impute the missing values by carrying forward the last observed value. However, this approach imposes a strong assumption that the covariate processes are step functions and thus is expected to introduce bias in the model estimation in both survival and recurrent event analysis (Prentice, 1982;Faucett et al., 1998;Cao et al., 2015;Li et al., 2016;Lyu et al., 2021). Another possible approach is to impute the missing values between two observation times by linear interpolation. The linear interpolation method assumes that the covariate value is a linear function of time between every two adjacent observation times. However, this assumption may not hold in practice, especially for binary covariates. As an alternative approach, estimation methods based on kernel smoothing to deal with intermittently observed time-dependent covariates have gained growing interest recently (Cao et al., 2015;Li et al., 2016). In what follows, we propose a novel semiparametric estimator by applying the kernel smoothing method to deal with intermittently observed time-dependent covariates in the additive-multiplicative rates model. In contrast with methods that impute the missing values of each individual, the proposed method estimates the mean covariate processes via kernel smoothing and is expected to yield lower bias.
First, we show that the estimating function in (6) can be written as a function of empirical processes. Given a vector a, we define the operator ⊗ such that a ⊗0 = 1, a ⊗1 = a, a ⊗2 = aa . For k = 0, 1, we define Then the estimating function in (6) can be re-expressed as In equation (7), the quantities (u) can be easily calculated because the covariates are observed when recurrent events occur, as shown in the motivating example. However, calculating the ratios x (t, β) requires the covariate process W(·) to be observed throughout the follow-up period. Nevertheless, in real applications, it is usually not feasible to continuously monitor the covariates and the covariates are only measured intermittently.
Next, we present how to approximate these ratios based on the observed data by applying the kernel smoothing method. It can be easily seen that the processes S (k) With intermitently observed covariate data, we approximate the unknown ratios using kernel smoothing estimators that converge to the same target functions. Consider the following smoothed processes, where K h (t) = K(t/ h)/ h, K is a second order kernel function with a bounded support on [−1, 1], and h is the bandwidth parameter. For bandwidth selection, as is shown in the Appendix, the bandwidth h = cn −v , where 1/4 < v < 1/2, so we chose h = cn −1/3 where the constant c was chosen by following the cross-validation approach as presented in Lyu et al. (2021). We define m(t) = E{O i (t)}, then it can be shown that the kernel smoothed processes, , β), and S zx,h (t), converge in probability to s (k) z (t)m(t), s (k) x (t, β)m(t), s z2x (t, β)m(t), s zx (t)m(t), respectively. Since m(t) cancels out in the ratios, the ratios of the kernel smoothed processes converge to the ratios of the corresponding limiting processes. Thus, the ratios S (1) (7) can be replaced by the kernel smoothed counterparts. Moreover, to account for estimation bias near the boundary t = 0 due to positive observation times, we set x,h (u, β) The proposed estimator θ h = ( γ h , β h ) can be obtained by solving U (θ ) = 0. The large sample properties of θ h are summarized below. A detailed proof is included in the Appendix.
Theorem 1. Under regularity conditions 1-10 in the Appendix, θ h converges in probability to θ 0 . Moreover, as n → ∞, √ n( θ h − θ 0 ) converges in distribution to a zero mean normal random variable with variance The asymptotic variance of θ h involves unknown nuisance functions that need to be estimated using kernel smoothing methods. Hence bootstrap is recommended for variance estimation because of its better finite-sample performance.
For the estimation of the baseline mean function, we show that the estimator in (4) can be written as Thus, the baseline mean function can be estimated by Following Liu et al. (2010), to ensure that the estimated baseline mean function is monotone, an alternative estimator is μ 0,h (t, θ h ) = max 0 u t μ 0,h (u, θ h ).

Simulation
Simulation studies were conducted to evaluate the performance of the proposed method. 1000 data replicates with sample size 100 and 200 were generated for each simulation scenario. The resampling size in the bootstrap method for variance estimation was set to be 50. Since the rate model does not fully specify the probability feature, we simulated the recurrent events based on the following intensity model where u i is the frailty variable with mean 1 and variance σ 2 . The frailty variable induces the within-subject correlations. We explored two distributions of the frailty variable: gamma distribution and log-normal distribution, and two values of the variance σ 2 = 0.2, 0.4, for different levels of correlations. The baseline intensity function λ 0 (t) = 0.3I(t 10) + 0.5I(10 < t 20). Note that the intensity model in (9) implies the rate model λ{t|W i (t)} = γ 0 Z i (t) + exp{β 0 X i (t)}λ 0 (t).
To evaluate the proposed method on different types of covariates, we let Z i (t) be a continuous covariate and X i (t) be a binary covariate. Z i (t) was simulated by a linear function of time t as Z i (t) = b 0i + b 1i t. The random intercept b 0i was generated from a normal distribution with mean 0.5 and variance 0.05. The random slope was simulated from a normal distribution with mean −0.05 and variance 5 · 10 −4 . With a negative mean of the random slope, the covariate Z(·) has a decreasing time trend at the population level. For the binary covariate X i (t), we first generated the baseline X i (0) from a Bernoulli distribution with probability 0.2. Then the binary covariate process was assumed to alternate between states 0 and 1. We assumed that the duration of state 0 of subject i followed an exponential distribution with rate function 1/(ξ i g(t)) and the duration of state 1 followed an exponential distribution with rate 1/ξ i , where ξ i was a subject-specific random effect which followed a gamma distribution with mean 1 and variance 0.25. The value of g(t) was 4 for t ∈ [0, 10] and changed to 6 afterwards, which indicates a decreasing time trend at the population level. For the values of regression coefficients, we considered three scenarios: (1) the true model included both an additive part and a multiplicative part: β 0 = 0.5, γ 0 = 0.2; (2) the additive-multiplicative model degenerated to the additive rates model: β 0 = 0, γ 0 = 0.2; (3) the additive-multiplicative degenerated to the proportional rates model: β 0 = 0.5, γ 0 = 0. In all scenarios, we assumed that the covariates of a subject were always observed at the event times of the same subject. For the regular visits, we assumed that the covariates were measured at the baseline visit (time 0) and at each pre-scheduled regular visit per unit time interval, which means that there were 20 regular visits in time period (0, 20] of each subject. The time of the regular visit in each unit time interval was simulated from a uniform distribution from 0 to 1. The censoring time was randomly generated from a uniform distribution from 0 to 20. Both the event times and regular visits happened after the censoring time were dropped. We applied the proposed method to the simulated data and compared the results with those from the LCCF and linear interpolation approaches. We used the Epanechnikov kernel function in the proposed kernel smoothing method. The results for simulated datasets with gamma or lognormal frailty are presented in Tables 1 and 2, respectively. We provide the relative bias (Bias) and Monte Carlo standard deviation (SD) of the point estimations for each estimation method that we compared. For the proposed method, we also report the average of the estimated standard errors by bootstrap method (ASE) and the coverage percentage (CP). When the true model is the additive-multiplicative model (β 0 = 0.5, γ 0 = 0.2), the LCCF method gives biased estimations for both β and γ . The linear interpolation method has small bias for the estimation of γ , which is likely due to the linear feature of covariate Z(·), but gives biased estimations for β. The proposed method provides virtually unbiased estimations for both regression coefficients. The ASEs are close to the empirical SDs and the coverage percentages are around 95%. As the variance of the frailty increases from 0.2 to 0.4, the Monte Carlo SDs of the estimations increases. When the true model degenerates to the additive model or the multiplicative model with one significant covariate (β = 0, γ = 0.2; β = 0.5, γ = 0), the proposed method provides virtually unbiased estimations for both coefficients, which indicates the robustness of the additive-multiplicative model.

Real Data Analysis
In this section, we analyzed the Indian pharyngitis data (Jose et al., 2018) using the proposed method. Pharyngitis is the infection of the back of the throat and it can be caused by viruses or bacteria. When the cause is group A streptococcus (GAS), pharyngitis is also known as strep throat. The symptoms of GAS pharyngitis include sore throat, fever, nausea and it may cause some rare but serious diseases including rheumatic heart disease if left untreated. GAS pharyngitis is common in children from age 5 to age 15 and can be transmitted through saliva or nasal secretions. Meanwhile, other bacteria, for example, group G streptococcus (GGS), may cause pharyngitis with similar clinical symptoms as well. In the motivating example, 307 school children aged 7 to 11 years old in a rural area in Velore, India were recruited to investigate Table 1: Simulation results: the frailty followed a gamma distribution; n is the sample size; σ 2 is the variance of the frailty distribution;m is the average number of recurrent events; Bias is the relative bias computed by dividing the difference of the mean of the 1000 estimated parameters and the true value by the true value (if the true value is 0, Bias is the mean of the 1000 estimated parameters); SD is the standard deviation of the 1000 estimated values; ASE is the mean of the 1000 estimated standard errors by bootstrap method; CP is the proportion of 95% confidence intervals covering the true value. the relationship between streptococcal infections and the risk of pharyngitis. Each child was examined weekly for the symptoms of pharyngitis. For those who were diagnosed with pharyngitis, throat cultures were obtained to test if GAS and GGS were positive. In the meantime, to monitor the streptococci status regularly, monthly regular visits were scheduled for each child.
Since the regular visits were pre-scheduled, it is reasonable to assume that they are independent of the recurrent event, censoring and covariate processes, as is required by the proposed method. Table 2: Simulation results: the frailty followed a lognormal distribution; n is the sample size; σ 2 is the variance of the frailty distribution;m is the average number of recurrent events; Bias is the relative bias computed by dividing the difference of the mean of the 1000 estimated parameters and the true value by the true value (if the true value is 0, Bias is the mean of the 1000 estimated parameters); SD is the standard deviation of the 1000 estimated values; ASE is the mean of the 1000 estimated standard errors by bootstrap method; CP is the proportion of 95% confidence intervals covering the true value. We considered the following four candidate models: (1) both covariates have multiplicative effects (Model MM); (2) both covariates have additive effects (Model AA); (3) GAS has an additive effect and GGS has a multiplicative effect (Model AM); (4) GAS has a multiplicative effect and GGS has an additive effect (Model MA). The same approach as described in Section 2 was applied to select the bandwidth used in the proposed kernel smoothing method. Table 3 shows the regression results of the four models. All four models suggest that the presence of Table 3: Analysis of Indian pharyngitis data. AM is the additive-multiplicative rates model which includes GAS in the additive part and GGS in the multiplicative part; MA is the additivemultiplicative rates model which includes GAS in the multiplicative part and GGS in the additive part; MM is the proportional rates model; AA is the additive rates model. Est is the estimated regression coefficient; SE is the standard error estimated by bootstrap with resampling size 100; CI is the 95% confidence interval. GAS increases the risk of pharyngitis, while the effect of GGS is not significant in any model. In light of the same direction of the effect and the similar statistical significance from the four candidate models, the choice would be based on study objectives and the interpretation of the effect. As discussed in Schaubel et al. (2006), in certain settings, the absolute covariate effect is of more interest than the relative covariate effect. For instance, the former can directly provide information for predicting change in event rate attributable to a covariate, while the latter would need information on the baseline rate.

Discussion
In this paper, we proposed a semiparametric estimator for the regression coefficients of the additive-multiplicative rates model to deal with the intermittently observed time-dependent covariates. The additive-multiplicative rates model generalizes the proportional rates model and additive rates model, and hence allows some covariates to have multiplicative effects on the risk of recurrent events and others to have additive effects. The proposed method applies the nonparametric kernel smoothing approach to estimate the mean processes of the timedependent covariates, thus it does not rely on any assumption of the covariate distribution or any specification of the covariate process and is expected to be more robust. The proposed method requires that the rate function m(t) of the observation time process is positive and bounded on [0, τ ]. As the observations of time-dependent covariates become more dense, we expect better performance of the proposed method.
In this paper, we assume that the observation process is independent of the covariates, the recurrent event process, and the censoring time. When such an independence assumption is violated, the proposed method may yield biased estimation. When the rate function of the observation process is determined by observed covariates, one can assume a proportional rate model on the observation process and construct an inverse-rate-weighting in the kernel smoothing estimator for the outcome model, following the arguments of Sun et al. (2021b). Extending their method to additive-multiplicative rates model is a future research direction.
In practice, a problem is to determine the covariates included in the additive part Z(t) and the multiplicative part X(t). As discussed in Liu et al. (2010), if a covariate is expected to greatly influence the risk difference or the researchers are most interested in the absolute risks, then it should be included in Z(t). Otherwise, if a covariate is expected to strongly influence the risk ratios or the researchers are interested in the relative risks, then it should be included in X(t). If the underlying biological process is not clear and the number of covariates is small, an alternative way is to consider all the possible candidate models and then develop rigorous model selection approach to determine the best model. The mean-square-type distance measure between the observed and expected recurrences implemented in Liu et al. (2010) could be a model selection criteria for the additive and/or multiplicative model structure. However, it cannot be directly applied to data with intermittently observed time-dependent covariates. Future research on model selection procedures for intermittently observed time-dependent covariates is warranted.
An R package has been developed to implement the proposed estimator for the additivemultiplicative rates model with intermittently observed time-dependent covariates, as well as the estimators for the proportional rates model and the additive rates model with such covariates. The R package rectime is available at www.github.com/TianmengL/rectime.

Proof of Consistency
To show the consistency of θ h , it is sufficient to prove that the processes that constitute the estimating function U (θ), including n −1 n i=1 τ 0 Z i (u) exp{−β X i (u)}dN i (u), n −1 n i=1 τ 0 X i (u)dN i (u), n −1 n i=1 dN i (u), S (k) z,h (t), S (k) x,h (t, β), k = 0, 1, S z2x,h (t, β) and S zx,h (t), converge to their limits uniformly. We know θ 0 = (γ 0 , β 0 ) where γ 0 is a m × 1 vector, β 0 is a q × 1 vector and m + q = p. Since θ 0 is in a compact set in R p by assumption 3, β 0 is contained in a compact set B in R q . The function classes F z,k = { According to Theorem 2.14.9 in Van Der Vaart and Wellner (1996) and following similar steps in the Appendix of Li et al. (2016)