simEvent  R Documentation 
The function simEvent
generates simulated recurrent events or
survival time (the first event time) from one stochastic process. The
function simEventData
provides a simple wrapper that calls
simEvent
internally and collects the generated survival data or
recurrent events into a data frame. More examples are available in one of
the package vignettes in addition to the function documentation.
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, ...)
z 
Timeinvariant or timevarying covariates. The default value is

zCoef 
Timeinvariant or timevarying coefficients of covariates. The
default value is 
rho 
Baseline rate (or intensity) function for the Poisson process.
The default is 
rhoCoef 
Coefficients of baseline rate function. The default value is

rhoMax 
A positive number representing an upper bound of the underlying rate function (excluding the frailty term but including the covariate effect) for the thinning method. If this argument is left unspecified, the function will try to determine an upper bound internally. 
origin 
The time origin set to be 
endTime 
The end of followup time set to be 
frailty 
A positive number or a function for frailty effect. The
default value is 
recurrent 
A logical value with default value 
interarrival 
A function object for randomly generating (positive)
interarrival time between two successive arrivals/events. The default
value is 
relativeRisk 
Relateive risk function for incorporating the covariates
and the covariate coefficients into the intensity function. The
applicable choices include 
method 
A character string specifying the method for generating simulated recurrent or survival data. The default method is thinning method (Lewis and Shedler 1979). Another available option is the inversion method (Cinlar 1975). When the rate function may go to infinite, the inversion method is used and a warning will be thrown out if the thinning method is initially specified. 
arguments 
A list that consists of named lists for specifying other
arguments in the corresponding functions. For example, if a function of
time named 
... 
Additional arguements passed from function 
nProcess 
Number of stochastic processes. If missing, the value will
be the number of row of the specified matrix 
For each process, a timeinvariant or timevarying baseline hazard rate
(intensity) function of failure can be specified. Covariates and their
coefficients can be specified and incorporated by the specified relative
risk functions. The default is the exponential relative risk function, which
corresponds to the Cox proportional hazard model (Cox 1972) for survival
data or AndersenGill model (Andersen and Gill 1982) for recurrent
events. Other relative risk function can be specified through the argument
relativeRisk
. 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 times
between two successive arrivals/events follow exponential distribution. A
general renewal process can be specified through interarrival
for
other distributions of the interarrival times in addition to the exponential
distribution.
The thinning method (Lewis and Shedler 1979) is applied for bounded hazard rate function by default. The inversion method (Cinlar 1975) is also available for possibly unbounded but integrable rate function over the given time period. The inversion method will be used when the rate function may go to infinite and a warning will be thrown out if the thinning method is specified originally.
For the covariates z
, the covariate coefficients zCoef
, and
the baseline hazard rate function rho
, a function of time can be
specified for timevarying effect. The first argument of the input function
has to be the time variable (not need to be named as "time" though). Other
arguments of the function can be specified through a named list in
arguments
, while the first argument should not be specified.
For the frailty effect frailty
, the starting point origin
, and
the end point of the process endTime
, functions that generate random
numbers can be specified. An argument n = 1
will be implicitly
specified if the function has an argument named n
, which is designed
for those common functions generating random numbers from stats
package. Similar to z
, zCoef
, and rho
, other arguments
of the function can be specified through a named list in arguments
.
For timevarying covariates, the function simEventData
assumes
covariates can be observed only at event times and censoring times. Thus,
covariate values are returned only at these time points. If we want other
observed covariate values to be recorded, we may write a simple wrapper
function for simEvent
similar to simEventData
.
The function simEvent
returns a simEvent
S4 class object and
the function simEventData
returns a data.frame
.
Andersen, P. K., & Gill, R. D. (1982). Cox's regression model for counting processes: A large sample study. The annals of statistics, 10(4), 1100–1120.
Cinlar, Erhan (1975). Introduction to stochastic processes. Englewood Cliffs, NJ: PrinticeHall.
Cox, D. R. (1972). Regression models and lifetables. Journal of the Royal Statistical Society. Series B (Methodological), 34(2), 187–220.
Lewis, P. A., & G. S. Shedler. (1979). Simulation of Nonhomogeneous Poisson Processes by Thinning. Naval Research Logistics Quarterly, 26(3), Wiley Online Library: 403–13.
library(reda) set.seed(123) ### timeinvariant covariates and coefficients ## one process simEvent(z = c(0.5, 1), zCoef = c(1, 0)) simEvent(z = 1, zCoef = 0.5, recurrent = FALSE) ## simulated data simEventData(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) ### timevarying covariates and timevarying coefficients zFun < 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 processes simEventData(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 arguments simDat < 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)))) ### nonnegative timevarying baseline hazard rate function simEvent(rho = function(timeVec) { sin(timeVec) + 1 }) simEventData(3, origin = rnorm(3), endTime = rnorm(3, 5), rho = function(timeVec) { sin(timeVec) + 1 }) ## specify other arguments simEvent(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 Bsplines 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 one simEvent(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 equivalently set.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 distribution rUnif < 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 one set.seed(123) simEvent(interarrival = function(n, rate) rgamma(n, shape = 1 / rate)) ## or equivalently set.seed(123) simEvent(interarrival = function(rate) rgamma(n = 1, shape = 1 / rate)) ### relative risk functioin set.seed(123) simEvent(relativeRisk = "linear") ## or equivalently rriskFun < function(z, zCoef, intercept) { as.numeric(z %*% zCoef) + intercept } set.seed(123) simEvent(relativeRisk = rriskFun, arguments = list(relativeRisk = list(intercept = 1)))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.