Movatterモバイル変換

Title:	Recurrent Event Data Analysis
Version:	0.5.6
Description:	Contains implementations of recurrent event data analysis routines including (1) survival and recurrent event data simulation from stochastic process point of view by the thinning method proposed by Lewis and Shedler (1979) <doi:10.1002/nav.3800260304> and the inversion method introduced in Cinlar (1975, ISBN:978-0486497976), (2) the mean cumulative function (MCF) estimation by the Nelson-Aalen estimator of the cumulative hazard rate function, (3) two-sample recurrent event responses comparison with the pseudo-score tests proposed by Lawless and Nadeau (1995) <doi:10.2307/1269617>, (4) gamma frailty model with spline rate function following Fu, et al. (2016) <doi:10.1080/10543406.2014.992524>.
Imports:	Rcpp, graphics, methods, splines2 (≥ 0.4.3), stats
LinkingTo:	Rcpp, RcppArmadillo
Suggests:	ggplot2, grDevices, knitr, rmarkdown, tinytest,
Depends:	R (≥ 3.2.3)
License:	GPL (≥ 3)
LazyData:	true
VignetteBuilder:	knitr
Collate:	'RcppExports.R' 'class.R' 'Recur.R' 'Survr.R' 'aic.R''baseline.R' 'coef.R' 'data.R' 'mcf-generic.R' 'mcf-formula.R''mcf-rateReg.R' 'mcfDiff.R' 'misc.R' 'plot.R' 'rateReg.R''reda-package.R' 'show.R' 'simEvent.R' 'summary.R' 'zzz.R'
URL:	https://wwenjie.org/reda,https://github.com/wenjie2wang/reda
BugReports:	https://github.com/wenjie2wang/reda/issues
Encoding:	UTF-8
RoxygenNote:	7.3.2
NeedsCompilation:	yes
Packaged:	2025-09-02 13:02:51 UTC; wenjie
Author:	Wenjie Wang [aut, cre], Haoda Fu [aut], Sy Han (Steven) Chiou [ctb], Jun Yan [ctb]
Maintainer:	Wenjie Wang <wang@wwenjie.org>
Repository:	CRAN
Date/Publication:	2025-09-02 14:00:07 UTC

Recurrent Event Data Analysis

Description

The R packagereda provides functions for simulating, exploring andmodeling recurrent event data.

Details

The main functions are summarized as follows:

• simEventData: Simulating survival, recurrent event, andmultiple event data from stochastic process point of view.

• mcf: Estimating the mean cumulative function (MCF) from afitted gamma frailty model, or from a sample recurrent event data by usingthe nonparametic MCF estimator (the Nelson-Aelen estimator of the cumulativehazard function).

• mcfDiff: Comparing two-sample MCFs by the pseudo-score testsand estimating their difference over time.

• rateReg: Fitting Gamma fraitly model with spline baseline ratefunction.

See the package vignettes for more introduction and demonstration.

Author(s)

Maintainer: Wenjie Wangwang@wwenjie.org (ORCID)

Authors:

• Haoda Fu

Other contributors:

• Sy Han (Steven) Chiou [contributor]

• Jun Yan (ORCID) [contributor]

See Also

Useful links:

• https://wwenjie.org/reda

• https://github.com/wenjie2wang/reda

• Report bugs athttps://github.com/wenjie2wang/reda/issues

Akaike Information Criterion (AIC)

Description

AIC,rateReg-method is an S4 class method calculating Akaikeinformation criterion (AIC) for one or severalrateReg objects,according to the formula - 2 * log-likelihood + 2 * nPar, where nParrepresents the number of parameters in the fitted model.

Usage

## S4 method for signature 'rateReg'AIC(object, ..., k = 2)

Arguments

Details

When comparing models fitted by maximum likelihood to the same data, thesmaller the AIC, the better the fit. A friendly warning will be thrown outif the numbers of observation were different in the model comparison.help(AIC, stats) for other details.

Value

If just one object is provided, a numeric value representingcalculated AIC. If multiple objects are provided, a data frame withrows corresponding to the objects and columnsdf andAIC,wheredf means degree of freedom, which is the number ofparameters in the fitted model.

See Also

rateReg for model fitting;summary,rateReg-method for summary of a fitted model;BIC,rateReg-method for BIC.

Examples

## See examples given in function rateReg.

Bayesian Information Criterion (BIC)

Description

BIC,rateReg-method is an S4 class method calculatingBayesian information criterion (BIC) or so-calledSchwarz's Bayesian criterion (SBC)for one or severalrateReg objects,according to the formula- 2 * log-likelihood + ln(nObs) * nPar,where nPar represents the number of parameters in the fitted modeland nObs is the number of observations.

Usage

## S4 method for signature 'rateReg'BIC(object, ...)

Arguments

Details

When comparing models fitted by maximum likelihood to the samedata, the smaller the BIC, the better the fit.help(BIC, stats) for other details.

Value

If just one object is provided, a numeric value representingcalculated BIC.If multiple objects are provided, a data frame with rowscorresponding to the objects and columnsdf andBIC,wheredf means degree of freedom,which is the number of parameters in the fitted model.

See Also

rateReg for model fitting;summary,rateReg-method for summary of a fitted model;AIC,rateReg-method for AIC.

Examples

## See examples given in function rateReg.

Formula Response for Recurrent Event Data

Description

Create an S4 class object that represents formula response for recurrentevent data with optional checking procedures embedded.

Usage

Recur(  time,  id,  event,  terminal,  origin,  check = c("hard", "soft", "none"),  ...)

Arguments

Details

This is a successor function of the deprecated functionSurvr. Seethe vignette by 'vignette("reda-Recur")' for details.

Value

AnRecur object.

Examples

library(reda)with(valveSeats, Recur(Days, ID))with(valveSeats, Recur(Days, ID, No.))with(valveSeats, Recur(Days, ID, No., terminal = 1))with(valveSeats, Recur(Days, ID, No., origin = 10))

An S4 Class Representing Formula Response for Recurrent Event Data

Description

The classRecur is an S4 that represents a formula response forrecurrent event data model. The functionRecur producesobjects of this class. See “Slots” for details.

Slots

.Data

A numeric matrix that consists of the following columns:

• time1: the beginning of time segements;

• time2: the end of time segements;

• id: Identificatorsof subjects;

• event: Event indicators;

• :terminal: Indicators of terminal events.

call

A function call producing the object.

ID

A character vector for unique original identificators of subjects.

ord

An integer vector for increasingly ordering data byid,time2, and- event. Sorting is often done in themodel-fitting steps, where the indices stored in this slot can be useddirectly.

rev_ord

An integer vector for reverting the ordering of the sorteddata (byord) to its original ordering. This slot is provided toeasily revert the sorting.

first_idx

An integer vector indicating the first record of eachsubject in the sorted matrix. It helps in the data checking produce andmay be helpful in model-fitting step, such as getting the origin time.

last_idx

An integer vector indicating the last record of each subjectin the sorted data. Similar tofirst_idx, it helps in the datachecking produce and may be helpful in the model-fitting step, such aslocating the terminal events.

check

A character string indicating how the data checking isperformed. It just records the option that users specified on datachecking.

time_class

A character vector preserving the class(es) of input times.

See Also

Recur

Recurrent Episodes

Description

Specify time segements or recurrent episodes by endpoints.

Usage

time1 %to% time2time1 %2% time2

Arguments

Details

This function is intended to be used for specifying the argumenttimein functionRecur.

Value

A list that consists of two elements named"time1" and"time2".

Formula Response for Recurrent Event Data

Description

Create an S4 class that represents formula response for recurrent event datamodeled by methods based on counts and rate function. Note that thefunction is deprecated since version 0.5.0 and will be removed in future.

Usage

Survr(ID, time, event, origin = 0, check = TRUE, ...)

Arguments

Details

This is a similar function toSurvr in packagesurvrec but with a more considerate checking procedure embedded forrecurrent event data modeled by methods based on counts and rate function.The checking rules apply to each subject respectively and include that

• Subject identification, event times, censoring time, and eventindicator cannot be missing or contain missing values.

• There has to be only one censoring time not earlier thanany event time.

• The time origin has to be the same and not later than any eventtime.

An S4 Class Representing Formula Response

Description

The classSurvr is an S4 that represents a formula response forrecurrent event data model. The functionSurvr producesobjects of this class. See “Slots” for details.

Slots

.Data

A numeric matrix object.

ID

A charactrer vector for original subject identificator.

check

A logical value indicating whether to performance data checking.

ord

An integer vector for increasingly ordering data byID,time, and1 - event.

See Also

Survr

Convert An Recur Object to A Character Vector

Description

Summarize and convert the recurrent episodes for each subjects intocharacter strings.

Usage

## S4 method for signature 'Recur'as.character(x, ...)

Arguments

Details

This function is intended to be a helper function for the 'show()' method of'Recur' objects. To be precise, the function set the maximum number ofrecurrent episodes for each subject to be 'max(2L,as.integer(getOption("reda.Recur.maxPrint")))'. By default, at most threerecurrent episodes will be summarized for printing. When subjects havingmore than three recurrent episodes, the first'getOption("reda.Recur.maxPrint") - 1' number of recurrent episodes and thelast one will be summarized. One may use 'options()' to adjust the setting.For example, the default value is equivalent to 'options(reda.Recur.maxPrint= 3)'.

Estimated Baseline Rate Function

Description

An S4 class generic function that returns the estimated baseline ratefunction.

Usage

baseRate(object, ...)## S4 method for signature 'rateReg'baseRate(object, level = 0.95, control = list(), ...)

Arguments

Value

AbaseRate object.

Functions

• baseRate(rateReg): Estiamted baseline rate from a fitted model.

See Also

rateReg for model fitting;summary,rateReg-method for summary of a fitted model;plot,baseRate.rateReg-method for ploting method.

Examples

## See examples given in function rateReg.

An S4 Class Representing Estimated Baseline Rate Function

Description

An S4 class that represents the estimated baseline rate function from model.The functionbaseRate produces objects of this class.

Slots

baseRate

A data frame.

level

A numeric value.

See Also

baseRate,rateReg-method

Checks for Recurrent Event Data

Description

Perform several checks for recurrent event data and update objectattributions if some rows of the contained data (in the.Data slot)have been removed by such asna.action.

Usage

check_Recur(x, check = c("hard", "soft", "none"))

Arguments

Value

AnRecur object invisibly.

Estimated Coefficients of Covariates

Description

coef,rateReg-method is an S4 class method that extracts estimatedcoefficients of covariates fromrateReg object produced by functionrateReg.

Usage

## S4 method for signature 'rateReg'coef(object, ...)

Arguments

Value

A named numeric vector.

See Also

rateReg for model fitting;confint,rateReg-method for confidence intervalsfor covariate coefficients;summary,rateReg-method for summary of a fitted model.

Examples

## See examples given in function rateReg.

Confidence Intervals for Covariate Coefficients

Description

confint,rateReg-method is an S4 class method forrateReg object, which returns approximate confidence intervalsfor all or specified covariates.

Usage

## S4 method for signature 'rateReg'confint(object, parm, level = 0.95, ...)

Arguments

Details

Under regularity condition (Shao 2003, Theorem 4.16 and Theorem 4.17, page287, 290), the approximate confidence intervals are constructed looselybased on Fisher information matrix and estimates of coefficients.

Value

A numeric matrix with row names and column names.

References

Shao, J. (2003),Mathematical statistics,Springer texts in statistics, New York: Springer, 2nd Edition.

See Also

rateReg for model fitting;coef,rateReg-method for point estimatesof covariate coefficients;summary,rateReg-method for summary of a fitted model.

Examples

## See examples given in function rateReg.

Is the xect from the Recur class?

Description

ReturnTRUE if the specified xect is from theRecurclass,FALSE otherwise.

Usage

is.Recur(x)

Arguments

Value

A logical value.

Mean Cumulative Function (MCF)

Description

An S4 class generic function that returns the mean cumulative function (MCF)estimates from a fitted model or returns the nonparametric MCF estimates (byNelson-Aalen estimator or Cook-Lawless cumulative sample mean estimator)from the sample data.

Usage

mcf(object, ...)## S4 method for signature 'formula'mcf(  object,  data,  subset,  na.action,  variance = c("LawlessNadeau", "Poisson", "bootstrap", "CSV", "none"),  logConfInt = FALSE,  adjustRiskset = TRUE,  level = 0.95,  control = list(),  ...)## S4 method for signature 'rateReg'mcf(  object,  newdata,  groupName,  groupLevels,  level = 0.95,  na.action,  control = list(),  ...)

Arguments

Details

Forformula object withRecur object as response, thecovariate specified at the right hand side of the formula should be either1 or any "linear" conbination of categorical variable in the data.The former computes the overall sample MCF. The latter computes the sampleMCF for each level of the combination of the categorical variable(s)specified, respectively.

The MCF estimates are computed on each unique time point of the sample data.By default, the size of risk set is adjusted over time based on the at-riskindicators, which results in the Nelson-Aalen nonparametric estimator(Nelson 2003). If the size of risk set remains a constant (total number ofprocesses) over time (specified byadjustRiskset = FALSE), thecumulative sample mean (CSM) function introduced in Chapter 1 of Cook andLawless (2007) will be computed instead. The point estimate of sample MCFat each time point does not assume any particular underlying model. Thevariance estimates at each time point is computed following the Lawless andNadeau method (LawLess and Nadeau 1995), the Poisson process method, or thebootstrap methods. The approximate confidence intervals are provided aswell, which are constructed based on the asymptotic normality of the MCFitself (by default) or the logarithm of MCF.

ForrateReg object,mcf estimates the baseline MCF and itsconfidence interval at each time grid if argumentnewdata is notspecified. Otherwise,mcf estimates MCF and its confidence intervalfor the givennewdata based on Delta-method.

Value

Amcf.formula ormcf.rateReg object.

A brief description of the slots of amcf.formula object is given asfollows:

• formula: Model Formula.

• data: Processed data based on the model formula or anempty data frame ifkeep.data is set to beFALSE.

• MCF: A data frame containing estimates for sample MCF.

• origin: Time origins.

• multiGroup: A logical value indicating whether MCFis estimated for different groups respectively.

• logConfInt: A logical value indicating whether thevariance estimates are based on the normality of logarithm ofthe MCF estimates.

• level: Confidence level specified.

Most slots of amcf.rateReg object are inherited from the inputrateReg object. A brief description of other slots is given asfollows:

• newdata: Given dataset used to estimate MCF.

• MCF: A data frame containing MCF estimates.

• level: Confidence level specified.

• na.action: The way handling missing values.

• control: The control list.

• multiGroup: A logical value indicating whether MCFis estimated for different groups respectively.

Functions

• mcf(formula): Sample MCF from data.

• mcf(rateReg): Estimated MCF from a fitted model.

References

Cook, R. J., and Lawless, J. (2007).The statistical analysis ofrecurrent events, Springer Science & Business Media.

Lawless, J. F. and Nadeau, C. (1995). Some Simple Robust Methods for theAnalysis of Recurrent Events.Technometrics, 37, 158–168.

Nelson, W. B. (2003).Recurrent Events Data Analysis for ProductRepairs, Disease Recurrences, and Other Applications (Vol. 10). SIAM.

See Also

rateReg for model fitting;mcfDiff for comparing two-sample MCFs.plot-method for plotting MCF.

Examples

library(reda)### sample MCF## Example 1. valve-seat data## the default variance estimates by Lawless and Nadeau (1995) methodvalveMcf0 <- mcf(Recur(Days, ID, No.) ~ 1, data = valveSeats)plot(valveMcf0, conf.int = TRUE, mark.time = TRUE, addOrigin = TRUE) +    ggplot2::xlab("Days") + ggplot2::theme_bw()## variance estimates following Poisson process modelvalveMcf1 <- mcf(Recur(Days, ID, No.) ~ 1,                 data = valveSeats, variance = "Poisson")## variance estimates by bootstrap method (with 1,000 bootstrap samples)set.seed(123)valveMcf2 <- mcf(Recur(Days, ID, No.) ~ 1,                 data = valveSeats, variance = "bootstrap",                 control = list(B = 200))## comparing the variance estimates from different methodslibrary(ggplot2)ciDat <- rbind(cbind(valveMcf0@MCF, Method = "Lawless & Nadeau"),               cbind(valveMcf1@MCF, Method = "Poisson"),               cbind(valveMcf2@MCF, Method = "Bootstrap"))ggplot(ciDat, aes(x = time, y = se)) +    geom_step(aes(color = Method, linetype = Method)) +    xlab("Days") + ylab("SE estimates") + theme_bw()## comparing the confidence interval estimates from different methodsggplot(ciDat, aes(x = time)) +    geom_step(aes(y = MCF)) +    geom_step(aes(y = lower, color = Method, linetype = Method)) +    geom_step(aes(y = upper, color = Method, linetype = Method)) +    xlab("Days") + ylab("Confidence intervals") + theme_bw()## Example 2. the simulated datasimuMcf <- mcf(Recur(time, ID, event) ~ group + gender,               data = simuDat, ID %in% 1 : 50)plot(simuMcf, conf.int = TRUE, lty = 1 : 4,     legendName = "Treatment & Gender")### estimate MCF difference between two groups## one sample MCF object of two groupsmcf0 <- mcf(Recur(time, ID, event) ~ group, data = simuDat)## two-sample pseudo-score testsmcfDiff.test(mcf0)## difference estimates over timemcf0_diff <- mcfDiff(mcf0, testVariance = "none")plot(mcf0_diff)## or explicitly ask for the difference of two sample MCFmcf1 <- mcf(Recur(time, ID, event) ~ 1, data = simuDat,            subset = group %in% "Contr")mcf2 <- mcf(Recur(time, ID, event) ~ 1, data = simuDat,            subset = group %in% "Treat")## perform two-sample tests and estimate difference at the same timemcf12_diff1 <- mcfDiff(mcf1, mcf2)mcf12_diff2 <- mcf1 - mcf2   # or equivalently using the `-` methodstopifnot(all.equal(mcf12_diff1, mcf12_diff2))mcf12_diff1plot(mcf12_diff1)### For estimated MCF from a fitted model,### see examples given in function rateReg.

An S4 Class Representing Sample MCF

Description

An S4 class that represents sample mean cumulative function (MCF) from data.The functionmcf produces objects of this class.

Slots

formula

Formula.

data

A data frame.

MCF

A data frame.

origin

A named numeric vector.

multiGroup

A logical value.

variance

A character vector.

logConfInt

A logical value.

level

A numeric value.

See Also

mcf,formula-method.

An S4 Class Respresenting Estimated MCF from a Fitted Model

Description

An S4 class that represents estimated mean cumulative function (MCF) fromModels. The functionmcf produces objects of this class.

Slots

call

Function call.

formula

Formula.

spline

A character.

knots

A numeric vector.

degree

A nonnegative integer.

Boundary.knots

A numeric vector.

newdata

A numeric matrix.

MCF

A data frame.

level

A numeric value between 0 and 1.

na.action

A length-one character vector.

control

A list.

multiGroup

A logical value.

See Also

mcf,rateReg-method

Comparing Two-Sample MCFs

Description

This function estimates the sample MCF difference between two groups. Boththe point estimates and the confidence intervals are computed (Lawless andNadeau 1995). The two-sample pseudo-score test proposed by Cook, Lawless,and Nadeau (1996) is also performed by default.

Usage

mcfDiff(mcf1, mcf2 = NULL, level = 0.95, ...)mcfDiff.test(  mcf1,  mcf2 = NULL,  testVariance = c("robust", "Poisson", "none"),  ...)

Arguments

Details

The functionmcfDiff estimates the two-sample MCFs' difference andinternally calls functionmcfDiff.test to perform the pseudo-scoretests by default. A- method is available as a simple wrapper for thefunctionmcfDiff for comparing two-sample MCFs from twomcf.formula objects. For instance, supposemcf1 andmcf2 aremcf.formula objects, each of which represents thesample MCF estimates for one group. The function callmcf1 - mcf2 isequivalent tomcfDiff(mcf1, mcf2).

The null hypothesis of the two-sample pseudo-score test is that there is nodifference between the two sample MCFs, while the alternative hypothesissuggests a difference. The test is based on a family of test statisticsproposed by Lawless and Nadeau (1995). The argumenttestVariancespecifies the method for computing the variance estimates of the teststatistics under different model assumption. See the document of argumenttestVariance for all applicable options. For the variance estimatesrobust to departures from Poisson process assumption, both constant weightand the linear weight function (with scaling) suggested in Cook, Lawless,and Nadeau (1996) are implemented. The constant weight is powerful in caseswhere the two MCFs are approximately proportional to each other. The linearweight function is originallya(u) = t - u, whereu representsthe time variable andt is the first time point when the risk set ofeither group becomes empty. It is further scaled by1 / t for teststatistics invariant to the unit of measurement of the time variable. Thelinear weight function puts more emphasis on the difference at earily timesthan later times and is more powerful for cases where the MCFs are no longerproportional to each other, but not crossing. Also see Cook and Lawless(2007, Section 3.7.5) for more details.

Value

The functionmcfDiff returns amcfDiff object (of S4 class)that contains the following slots:

• call: Function call.

• MCF: Estimated Mean cumulative function Difference at each timepoint.

• origin: Time origins of the two groups.

• variance: The method used for variance estimates.

• logConfInt: A logical value indicating whether normality isassumed forlog(MCF) instead of MCF itself. FormcfDiffobject, it is alwaysFALSE.

• level: Confidence level specified.

• test: AmcfDiff.test object for the hypothesis testresults.

The functionmcfDiff.test returns amcfDiff.test object (of S4class) that contains the following slots:

• .Data: A numeric matrix (of two rows and five columns) forhypothesis testing results.

• testVariance: A character string (or vector of length one)indicating the method used for the variance estimates of the test statistic.

References

Lawless, J. F., & Nadeau, C. (1995). Some Simple Robust Methods for theAnalysis of Recurrent Events.Technometrics, 37(2), 158–168.

Cook, R. J., Lawless, J. F., & Nadeau, C. (1996). Robust Tests for TreatmentComparisons Based on Recurrent Event Responses.Biometrics, 52(2),557–571.

Cook, R. J., & Lawless, J. (2007).The Statistical Analysis ofRecurrent Events. Springer Science & Business Media.

Examples

## See examples given for function mcf.

An S4 Class Representing Sample MCF Difference

Description

An S4 class that represents the difference between two sample meancumulative functions from data. The functionmcfDiffproduces objects of this class.

Slots

call

A function call.

MCF

A data frame.

origin

A named numeric vector.

variance

A character vector.

logConfInt

A logical value.

level

A numeric value.

test

AmcfDiff.test class object.

See Also

mcfDiff

An S4 Class Representing the Two-Sample Pseudo-Score Test Results

Description

An S4 class that represents the results of the two-sample pseudo-score testsbetween two sample mean cumulative functions. The functionmcfDiff.test produces objects of this class.

Slots

.Data

A numeric matrix.

testVariance

A character vector.

See Also

mcfDiff.test

Parametrizations of Covariates and Covariate Coefficients

Description

This function helps the parametrizations of covariates and covariatecoeffcients when users specify a general hazard rate function in functionsimEvent andsimEventData. It applies the specified function(or the built-in option)FUN to thei_{th} row of the covariatematrixz and thei_{th} row of the coefficient matrix,iteratively, fori from one to the number of rows of the covariatematrixz.

Usage

parametrize(z, zCoef, FUN = c("exponential", "linear", "excess"), ...)

Arguments

Value

A numeric vector.

See Also

simEvent

Examples

## time pointstimeVec <- c(0.5, 2)## time-variant covariateszMat <- cbind(0.5, ifelse(timeVec > 1, 1, 0))## time-varying coefficientszCoefMat <- cbind(sin(timeVec), timeVec)## the following three ways are equivalent for the exponential form,## where the first one (using the built-in option) has the best performanceparametrize(zMat, zCoefMat, FUN = "exponential")parametrize(zMat, zCoefMat, function(z, zCoef) exp(z %*% zCoef))sapply(1 : 2, function(i) as.numeric(exp(zMat[i, ] %*% zCoefMat[i, ])))

Plot Baseline Rate or Mean Cumulative Function (MCF)

Description

S4 class methods plotting sample MCF from data, estimated MCF, or estimatedbaseline hazard rate function from a fitted model by usingggplot2plotting system. The plots generated are thus able to be further customizedproperly.

Usage

## S4 method for signature 'mcf.formula,missing'plot(  x,  y,  lty,  col,  legendName,  legendLevels,  conf.int = FALSE,  mark.time = FALSE,  addOrigin = FALSE,  ...)## S4 method for signature 'mcf.rateReg,missing'plot(x, y, conf.int = FALSE, lty, col, ...)## S4 method for signature 'baseRate.rateReg,missing'plot(x, y, conf.int = FALSE, lty, col, ...)## S4 method for signature 'mcfDiff,missing'plot(  x,  y,  lty,  col,  legendName,  legendLevels,  conf.int = TRUE,  addOrigin = FALSE,  ...)

Arguments

Value

Aggplot object.

See Also

mcf for estimation of MCF;rateReg for model fitting.

Examples

## See examples given in function mcf and rateReg.

Recurrent Events Regression Based on Counts and Rate Function

Description

This function fits recurrent event data (event counts) by gamma frailtymodel with spline rate function. The default model is the gamma frailtymodel with one piece constant baseline rate function, which is equivalent tonegative binomial regression with the same shape and rate parameter in thegamma prior. Spline (including piecewise constant) baseline hazard ratefunction can be specified for the model fitting.

Usage

rateReg(  formula,  data,  subset,  na.action,  start = list(),  control = list(),  contrasts = NULL,  ...)rateReg.control(  df = NULL,  degree = 0L,  knots = NULL,  Boundary.knots = NULL,  periodic = FALSE,  verbose = TRUE,  ...)

Arguments

Details

FunctionRecur in the formula response by default first checksthe dataset and will report an error if the dataset does not fall intorecurrent event data framework. Subject's ID will be pinpointed if itsobservation violates any checking rule. SeeRecur for all thechecking rules.

FunctionrateReg first constructs the design matrix fromthe specified arguments:formula,data,subset,na.action andconstrasts before model fitting.The constructed design matrix will be checked again tofit the recurrent event data frameworkif any observation with missing covariates is removed.

The model fitting process involves minimization of negative loglikelihood function, which calls functionconstrOptiminternally.help(constrOptim) for more details.

The argumentstart is an optional listthat allows users to specify the initial guess forthe parameter values for the minimization ofnegative log likelihood function.The available numeric vector elements in the list include

• beta: Coefficient(s) of covariates,set to be all 0.1 by default.

• theta: Parameter in Gamma(theta, 1 / theta) forfrailty random effect, set to be 0.5 by default.

• alpha: Coefficient(s) of baseline rate function,set to be all 0.05 by default.

The argumentcontrol allows users to control the process ofminimization of negative log likelihood function passed toconstrOptim and specify the boundary knots of baseline rate function.

Value

ArateReg object, whose slots include

• call: Function call ofrateReg.

• formula: Formula used in the model fitting.

• nObs: Number of observations.

• spline: A list contains

  • spline: The name of splines used.

  • knots: Internal knots specified for the baselinerate function.

  • Boundary.knots: Boundary knots specified for thebaseline rate function.

  • degree: Degree of spline bases specified inbaseline rate function.

  • df: Degree of freedom of the model specified.

• estimates: Estimated coefficients of covariates andbaseline rate function, and estimated rate parameter ofgamma frailty variable.

• control: The control list specified for model fitting.

• start: The initial guess specified for the parametersto be estimated.

• na.action: The procedure specified to deal withmissing values in the covariate.

• xlevels: A list that records the levels ineach factor variable.

• contrasts: Contrasts specified and used for eachfactor variable.

• convergCode:code returned by functionoptim, which is an integer indicating why theoptimization process terminated.help(optim) for details.

• logL: Log likelihood of the fitted model.

• fisher: Observed Fisher information matrix.

References

Fu, H., Luo, J., & Qu, Y. (2016).Hypoglycemic events analysis via recurrent time-to-event (HEART) models.Journal Of Biopharmaceutical Statistics, 26(2), 280–298.

See Also

summary,rateReg-method for summary of fitted model;coef,rateReg-method for estimated covariate coefficients;confint,rateReg-method for confidence interval ofcovariate coefficients;baseRate,rateReg-method for estimated coefficients of baselinerate function;mcf,rateReg-method for estimated MCF from a fitted model;plot,mcf.rateReg-method for plotting estimated MCF.

Examples

library(reda)## constant rate function(constFit <- rateReg(Recur(time, ID, event) ~ group + x1, data = simuDat))## six pieces' piecewise constant rate function(piecesFit <- rateReg(Recur(time, ID, event) ~ group + x1,                      data = simuDat, subset = ID %in% 1:50,                      knots = seq.int(28, 140, by = 28)))## fit rate function with cubic spline(splineFit <- rateReg(Recur(time, ID, event) ~ group + x1, data = simuDat,                      knots = c(56, 84, 112), degree = 3))## more specific summarysummary(constFit)summary(piecesFit)summary(splineFit)## model selection based on AIC or BICAIC(constFit, piecesFit, splineFit)BIC(constFit, piecesFit, splineFit)## estimated covariate coefficientscoef(piecesFit)coef(splineFit)## confidence intervals for covariate coefficientsconfint(piecesFit)confint(splineFit, "x1", 0.9)confint(splineFit, 1, 0.975)## estimated baseline rate functionsplinesBase <- baseRate(splineFit)plot(splinesBase, conf.int = TRUE)## estimated baseline mean cumulative function (MCF) from a fitted modelpiecesMcf <- mcf(piecesFit)plot(piecesMcf, conf.int = TRUE, col = "blueviolet")## estimated MCF for given new datanewDat <- data.frame(x1 = rep(0, 2), group = c("Treat", "Contr"))splineMcf <- mcf(splineFit, newdata = newDat, groupName = "Group",                 groupLevels = c("Treatment", "Control"))plot(splineMcf, conf.int = TRUE, lty = c(1, 5))## example of further customization by ggplot2library(ggplot2)plot(splineMcf) +    geom_ribbon(aes(x = time, ymin = lower,                    ymax = upper, fill = Group),                data = splineMcf@MCF, alpha = 0.2) +    xlab("Days")

An S4 Class Representing a Fitted Model

Description

The classrateReg is an S4 class that represents a fitted model. ThefunctionrateReg produces objects of this class. See“Slots” for details.

Slots

call

Function call.

formula

Formula.

nObs

A positive integer

spline

A list.

estimates

A list.

control

A list.

start

A list.

na.action

A character vector (of length one).

xlevels

A list.

contrasts

A list.

convergCode

A nonnegative integer.

logL

A numeric value.

fisher

A numeric matrix.

See Also

rateReg

Show an object.

Description

S4 class methods that display objects produced from this package (similar toS3 classprint methods).

Usage

## S4 method for signature 'Recur'show(object)## S4 method for signature 'rateReg'show(object)## S4 method for signature 'summary.rateReg'show(object)## S4 method for signature 'summary.Recur'show(object)## S4 method for signature 'mcf.formula'show(object)## S4 method for signature 'mcf.rateReg'show(object)## S4 method for signature 'simEvent'show(object)## S4 method for signature 'mcfDiff'show(object)## S4 method for signature 'mcfDiff.test'show(object)

Arguments

Simulated Survival times or Recurrent Events

Description

The functionsimEvent generates simulated recurrent events orsurvival time (the first event time) from one stochastic process. ThefunctionsimEventData provides a simple wrapper that callssimEvent internally and collects the generated survival data orrecurrent events into a data frame. More examples are available in one ofthe package vignettes in addition to the function documentation.

Usage

simEvent(  z = 0,  zCoef = 1,  rho = 1,  rhoCoef = 1,  rhoMax = NULL,  origin = 0,  endTime = 3,  frailty = 1,  recurrent = TRUE,  interarrival = "rexp",  relativeRisk = c("exponential", "linear", "excess", "none"),  method = c("thinning", "inversion"),  arguments = list(),  ...)simEventData(nProcess = 1, z = 0, origin = 0, endTime = 3, frailty = 1, ...)

Arguments

Details

For each process, a time-invariant or time-varying baseline hazard rate(intensity) function of failure can be specified. Covariates and theircoefficients can be specified and incorporated by the specified relativerisk functions. The default is the exponential relative risk function, whichcorresponds to the Cox proportional hazard model (Cox 1972) for survivaldata or Andersen-Gill model (Andersen and Gill 1982) for recurrentevents. Other relative risk function can be specified through the argumentrelativeRisk. In addition, a frailty effect can be considered.Conditional on predictors (or covariates) and the unobserved frailty effect,the process is by default a Poisson process, where the interarrival timesbetween two successive arrivals/events follow exponential distribution. Ageneral renewal process can be specified throughinterarrival forother distributions of the interarrival times in addition to the exponentialdistribution.

The thinning method (Lewis and Shedler 1979) is applied for bounded hazardrate function by default. The inversion method (Cinlar 1975) is alsoavailable for possibly unbounded but integrable rate function over the giventime period. The inversion method will be used when the rate function may goto infinite and a warning will be thrown out if the thinning method isspecified originally.

For the covariatesz, the covariate coefficientszCoef, andthe baseline hazard rate functionrho, a function of time can bespecified for time-varying effect. The first argument of the input functionhas to be the time variable (not need to be named as "time" though). Otherarguments of the function can be specified through a named list inarguments, while the first argument should not be specified.

For the frailty effectfrailty, the starting pointorigin, andthe end point of the processendTime, functions that generate randomnumbers can be specified. An argumentn = 1 will be implicitlyspecified if the function has an argument namedn, which is designedfor those common functions generating random numbers fromstatspackage. Similar toz,zCoef, andrho, other argumentsof the function can be specified through a named list inarguments.

For time-varying covariates, the functionsimEventData assumescovariates can be observed only at event times and censoring times. Thus,covariate values are returned only at these time points. If we want otherobserved covariate values to be recorded, we may write a simple wrapperfunction forsimEvent similar tosimEventData.

Value

The functionsimEvent returns asimEvent S4 class object andthe functionsimEventData returns adata.frame.

References

Andersen, P. K., & Gill, R. D. (1982). Cox's regression model for countingprocesses: A large sample study.The annals of statistics, 10(4),1100–1120.

Cinlar, Erhan (1975).Introduction to stochastic processes. EnglewoodCliffs, NJ: Printice-Hall.

Cox, D. R. (1972). Regression models and life-tables.Journal of the Royal Statistical Society. Series B(Methodological), 34(2), 187–220.

Lewis, P. A., & G. S. Shedler. (1979). Simulation ofNonhomogeneous Poisson Processes by Thinning.Naval Research Logistics Quarterly,26(3), Wiley Online Library: 403–13.

Examples

library(reda)set.seed(123)### time-invariant covariates and coefficients## one processsimEvent(z = c(0.5, 1), zCoef = c(1, 0))simEvent(z = 1, zCoef = 0.5, recurrent = FALSE)## simulated datasimEventData(z = c(0.5, 1), zCoef = c(1, 0), endTime = 2)simEventData(z = cbind(rnorm(3), 1), zCoef = c(1, 0))simEventData(z = matrix(rnorm(5)), zCoef = 0.5, recurrent = FALSE)### time-varying covariates and time-varying coefficientszFun <- function(time, intercept) {    cbind(time / 10 + intercept, as.numeric(time > 1))}zCoefFun <- function(x, shift) {    cbind(sqrt(x + shift), 1)}simEvent(z = zFun, zCoef = zCoefFun,         arguments = list(z = list(intercept = 0.1),                          zCoef = list(shift = 0.1)))## same function of time for all processessimEventData(3, z = zFun, zCoef = zCoefFun,             arguments = list(z = list(intercept = 0.1),                              zCoef = list(shift = 0.1)))## same function within one process but different between processes## use quote function in the argumentssimDat <- simEventData(3, z = zFun, zCoef = zCoefFun,                       arguments = list(                           z = list(intercept = quote(rnorm(1) / 10)),                           zCoef = list(shift = 0.1)                       ))## check the intercept randomly generated,## which should be the same within each ID but different between IDs.unique(with(simDat, cbind(ID, intercept = round(X.1 - time / 10, 6))))### non-negative time-varying baseline hazard rate functionsimEvent(rho = function(timeVec) { sin(timeVec) + 1 })simEventData(3, origin = rnorm(3), endTime = rnorm(3, 5),             rho = function(timeVec) { sin(timeVec) + 1 })## specify other argumentssimEvent(z = c(rnorm(1), rbinom(1, 1, 0.5)) / 10,         rho = function(a, b) { sin(a + b) + 1 },         arguments = list(rho = list(b = 0.5)))simEventData(z = cbind(rnorm(3), rbinom(3, 1, 0.5)) / 10,             rho = function(a, b) { sin(a + b) + 1 },             arguments = list(rho = list(b = 0.5)))## quadratic B-splines with one internal knot at "time = 1"## (using function 'bSpline' from splines2 package)simEvent(rho = splines2::bSpline, rhoCoef = c(0.8, 0.5, 1, 0.6),         arguments = list(rho = list(degree = 2, knots = 1,                                     intercept = TRUE,                                     Boundary.knots = c(0, 3))))### frailty effect## Gamma distribution with mean onesimEvent(z = c(0.5, 1), zCoef = c(1, 0), frailty = rgamma,         arguments = list(frailty = list(shape = 2, scale = 0.5)))## lognormal with mean zero (on the log scale)set.seed(123)simEvent(z = c(0.5, 1), zCoef = c(1, 0), frailty = "rlnorm",         arguments = list(frailty = list(sdlog = 1)))## or equivalentlyset.seed(123)logNorm <- function(a) exp(rnorm(n = 1, mean = 0, sd = a))simEvent(z = c(0.5, 1), zCoef = c(1, 0), frailty = logNorm,         arguments = list(frailty = list(a = 1)))### renewal process## interarrival times following uniform distributionrUnif <- function(n, rate, min) runif(n, min, max = 2 / rate)simEvent(interarrival = rUnif,         arguments = list(interarrival = list(min = 0)))## interarrival times following Gamma distribution with scale oneset.seed(123)simEvent(interarrival = function(n, rate) rgamma(n, shape = 1 / rate))## or equivalentlyset.seed(123)simEvent(interarrival = function(rate) rgamma(n = 1, shape = 1 / rate))### relative risk functioinset.seed(123)simEvent(relativeRisk = "linear")## or equivalentlyrriskFun <- function(z, zCoef, intercept) {    as.numeric(z %*% zCoef) + intercept}set.seed(123)simEvent(relativeRisk = rriskFun,         arguments = list(relativeRisk = list(intercept = 1)))

An S4 Class for Simulated Recurrent Event or Survival Times

Description

An S4 class that represents the simulated recurrent event or survival timefrom one stochastic process. The functionsimEvent producesobjects of this class.

Slots

.Data

A numerical vector of possibly length zero.

call

A function call.

z

A list.

zCoef

A list.

rho

A list.

rhoCoef

A numerical vector.

frailty

A list.

origin

A list.

endTime

A list.

censoring

A list.

recurrent

A logical vector.

interarrival

A list.

relativeRisk

A list.

method

A character vector.

See Also

simEvent

Simulated Sample Dataset for Demonstration

Description

A simulated data frame with covariates namedID,time,event,group,x1,andgender, where

• ID: Subjects identification;

• time: Event or censoring time;

• event: Event indicator, 1 = event, 0 = censored;

• group: Treatment group indicator;

• x1: Continuous variable.

• gender: Gender of subjects.

Format

A data frame with 500 rows and 6 variables.

Details

The sample dataset is originally simulated by the thinningmethod developed by Lewis and Shedler (1979) andfurther processed for a better demonstration purpose.See Fu et al. (2016) for details also.

References

Lewis, P. A., & Shedler, G. S. (1979).Simulation of nonhomogeneous Poisson processes by thinning.Naval Research Logistics Quarterly, 26(3), 403–413.

Fu, H., Luo, J., & Qu, Y. (2016).Hypoglycemic events analysis via recurrent time-to-event (HEART) models.Journal Of Biopharmaceutical Statistics, 26(2), 280–298.

Summarize anRecur object

Description

Summarize anRecur object

Usage

## S4 method for signature 'Recur'summary(object, ...)

Arguments

Value

summary.Recur object.

Summarizing a Fitted Model

Description

Summary of estimated coefficients of covariates, rate function bases,and estimated rate parameter of frailty random variable, etc.

Usage

## S4 method for signature 'rateReg'summary(object, showCall = TRUE, showKnots = TRUE, ...)

Arguments

Details

summary,rateReg-method returns asummary.rateReg object,whose slots include

• covarCoef: Estimated covariate coefficients.

• frailtyPar: Estimated rate parameter of gamma frailty.

• baseRateCoef: Estimated coeffcients of baselinerate function.

For the meaning of other slots, seerateReg.

Value

summary.rateReg object

See Also

rateReg for model fitting;coef,rateReg-method for point estimates ofcovariate coefficients;confint,rateReg-method for confidence intervalsof covariate coeffcients;baseRate,rateReg-method for coefficients of baselinerate function.

Examples

## See examples given in function rateReg.

An S4 Class for Summarized Recur Object

Description

An S4 Class for Summarized Recur Object

Slots

call

A function call.

sampleSize

An integer representing the sample size (number ofsubjects).

reSize

An integer representing the number of recurrent events.

avgReSize

A numeric value representing the average number of recurrentevents per subject.

propTem

A numeric value representing the proportion of subjects havingterminal event.

medFU

A numeric value for median follow-up time.

medTem

A numeric value for median survival time of the terminalevents.

An S4 Class Representing Summary of a Fitted Model

Description

The classsummary.rateReg is an S4 class with selective slots ofrateReg object. See “Slots” for details. The functionsummary,rateReg-method produces objects of this class.

Slots

call

Function call.

spline

A character.

knots

A numeric vector.

Boundary.knots

A numeric vector.

covarCoef

A numeric matrix.

frailtyPar

A numeric matrix.

degree

A nonnegative integer.

baseRateCoef

A numeric matrix.

logL

A numeric value.

See Also

summary,rateReg-method

Valve Seats Dataset

Description

Valve seats wear out in certain diesel engines, each with 16 valve seats.The dataset served as an example of recurrence data in Nelson (1995),which consists of valve-seat replacements on 41 engines in a fleet.The covariates are namedID,Days, andNo., where

• ID: The engine number;

• Days: Engine age in days;

• No.: Event indicator, '1' for a valve-seat replacementand, '0' for the censoring age of an engine.

Format

A data frame with 89 rows and 3 variables.

References

Nelson, W. (1995), Confidence Limits for RecurrenceData-Applied to Cost or Number of Product Repairs,Technometrics,37, 147–157.