Fitting semiparametric regression model to recurrent event data

Sy Han (Steven) Chiou

2018-11-01

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(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.cox)
> summary(fit)
Call: reReg(formula = reSurv(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.960  0.032  30.025 < 2.2e-16 ***
x2   -1.010  0.022 -45.580 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ghosh and Lin (2002) cox.GL 

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

Coefficients (rate):
   Estimate StdErr z.value p.value
x1    0.905     NA      NA      NA
x2   -0.973     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    1.193     NA      NA      NA
x2   -0.993     NA      NA      NA

Huang and Wang (2004) cox.HW 

> fit <- reReg(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.cox, method = "cox.HW")
> summary(fit)
Call: reReg(formula = reSurv(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.937     NA      NA      NA
x2   -1.077     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    1.168     NA      NA      NA
x2   -1.064     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(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.am, method = "am.GL")
> summary(fit)
Call: reReg(formula = reSurv(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    0.906     NA      NA      NA
x2   -0.912     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.496     NA      NA      NA
x2   -1.556     NA      NA      NA

Xu et al. (2017) am.XCHWY 

> fit <- reReg(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.am, method = "am.XCHWY")
> summary(fit)
Call: reReg(formula = reSurv(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.075     NA      NA      NA
x2   -0.844     NA      NA      NA

Coefficients (hazard):
   Estimate StdErr z.value p.value
x1    0.754     NA      NA      NA
x2   -1.475     NA      NA      NA

Xu et al. (2018) sc.XCYH

Fitting recurrent event data generated from the Cox model:

> fit1 <- reReg(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.cox, method = "sc.XCYH")
> summary(fit1)
Call: reReg(formula = reSurv(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.120     NA      NA      NA
x2   -0.003     NA      NA      NA

Multiplicative coefficients:
   Estimate StdErr z.value p.value
x1    0.966     NA      NA      NA
x2   -1.052     NA      NA      NA

Fitting recurrent event data generated from the accelerated model:

> fit2 <- reReg(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.am, method = "sc.XCYH")
> summary(fit2)
Call: reReg(formula = reSurv(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.917     NA      NA      NA
x2   -0.979     NA      NA      NA

Multiplicative coefficients:
   Estimate StdErr z.value p.value
x1    0.988     NA      NA      NA
x2   -0.911     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(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
+     method = "cox.HW", se = "bootstrap", B = 100))
   user  system elapsed 
528.482   0.096 528.805 
> summary(fit)
Call: reReg(formula = reSurv(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.937  0.066  14.193 < 2.2e-16 ***
x2   -1.077  0.057 -18.821 < 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.168  0.158   7.412 < 2.2e-16 ***
x2   -1.064  0.114  -9.314 < 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(reSurv(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 
 68.582   0.479 131.577 
> summary(fit)
Call: reReg(formula = reSurv(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.937  0.075  12.543 < 2.2e-16 ***
x2   -1.077  0.065 -16.578 < 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.168  0.152   7.692 < 2.2e-16 ***
x2   -1.064  0.109  -9.743 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sandwich (resampling) variance estimator

> system.time(fit <- reReg(reSurv(Time, id, event, status) ~ x1 + x2, data = dat.cox, 
+     method = "cox.HW", se = "resampling", B = 100))
   user  system elapsed 
100.790   0.000 100.813 
> summary(fit)
Call: reReg(formula = reSurv(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.937  0.051  18.391 < 2.2e-16 ***
x2   -1.077  0.051 -21.274 < 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.168  0.124   9.386 < 2.2e-16 ***
x2   -1.064  0.092 -11.540 < 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.