Fitting semiparametric regression model to recurrent event data

Sy Han (Steven) Chiou

2019-11-30

Source: vignettes/reReg-reg.Rmd

In this vignette, we demonstrate how to use the reReg function in reReg package to fit semiparametric regression model to recurrent event data.

> library(reReg)
> args(reReg)
function (formula, data, B = 200, method = c("cox.LWYY", "cox.GL", 
    "cox.HW", "am.GL", "am.XCHWY", "sc.XCYH"), se = c("NULL", 
    "bootstrap", "resampling"), contrasts = NULL, control = list()) 
NULL

The arguments are as follows

  • 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 reSurv. See the vignette on Visualization of recurrent event data for examples in creating reSurv objects.
  • data an optional data frame in which to interpret the variables occurring in the formula.
  • B a numeric value specifies the number of resampling for variance estimation. When B = 0, variance estimation will not be performed.
  • method a character string specifying the underlying model.
  • se a character string specifying the method for standard error estimation.
  • contrasts an optional list.
  • control a list of control parameters.

We will illustrate the usage of reReg with simulated data generated from simSC. Readers are referred to the vignette on Simulating recurrent event data for using simSC to generate recurrent event data.

Point estimation

We give brief description on each of the available methods below.

D. Y. Lin, Wei, and Ying (1998) cox.LWYY 

> set.seed(123)
> dat.cox <- simSC(n = 500, a = c(1, -1), b = c(1, -1))
> fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox)
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox)

Method: Lin-Wei-Yang-Ying method (fitted with coxph with robust variance)

Coefficients effect:
   Estimate StdErr z.value   p.value    
x1    0.965  0.031  31.109 < 2.2e-16 ***
x2   -1.004  0.018 -57.072 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ghosh and Lin (2002) cox.GL 

> fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, method = "cox.GL")
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
    method = "cox.GL")

Coefficients (rate):
   Estimate StdErr z.value p.value
x1    0.974     NA      NA      NA
x2   -0.963     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    1.009     NA      NA      NA
x2   -0.951     NA      NA      NA

Huang and Wang (2004) cox.HW 

> fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, method = "cox.HW")
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
    method = "cox.HW")

Method: Huang-Wang Model 

Coefficients (rate):
   Estimate StdErr z.value p.value
x1    0.975     NA      NA      NA
x2   -1.016     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    1.021     NA      NA      NA
x2   -0.962     NA      NA      NA

Ghosh and Lin (2003) am.GL 

> set.seed(123)
> dat.am <- simSC(n = 500, a = c(1, -1), b = c(1, -1), type = "am")
> fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.am, method = "am.GL")
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.am, 
    method = "am.GL")

Method: Ghosh-Lin Model 

Coefficients (rate):
   Estimate StdErr z.value p.value
x1    1.047     NA      NA      NA
x2   -0.965     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.825     NA      NA      NA
x2   -0.798     NA      NA      NA

Xu et al. (2017) am.XCHWY 

> fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.am, method = "am.XCHWY")
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.am, 
    method = "am.XCHWY")

Method: Xu et al. (2016) Model 

Coefficients (rate):
   Estimate StdErr z.value p.value
x1    1.196     NA      NA      NA
x2   -0.845     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    1.006     NA      NA      NA
x2   -0.605     NA      NA      NA

Xu et al. (2018) sc.XCYH

Fitting recurrent event data generated from the Cox model:

> fit1 <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, method = "sc.XCYH")
> summary(fit1)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
    method = "sc.XCYH")

Method: Generalized Scale-Change Model 

Scale-change effect:
   Estimate StdErr z.value p.value
x1   -0.180     NA      NA      NA
x2   -0.029     NA      NA      NA

Multiplicative coefficients:
   Estimate StdErr z.value p.value
x1    0.875     NA      NA      NA
x2   -1.011     NA      NA      NA

Fitting recurrent event data generated from the accelerated model:

> fit2 <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.am, method = "sc.XCYH")
> summary(fit2)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.am, 
    method = "sc.XCYH")

Method: Generalized Scale-Change Model 

Scale-change effect:
   Estimate StdErr z.value p.value
x1    0.884     NA      NA      NA
x2   -0.971     NA      NA      NA

Multiplicative coefficients:
   Estimate StdErr z.value p.value
x1    1.031     NA      NA      NA
x2   -0.898     NA      NA      NA

Variance estimation

The variance estimation was not performed in the examples above. However, the package provides several approaches for variance estimation. We will take cox.HW for example.

Bootstrap

> system.time(fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
+     method = "cox.HW", se = "bootstrap", B = 100))
   user  system elapsed 
265.813   0.032 265.963 
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
    B = 100, method = "cox.HW", se = "bootstrap")

Method: Huang-Wang Model 

Coefficients (rate):
   Estimate StdErr z.value   p.value    
x1    0.975  0.054  18.134 < 2.2e-16 ***
x2   -1.016  0.024 -41.690 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (hazard):
   Estimate StdErr z.value   p.value    
x1    1.021  0.142   7.171 < 2.2e-16 ***
x2   -0.962  0.066 -14.598 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bootstrap with parallel computing

Detect the number of CPU cores on the current host.

> (core <- parallel::detectCores())
[1] 16

The default number of CPU cores used in the parallel computing is half of the available CPU cores. We will use all the CPU cores available in this example.

> system.time(fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
+     method = "cox.HW", se = "bootstrap", B = 100, control = list(parallel = TRUE, 
+         parCl = core)))
   user  system elapsed 
 81.544   0.373 113.792 
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
    B = 100, method = "cox.HW", se = "bootstrap", control = list(parallel = TRUE, 
        parCl = core))

Method: Huang-Wang Model 

Coefficients (rate):
   Estimate StdErr z.value   p.value    
x1    0.975  0.049  19.872 < 2.2e-16 ***
x2   -1.016  0.026 -38.700 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (hazard):
   Estimate StdErr z.value   p.value    
x1    1.021  0.131   7.785 < 2.2e-16 ***
x2   -0.962  0.077 -12.497 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sandwich (resampling) variance estimator

> system.time(fit <- reReg(Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
+     method = "cox.HW", se = "resampling", B = 100))
   user  system elapsed 
111.568   0.172 111.781 
> summary(fit)
Call: reReg(formula = Recur(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
    B = 100, method = "cox.HW", se = "resampling")

Method: Huang-Wang Model 

Coefficients (rate):
   Estimate StdErr z.value   p.value    
x1    0.975  0.040  24.133 < 2.2e-16 ***
x2   -1.016  0.027 -37.199 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (hazard):
   Estimate StdErr z.value   p.value    
x1    1.021  0.156   6.541 < 2.2e-16 ***
x2   -0.962  0.088 -10.970 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Reference

Ghosh, Debashis, and D Y Lin. 2003. “Semiparametric Analysis of Recurrent Events Data in the Presence of Dependent Censoring.” Biometrics 59 (4): 877–85.

Ghosh, Debashis, and DY Lin. 2002. “Marginal Regression Models for Recurrent and Terminal Events.” Statistica Sinica 12 (3): 663–88.

Huang, Chiung-Yu, and Mei-Cheng Wang. 2004. “Joint Modeling and Estimation for Recurrent Event Processes and Failure Time Data.” Journal of the American Statistical Association 99 (468): 1153–65.

Lin, D Y, L J Wei, and Zhiliang Ying. 1998. “Accelerated Failure Time Models for Counting Processes.” Biometrika 85 (3): 605–18.

Xu, Gongjun, Sy Han Chiou, Chiung-Yu Huang, Mei-Cheng Wang, and Jun Yan. 2017. “Joint Scale-Change Models for Recurrent Events and Failure Time.” Journal of the American Statistical Association 112 (518). Taylor & Francis: 794–805.

Xu, Gongjun, Sy Han Chiou, Jun Yan, Kieren Marr, and Chiung-Yu Huang. 2018. “Generalized Scale-Change Models for Recurrent Event Processes Under Informative Censoring.” Statistica Sinica.