Movatterモバイル変換

reReg: Recurrent Event Regression

Description

The package offers a comprehensive collection of practical and easy-to-use tools for analyzingrecurrent event data, with or without the presence of a (possibly) correlated terminal event.The modeling framework is based on a joint frailty scale-change model,that encompasses many existing models, including the popular Cox-type models,as special cases and accommodates informative censoring through a subject-specific frailty.The implemented estimating procedure does not require any parametric assumption on the frailtydistribution.The package allows the users to specify different model forms for both the recurrent event processand the terminal event.The package also includes tools for visualization of recurrent events andsimulation from the regression models.

Author(s)

Maintainer: Sy Han (Steven) Chiouschiou@smu.edu

Authors:

• Chiung-Yu HuangChiungYu.Huang@ucsf.edu

References

Chiou, S.H., Xu, G., Yan, J. and Huang, C.-Y. (2023). Regression Modeling for Recurrent Events Possibly with an Informative Terminal Event Using R Package reReg.Journal of Statistical Software,105(5): 1–34.

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 EventProcesses under Informative Censoring.Statistica Sinica,30: 1773–1795.

See Also

Useful links:

• https://github.com/stc04003/reReg

• Report bugs athttps://github.com/stc04003/reReg/issues

Function used to combine baseline functions in one plot

Description

Combine different plots into one.

Usage

basebind(..., legend.title, legend.labels, control = list())

Arguments

Examples

data(simDat)fm <- Recur(t.stop, id, event, status) ~ x1 + x2fit1 <- reReg(fm, subset = x1 == 0, data = simDat, B = 200)fit2 <- reReg(fm, subset = x1 == 1, data = simDat, B = 200)basebind(plot(fit1), plot(fit2))

Produce Event Plot or Mean Cumulative Function Plot

Description

Plot the event plot or the mean cumulative function (MCF) from anRecur object.

Usage

## S3 method for class 'Recur'plot(  x,  mcf = FALSE,  event.result = c("increasing", "decreasing", "asis"),  event.calendarTime = FALSE,  mcf.adjustRiskset = TRUE,  mcf.conf.int = FALSE,  control = list(),  ...)

Arguments

Details

The argumentcontrol consists of options with argument defaults to a list withthe following values:

xlab

customizable x-label, default value is "Time".

ylab

customizable y-label, default value is "Subject" for event plot and"Cumulative mean" for MCF plot.

main

customizable title, the default value is "Recurrent event plot"whenmcf = FALSE and"Sample cumulative mean function plot" whenmcf = TRUE.

terminal.name

customizable label for terminal event,the default value is "Terminal event".

recurrent.name

customizable legend title for recurrent event,the default value is "Recurrent events".

recurrent.types

customizable label for recurrent event type,the default value isNULL.

alpha

between 0 and 1, controls the transparency of points.

Thexlab,ylab andmain parameters can be specifiedoutside of thecontrol list.

Value

Aggplot object.

References

Nelson, W. B. (1995) Confidence Limits for Recurrence Data-Applied to Costor Number of Product Repairs.Technometrics,37(2): 147–157.

See Also

Recur,plotEvents,mcf

Examples

data(simDat)reObj <- with(simDat, Recur(t.start %to% t.stop, id, event, status))## Event plots:plot(reObj)plot(reObj, event.result = "decreasing")## With (hypothetical) multiple event typessimDat$event2 <- with(simDat, ifelse(t.stop > 10 & event > 0, 2, event))reObj2 <- with(simDat, Recur(t.start %to% t.stop, id, event2, status))plot(reObj2)## With (hypothetical) calendar timessimDat2 <- simDatsimDat2$t.start <- as.Date(simDat2$t.start + simDat2$x2 * 5, origin = "20-01-01")simDat2$t.stop <- as.Date(simDat2$t.stop + simDat2$x2 * 5, origin = "20-01-01")reObj3 <- with(simDat2, Recur(t.start %to% t.stop, id, event, status))plot(reObj3, event.calendarTime = TRUE)## MCF plotsplot(reObj, mcf = TRUE)plot(reObj, mcf = TRUE, mcf.adjustRiskset = FALSE)library(reReg)data(simDat)reObj <- with(simDat, Recur(t.start %to% t.stop, id, event, status))summary(reObj)

Plot the Baseline Cumulative Rate Function and the Baseline Cumulative Hazard Function

Description

Plot the baseline cumulative rate function and the baseline cumulative hazard function(if applicable) for anreReg object.

Usage

## S3 method for class 'reReg'plot(  x,  baseline = c("both", "rate", "hazard"),  smooth = FALSE,  newdata = NULL,  frailty = NULL,  showName = FALSE,  control = list(),  ...)

Arguments

Details

The argumentcontrol consists of options with argument defaults to a listwith the following values:

xlab

customizable x-label, default value is "Time".

ylab

customizable y-label, default value is empty.

main

customizable title, default value are "Baseline cumulative rate andhazard function" whenbaseline = "both","Baseline cumulative rate function" whenbaseline = "rate",and "Baseline cumulative hazard function" whenbaseline = "hazard".

Value

Aggplot object.

See Also

reReg

Examples

data(simDat)fm <- Recur(t.start %to% t.stop, id, event, status) ~ x1 + x2fit <- reReg(fm, data = simDat, B = 0)plot(fit)plot(fit, xlab = "Time (days)", smooth = TRUE)## Predicted cumulative rate and hazard given covariatesnewdata <- expand.grid(x1 = 0:1, x2 = mean(simDat$x2))plot(fit, newdata = newdata, showName = TRUE)

Produce Event Plots

Description

Plot the event plot for anRecur object.The usage of the function is similar to that ofplot.Recur() but with more flexible options.

Usage

plotEvents(  formula,  data,  result = c("increasing", "decreasing", "asis"),  calendarTime = FALSE,  control = list(),  ...)

Arguments

Details

The argumentcontrol consists of options with argument defaults to a list withthe following values:

xlab

customizable x-label, default value is "Time".

ylab

customizable y-label, default value is "Subject" for event plot and"Cumulative mean" for MCF plot.

main

customizable title, the default value is "Recurrent event plot"whenmcf = FALSE and"Sample cumulative mean function plot" whenmcf = TRUE.

terminal.name

customizable label for terminal event,the default value is "Terminal event".

recurrent.name

customizable legend title for recurrent event,the default value is "Recurrent events".

recurrent.types

customizable label for recurrent event type,the default value isNULL.

alpha

between 0 and 1, controls the transparency of points.

Thexlab,ylab andmain parameters can be specifiedoutside of thecontrol list.

Value

Aggplot object.

See Also

Recur,plot.Recur

Examples

data(simDat)plotEvents(Recur(t.start %to% t.stop, id, event, status) ~ 1, data = simDat,           xlab = "Time in days", ylab = "Subjects arranged by terminal time")## Separate plots by x1plotEvents(Recur(t.start %to% t.stop, id, event, status) ~ x1, data = simDat)## For multiple recurrent eventssimDat$x3 <- ifelse(simDat$x2 < 0, "x2 < 0", "x2 > 0")simDat$event <- simDat$event * sample(1:3, nrow(simDat), TRUE)plotEvents(Recur(t.start %to% t.stop, id, event, status) ~ x1 + x3, data = simDat)

Plot options for plotEvents

Description

This function provides the plotting options for theplotEvents() function.

Usage

plotEvents.control(  xlab = NULL,  ylab = NULL,  main = NULL,  terminal.name = NULL,  recurrent.name = NULL,  recurrent.type = NULL,  legend.position = NULL,  base_size = 12,  cex = NULL,  alpha = 0.7,  width = NULL,  bar.color = NULL,  recurrent.color = NULL,  recurrent.shape = NULL,  recurrent.stroke = NULL,  terminal.color = NULL,  terminal.shape = NULL,  terminal.stroke = NULL,  not.terminal.color = NULL,  not.terminal.shape = NULL)

Arguments

See Also

plotEvents

Plot the Baseline Cumulative Hazard Function for the Terminal Time

Description

Plot the baseline cumulative hazard function for anreReg object.The 95% confidence interval on the baseline cumulative rate function

Usage

plotHaz(  x,  newdata = NULL,  frailty = NULL,  showName = FALSE,  type = c("unrestricted", "bounded", "scaled"),  smooth = FALSE,  control = list(),  ...)

Arguments

Details

The argumentcontrol consists of options with argumentdefaults to a list with the following values:

xlab

customizable x-label, default value is "Time".

ylab

customizable y-label, default value is empty.

main

customizable title, default value is "Baseline cumulative hazard function".

These arguments can also be passed down without specifying acontrol list.

Value

Aggplot object.

See Also

reRegplot.reReg

Examples

data(simDat)fm <- Recur(t.start %to% t.stop, id, event, status) ~ x1 + x2fit <- reReg(fm, data = simDat, model = "cox|cox", B = 0)## Plot both the baseline cumulative rate and hazard functionplot(fit)## Plot baseline cumulative hazard functionplotHaz(fit)plotHaz(fit, smooth = TRUE)

Plotting the Baseline Cumulative Rate Function for the Recurrent Event Process

Description

Plot the baseline cumulative rate function for anreReg object.

Usage

plotRate(  x,  newdata = NULL,  frailty = NULL,  showName = FALSE,  type = c("unrestricted", "bounded", "scaled"),  smooth = FALSE,  control = list(),  ...)

Arguments

Details

TheplotRate() plots the estimated baseline cumulative rate function depending on the identifiability assumption.Whentype = "unrestricted" (default), the baseline cumulative rate functionis plotted under the assumptionE(Z) = 1.Whentype = "scaled", the baseline cumulative rate function is plottedunder the assumption\Lambda(\min(Y^\ast, \tau)) = 1.Whentype = "bounded", the baseline cumulative rate function is plottedunder the assumption\Lambda(\tau) = 1.See?reReg for the specification of the notations and underlying models.

The argumentcontrol consists of options with argument defaultsto a list with the following values:

xlab

customizable x-label, default value is "Time".

ylab

customizable y-label, default value is empty.

main

customizable title, default value is "Baseline cumulative rate function".

These arguments can also be specified outside of thecontrol list.

Value

Aggplot object.

See Also

reRegplot.reReg

Examples

data(simDat)fm <- Recur(t.start %to% t.stop, id, event, status) ~ x1 + x2fit <- reReg(fm, data = simDat, model = "cox|cox", B = 0)## Plot both the baseline cumulative rate and hazard functionplot(fit)## Plot baseline cumulative rate functionplotRate(fit)plotRate(fit, smooth = TRUE)

Fits Semiparametric Regression Models for Recurrent Event Data

Description

Fits a general (joint) semiparametric regression model for the recurrent event data,where the rate function of the underlying recurrent event process andthe 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.

Usage

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

Arguments

Details

Model specification:

Suppose the recurrent event process and the failure events areobserved in the time intervalt\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 an byp 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 themodel argument with the rate functionand the hazard function separated by "|".For examples,Wang, Qin and Chiang (2001) (\alpha = \eta = \theta = 0)can be called withmodel = "cox";Huang and Wang (2004) (\alpha = \eta = 0)can be called withmodel = "cox|cox";Xu et al. (2017) (\alpha = \beta and\eta = \theta)can be called withmodel = "am|am";Xu et al. (2019) (\eta = \theta = 0) can be called withmodel = "gsc".Users can mix the models depending on the application. For example,model = "cox|ar" postulate a Cox proportional model for therecurrent event rate function and an accelerated rate model forthe 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, specifyingmodel = "cox" is equivalent tomodel = "cox|cox".Some models that assumesZ = 1 and requires independentcensoring are also implemented inreReg;these includesmodel = "cox.LWYY" for Lin et al. (2000),model = "cox.GL" for Ghosh and Lin (2002),andmodel = "am.GL" for Ghosh and Lin (2003).Additionally, an improved estimation of the proportional rate model(Huang and Huang 2022) can be called bymodel = "cox.HH" withadditionalcontrol options to specify the underlying procedure.Seeonline vignettefor a detailed discussion of the implemented regression models.

Variance estimation:

The available methods for variance estimation are:

boot

performs nonparametric bootstrap.

sand

performs the efficient resampling-based variance estimation.

Improving proportional rate model:A common semiparametric regression model for recurrent event processunder the noninformative censoring assumption is the Cox-type proportional rate model(available inreReg() viamodel = "cox.LWYY").However, the construction of the pseudo-partial score function ignores thedependency among recurrent events and thus could be inefficient. To improve upon this popular method, Huang and Huang (2022) proposed to combinea system of weighted pseudo-partial score equations via the generalized method of moments (GMM)and empirical likelihood (EL) estimation.The proposed GMM and EL procedures are available inreReg viamodel = "cox.HH"with additional control specifications.Seeonline vignettefor an illustration of this feature.

Control options:

Thecontrol list consists of the following parameters:

tol

absolute error tolerance.

init

a list contains initial guesses used for root search.

solver

the equation solver used for root search.The available options areBB::BBsolve,BB::dfsane,BB::BBoptim,optimx::optimr,dfoptim::hjk,dfoptim::mads,optim,andnleqslv::nleqslv.

eqType

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

boot.parallel

an logical value indicating whether parallel computationwill be applied whense = "boot" is called.

boot.parCl

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

cppl

A character string indicating either to improve the proportional rate model viathe generalized method of moments (cppl = "GMM") or empirical likelihood estimation (cppl = "EL").This option is only used whenmodel = "cox.HH".

cppl.wfun

A list of (up to two) weight functions to be combined with the weighted pseudo-partial likelihood scores.Available options are"Gehan" and"cumbase",which correspond to the Gehan's weight and the cumulative baseline hazard function, respectively.Alternatively, the weight functions can be specified with function formulas.This option is only used whenmodel = "cox.HH".

trace

A logical variable denoting whether some of theintermediate results of iterations should be displayed to the user. Default isFALSE.

References

Chiou, S.H., Xu, G., Yan, J. and Huang, C.-Y. (2023). Regression Modeling for Recurrent Events Possibly with an Informative Terminal Event Using R Package reReg.Journal of Statistical Software,105(5): 1–34.

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 EventProcesses under Informative Censoring.Statistica Sinica,30: 1773–1795.

Huang, M.-Y. and Huang, C.-Y. (2022). Improved semiparametric estimation of the proportional rate model with recurrent event data.Biometrics,79 3: 1686–1700.

See Also

Recur,simGSC

Examples

data(simDat)## Nonparametric estimateplot(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 modelsummary(reReg(fm, data = simDat, model = "cox", B = 50))## Fit the joint Cox/Cox modelsummary(reReg(fm, data = simDat, model = "cox|cox", B = 50))## Fit the scale-change rate modelsummary(reReg(fm, data = simDat, model = "gsc", B = 50, se = "sand"))

Package options for reReg

Description

This function provides the fitting options for thereReg() function.

Usage

reReg.control(  eqType = c("logrank", "gehan", "gehan_s"),  solver = c("BB::dfsane", "BB::BBsolve", "BB::BBoptim", "optimx::optimr",    "dfoptim::hjk", "dfoptim::mads", "optim", "nleqslv::nleqslv"),  tol = 1e-07,  cppl = NULL,  cppl.wfun = list(NULL, NULL),  init = list(alpha = 0, beta = 0, eta = 0, theta = 0),  boot.parallel = FALSE,  boot.parCl = NULL,  maxit1 = 100,  maxit2 = 10,  trace = FALSE,  numAdj = 1e-07)

Arguments

See Also

reReg

Create anreSurv Object

Description

Create a recurrent event survival object, used as a response variable inreReg.This function is deprecated in Version 1.1.6.A recurrent event object is now being created withRecur().See '?Recur()' for details.

Usage

reSurv(time1, time2, id, event, status, origin = 0)

Arguments

Examples

## Not run:   data(simDat)  ## being deprecated in Verson 1.1.7  with(dat, reSurv(Time, id, event, status))  ## Use Recur() instead  with(dat, Recur(Time, id, event, status))## End(Not run)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the linksbelow to see their documentation.

reda

%2%,%to%,check_Recur,is.Recur,mcf,Recur

Calculate Residuals for a ‘reReg’ Fit

Description

Calculates residuals for a joint frailty scale-change model fitted by 'reReg'.Under the recurrent event model, at each observation time,t,the residual is calculated as

\mbox{observed number of recurrent events at } t- \mbox{expected number of recurrent events at} t.

The expected number of recurrent events att is calculated by thecumulative rate function att.Under the failure time model, the residual is calculated as

\Delta - H(t),

where\Delta is the terminal event indicator andH(t) is the cumulative hazard function att.

Usage

## S3 method for class 'reReg'residuals(object, model = c("recurrent", "failure"), ...)

Arguments

Simulated dataset for demonstration

Description

A simulated data frame with the following variables

id

subjects identification

t.start

start of the interval

t.stop

endpoint of the interval; when time origin is 0 this variable also marks the recurrence or terminal/censoring time

status

terminal event indicator; 1 if a terminal event is recorded

event

recurrent event indicator; 1 if a recurrent event is recorded

x1

baseline covariate generated from a standard uniform distribution

x2

baseline covariate generated from a standard uniform distribution (independent fromz1

Usage

data(simDat)

Format

A data frame with 874 rows and 7 variables.

Details

SeesimGSC for instruction on simulating recurrent event data fromscale-change models.

Function to generate simulated recurrent event data

Description

The functionsimGSC() generates simulated recurrent event data from eithera Cox-type model, an accelerated mean model, an accelerated rate model, or a generalized scale-change model.

Usage

simGSC(  n,  summary = FALSE,  para,  xmat,  censoring,  frailty,  tau,  origin,  Lam0,  Haz0)

Arguments

Details

The functionsimGSC() generates simulated recurrent event data overthe interval(0, \tau) based on the specification of the recurrent process andthe terminal events.Specifically, the rate function,\lambda(t), of the recurrent processcan be specified as one of the following model:

\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 an byp covariate matrix and\alpha,Z is an unobserved shared frailty variable, and(\alpha, \eta) and(\beta, \theta) correspond to the shape and size parameters of therate function and the hazard function, respectively.

Under the default settings, thesimGSC() function assumesp = 2and the regression parameters to be\alpha = \eta = (0, 0)^\top,and\beta = \theta = (1, 1)^\top.When thexmat argument is not specified, thesimGSC() functionassumesX_i is a two-dimensional vectorX_i = (X_{i1}, X_{i2}), i = 1, \ldots, n,whereX_{i1} is a Bernoulli variable with rate 0.5 andX_{i2} is a standard normal variable.With the defaultxmat, the censoring time $C$ is generated froman exponential distribution with mean\tau X_{i1} + Z^2\tau(1 - X_{i1}).Thus, the censoring distribution is covariate dependent andis informative whenZ is not a constant.When thefrailty argument is not specified, the frailty variableZ is generatedfrom a gamma distribution with a unit mean and a variance of 0.25.The default values fortau andorigin are 60 and 0, respectively.When argumentsLam0 andHaz0 are left unspecified,thesimGSC() function uses\Lambda_0(t) = 2\log(1 + t)andH_0(t) = \log(1 + t) / 5, respectively.This is equivalent to settingLam0 = function(x) 2 * log(1 + x) andHaz0 = function(x) log(1 + x) / 5.Overall, the default specifications generate the recurrent events and the terminal eventsfrom the model:

\lambda(t) = \displaystyle \frac{2Z}{1 + te^{-X_{i1} - X_{i2}}},h(t) = \displaystyle \frac{Z}{5(1 + te^{X_{i1} + X_{i2}})}, t\in[0, 60].

Seeonline vignettefor more examples.

See Also

reReg

Examples

set.seed(123)simGSC(100, summary = TRUE)