Fits a general (joint) semiparametric regression model for the recurrent event data, where the rate function of the underlying recurrent event process and the hazard function of the terminal event can be specified as a Cox-type model, an accelerated mean model, an accelerated rate model, or a generalized scale-change model. See details for model specifications.

reReg(
formula,
data,
subset,
model = "cox",
B = 0,
se = c("boot", "sand"),
control = list()
)

## Arguments

formula a formula object, with the response on the left of a "~" operator, and the predictors on the right. The response must be a recurrent event survival object as returned by function Recur. an optional data frame in which to interpret the variables occurring in the "formula". n optional logical vector specifying a subset of observations to be used in the fitting process. a character string specifying the underlying model. See Details. a numeric value specifies the number of bootstraps for variance estimation. When B = 0, variance estimation will not be performed. a character string specifying the method for the variance estimation. See Details. bootnonparametric bootstrap approach sandresampling-based sandwich estimator a list of control parameters.

## Details

Suppose the recurrent event process and the failure events are observed in the time interval $$t\in[0,\tau]$$, for some constant $$\tau$$. We formulate the recurrent event rate function, $$\lambda(t)$$, and the terminal event hazard function, $$h(t)$$, in the form of $$\lambda(t) = Z \lambda_0(te^{X^\top\alpha}) e^{X^\top\beta}, h(t) = Z h_0(te^{X^\top\eta})e^{X^\top\theta},$$ where $$\lambda_0(t)$$ is the baseline rate function, $$h_0(t)$$ is the baseline hazard function, $$X$$ is a $$n$$ by $$p$$ covariate matrix and $$\alpha$$, $$Z$$ is an unobserved shared frailty variable, and $$(\alpha, \eta)$$ and $$(\beta, \theta)$$ correspond to the shape and size parameters, respectively. The model includes several popular semiparametric models as special cases, which can be specified via the model argument with the rate function and hazard function separated by "|". For examples, Wang, Qin and Chiang (2001) ($$\alpha = \eta = \theta = 0$$) can be called with model = "cox"; Huang and Wang (2004) ($$\alpha = \eta = 0$$) can be called with model = "cox|cox"; Xu et al. (2017) ($$\alpha = \beta$$ and $$\eta = \theta$$) can be called with model = "am|am"; Xu et al. (2019) ($$\eta = \theta = 0$$) can be called with model = "gsc". Users can mix the models depending on the application. For example, model = "cox|ar" postulate a Cox proportional model for the recurrent event rate function and an accelerated rate model for the terminal event hazard function ($$\alpha = \theta = 0$$). If only one model is specified without an "|", it is used for both the rate function and the hazard function. For example, specifying model = "cox" is equivalent to model = "cox|cox". Some models that assumes Z = 1 and requires independent censoring are also implemented in reReg; these includes model = "cox.LWYY" for Lin et al. (2000), model = "cox.GL" for Ghosh and Lin (2002), and model = "am.GL" for Ghosh and Lin (2003).

The available methods for variance estimation are:

boot

performs nonparametric bootstrap.

sand

performs the efficient resampling-based variance estimation.

The control list consists of the following parameters:

tol

absolute error tolerance.

alpha, beta, eta, theta

initial guesses used for root search.

solver

the equation solver used for root search. The available options are BB::BBsolve, BB::dfsane, BB:BBoptim, and optim.

eqType

a character string indicating whether the log-rank type estimating equation or the Gehan-type estimating equation (when available) will be used.

boot.parallel

an logical value indicating whether parallel computation will be applied when se = "boot" is called.

boot.parCl

an integer value specifying the number of CPU cores to be used when parallel = TRUE. The default value is half the CPU cores on the current host.

Lin, D., Wei, L., Yang, I. and Ying, Z. (2000). Semiparametric Regression for the Mean and Rate Functions of Recurrent Events. Journal of the Royal Statistical Society: Series B (Methodological), 62: 711--730.

Wang, M.-C., Qin, J., and Chiang, C.-T. (2001). Analyzing Recurrent Event Data with Informative Censoring. Journal of the American Statistical Association, 96(455): 1057--1065.

Ghosh, D. and Lin, D.Y. (2002). Marginal Regression Models for Recurrent and Terminal Events. Statistica Sinica: 663--688.

Ghosh, D. and Lin, D.Y. (2003). Semiparametric Analysis of Recurrent Events Data in the Presence of Dependent Censoring. Biometrics, 59: 877--885.

Huang, C.-Y. and Wang, M.-C. (2004). Joint Modeling and Estimation for Recurrent Event Processes and Failure Time Data. Journal of the American Statistical Association, 99(468): 1153--1165.

Xu, G., Chiou, S.H., Huang, C.-Y., Wang, M.-C. and Yan, J. (2017). Joint Scale-change Models for Recurrent Events and Failure Time. Journal of the American Statistical Association, 112(518): 796--805.

Xu, G., Chiou, S.H.,Yan, J., Marr, K., and Huang, C.-Y. (2019). Generalized Scale-Change Models for Recurrent Event Processes under Informative Censoring. Statistica Sinica: pre-print.

## See also

Recur, simGSC

## Examples

data(simDat)

## Nonparametric estimate
plot(reReg(Recur(t.start %to% t.stop, id, event, status) ~ 1, data = simDat, B = 50)) fm <- Recur(t.start %to% t.stop, id, event, status) ~ x1 + x2
## Fit the Cox rate model
summary(reReg(fm, data = simDat, model = "cox", B = 50))
#> Call:
#> reReg(formula = fm, data = simDat, model = "cox", B = 50)
#>
#> Recurrent event process:
#>    Estimate   StdErr  z.value p.value
#> x1 -1.00483  0.16506 -6.08771       0
#> x2 -0.97517  0.14748 -6.61233       0
#> ## Fit the joint Cox/Cox model
summary(reReg(fm, data = simDat, model = "cox|cox", B = 50))
#> Call:
#> reReg(formula = fm, data = simDat, model = "cox|cox", B = 50)
#>
#> Recurrent event process:
#>    Estimate   StdErr  z.value p.value
#> x1 -1.00483  0.14951 -6.72088       0
#> x2 -0.97517  0.13883 -7.02440       0
#>
#> Terminal event:
#>    Estimate  StdErr z.value p.value
#> x1  1.05295 0.26785 3.93119  0.0001
#> x2  0.85086 0.26255 3.24076  0.0012
#> ## Fit the scale-change rate model
summary(reReg(fm, data = simDat, model = "gsc", B = 50, se = "sand"))
#> Call:
#> reReg(formula = fm, data = simDat, model = "gsc", B = 50, se = "sand")
#>
#> Recurrent event process (shape):
#>    Estimate   StdErr z.value p.value
#> x1 -0.02252  0.38372 -0.0587  0.9532
#> x2 -0.11480  0.24091 -0.4765  0.6337
#>
#> Recurrent event process (size):
#>    Estimate   StdErr z.value p.value
#> x1 -1.02816  0.29776 -3.4529 0.00055 ***
#> x2 -1.04991  0.25580 -4.1044   4e-05 ***
#>