Fitting semiparametric regression model to recurrent event data

Sy Han (Steven) Chiou

2020-06-11

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)
Registered S3 methods overwritten by 'ggplot2':
  method         from 
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang
> 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.814  0.067  12.129 < 2.2e-16 ***
x2   -0.861  0.040 -21.632 < 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.531     NA      NA      NA
x2   -0.795     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.772     NA      NA      NA
x2   -0.867     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.874     NA      NA      NA
x2   -0.951     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.831     NA      NA      NA
x2   -0.951     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.386     NA      NA      NA
x2   -1.040     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.729     NA      NA      NA
x2   -1.406     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.122     NA      NA      NA
x2   -1.040     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.531     NA      NA      NA
x2   -1.369     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.044     NA      NA      NA
x2   -0.015     NA      NA      NA

Multiplicative coefficients:
   Estimate StdErr z.value p.value
x1    0.878     NA      NA      NA
x2   -0.944     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    1.015     NA      NA      NA
x2   -0.899     NA      NA      NA

Multiplicative coefficients:
   Estimate StdErr z.value p.value
x1    1.078     NA      NA      NA
x2   -0.965     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 
277.988   0.180 278.316 
> 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.874  0.091   9.639 < 2.2e-16 ***
x2   -0.951  0.045 -21.224 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (hazard):
   Estimate StdErr z.value   p.value    
x1    0.831  0.161   5.176 < 2.2e-16 ***
x2   -0.951  0.091 -10.420 < 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] 8

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 
 47.932   0.132 103.245 
> 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.874  0.084  10.387 < 2.2e-16 ***
x2   -0.951  0.043 -21.923 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (hazard):
   Estimate StdErr z.value   p.value    
x1    0.831  0.146   5.699 < 2.2e-16 ***
x2   -0.951  0.093 -10.172 < 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 
 66.460   0.024  66.494 
> 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.874  0.062  14.019 < 2.2e-16 ***
x2   -0.951  0.038 -25.057 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficients (hazard):
   Estimate StdErr z.value   p.value    
x1    0.831  0.148   5.595 < 2.2e-16 ***
x2   -0.951  0.080 -11.874 < 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.