The Cox-type proportional rate model (available in reReg() via model = "cox.LWYY") is a common semiparametric regression model for recurrent event processes under the noninformative censoring assumption. Since the construction of the pseudo-partial score function ignores the dependency among recurrent events, the efficiency of the Cox-type proportional rate model can be improved by combining a system of weighted pseudo-partial score equations via the generalized method of moments (GMM) and empirical likelihood (EL) estimation (Huang and Huang 2023). In this vignette, we demonstrate the proposed GMM and EL procedures in the reReg package. We will illustrate the idea with simulated data generated from the simGSC() function.

> library(reReg)
> packageVersion("reReg")
[1] '1.4.6'

Loading data set

We first generate a simulated data where the recurrent event process is generated from a Cox-type proportional rate model. This can be achieved with the default settings of simGSC(); see the vignette on simulating recurrent event data for illustrations of simGSC().

> set.seed(1); head(dat <- simGSC(50, frailty = rep(1, 50)))
  id   t.start    t.stop event status x1          x2
1  1 0.0000000 0.1363067     1      0  0 -0.02806437
2  1 0.1363067 1.0900133     1      0  0 -0.02806437
3  1 1.0900133 1.5131840     1      0  0 -0.02806437
4  1 1.5131840 2.0214052     1      0  0 -0.02806437
5  1 2.0214052 2.1870742     1      0  0 -0.02806437
6  1 2.1870742 3.9498226     1      0  0 -0.02806437

Fitting Cox-type model

We can fit the Cox-type proportional rate model (Lin et al. 2000) with reReg() via the following specification.

> fit <- reReg(Recur(t.stop, id, event) ~ x1 + x2, model = "cox.LWYY", data = dat)
> summary(fit)
Call: 
reReg(formula = Recur(t.stop, id, event) ~ x1 + x2, data = dat, 
    model = "cox.LWYY")

Fitted with the Cox model of Lin et al. (2000):
Recurrent event process:
   Estimate   StdErr z.value   p.value    
x1 -0.99571  0.12516 -7.9555 1.785e-15 ***
x2 -1.08198  0.12847 -8.4221 < 2.2e-16 ***

Huang and Huang (2023) proposed improved semiparametric estimation methods of the above proportional rate model. By combining different weighted pseudo-partial score functions through generalized methods of moments or empirical likelihood methods, Huang and Huang (2023) showed that substantial efficiency gain can be achieved without imposing additional model assumptions than those assumed in Lin et al. (2000). Specifically, the efficiency gain is a result of combining a set of asymptotically independently and identically distributed estimating equations derived from the weighted pseudo-partial score equations of Lin et al. (2000).

Improved estimation

Two available approaches to combine estimating equations are the generalized method of moments (GMM) and the empirical likelihood (EL) estimations, both of which can be specified with the cppl argument within the control list. In addition to the combination approach, weight functions to be combined with the weighted pseudo-partial score functions can also be specified with the cppl.wfun argument within the control list. In theory, combining more weight functions leads to better asymptotic efficiency, especially when the sample size is large. However, as observed in simulation studies of Huang and Huang (2023), combining one or two weighted functions can achieve substantial efficiency gain under scenarios with small to moderate sample sizes. When the sample size is small, Huang and Huang (2023) suggests using a combination of the unit weight and a decreasing function, since the decreasing function can downweight the later time period, in which the events are more effectively censored. For those reason, the reReg() currently only allows user to specify up to two weight functions with cppl.wfun. The reReg() allows users choose either the cumulative baseline rate function or the Gehan’s weight as two possible options for the weight functions. Additionally, the reReg() also allows users to specify arbitrary function in cppl.wfun via function formulas.

The following gives an example of the improved estimation by combining estimating equations via GMM with the cumulative baseline rate function and the Gehan’s weight.

> fit2 <- reReg(Recur(t.stop, id, event) ~ x1 + x2, model = "cox.LWYY", data = dat, 
+               control = list(cppl = "GMM", cppl.wfun = list("Gehan", "cumbase"))) 
> summary(fit2) 
Call: 
reReg(formula = Recur(t.stop, id, event) ~ x1 + x2, data = dat, 
    model = "cox.LWYY", control = list(cppl = "GMM", cppl.wfun = list("Gehan", 
        "cumbase")))

Recurrent event process:
   Estimate   StdErr z.value   p.value    
x1 -0.99767  0.12536 -7.9582 1.746e-15 ***
x2 -1.07355  0.12585 -8.5305 < 2.2e-16 ***

For this data set, the standard errors changed from 0.1251605 and 0.1284696 to 0.1253645 and 0.1258481, respectively. This achieves a -0.16% and 2.04% efficiency gain, respectively. The above result can be called by updating the original fit, fit, with the update() function, as illustrated in the following.

> fit2 <- update(fit, control = list(cppl = "GMM", cppl.wfun = list("Gehan", "cumbase")))

The following gives an example of the improved estimation by combining estimating equations via EL with the user-specified function \(1 / (x + 1)\) and the Gehan’s weight.

> fit3 <- reReg(Recur(t.stop, id, event) ~ x1 + x2, model = "cox.LWYY", data = dat, 
+               control = list(cppl = "EL", cppl.wfun = list(function(x) 1 / (x + 1), "Gehan"))) 
> summary(fit3) 
Call: 
reReg(formula = Recur(t.stop, id, event) ~ x1 + x2, data = dat, 
    model = "cox.LWYY", control = list(cppl = "EL", cppl.wfun = list(function(x) 1/(x + 
        1), "Gehan")))

Recurrent event process:
   Estimate   StdErr z.value   p.value    
x1 -0.94918  0.12057 -7.8725 3.475e-15 ***
x2 -1.04422  0.12208 -8.5536 < 2.2e-16 ***

For this data set, the standard errors changed from 0.1251605 and 0.1284696 to 0.1205686 and 0.1220787, respectively. This achieves a 3.67% and 4.97% efficiency gain, respectively. Although the above example assumes the covariates are time-independent, the proposed approach and the implementation allow time-dependent covariates.

Conclusion

Although this vignette is focus on improving the proportional rate model, the proposed approaches by Huang and Huang (2023) can be extended to handle other types of rate function. Such an extension will be investigated and the reReg package will be expanded accordingly.

Reference

Huang, Ming-Yueh, and Chiung-Yu Huang. 2023. “Improved Semiparametric Estimation of the Proportional Rate Model with Recurrent Event Data.” Biometrics 79 (3): 1686–1700.

Lin, D. Y., L. J. Wei, I. Yang, and Z. Ying. 2000. “Semiparametric Regression for the Mean and Rate Functions of Recurrent Events.” Journal of the Royal Statistical Society: Series B 62 (4): 711–30.