simEvent: Simulated Survival times or Recurrent Events

In reda: Recurrent Event Data Analysis

Simulated Survival times or Recurrent Events

Description

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.

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

z	Time-invariant or time-varying covariates. The default value is 0 for no covariate effect. This argument should be a numeric vector for time-invariant covariates or a function of times that returns a numeric matrix for time-varying covariates, where each row represents the covariate vector at one perticular time point.
zCoef	Time-invariant or time-varying coefficients of covariates. The default value is 1. Similar to the argument z, this argument should be a numeric vector for time-invariant coefficients or a function of times that returns a numeric matrix for time-varying coefficients, where each row represents the coefficient vector at one perticular time point. The dimension of the z and zCoef (either specified or generated) has to match with each other.
rho	Baseline rate (or intensity) function for the Poisson process. The default is 1 for a homogenous process of unit intensity. This argument can be either a non-negative numeric value for a homogenous process or a function of times for a non-homogenous process. In the latter case, the function should be able to take a vector of time points and return a numerical matrix (or vector) with each row representing the baseline hazard rate vector (or scalar value) at each time point.
rhoCoef	Coefficients of baseline rate function. The default value is 1. It can be useful when rho is a function generating spline bases.
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 0 by default. It should be either a numeric value less than endTime or a function that returns such a numeric value.
endTime	The end of follow-up time set to be 3 by default. Similar to origin, endTime should be either a numeric value greater than origin or a function that returns such a numeric value.
frailty	A positive number or a function for frailty effect. The default value is 1 for no frailty effect. Other positive value can be specified directly for a shared frailty effect within a cluster. Similar to z, zCoef, and rho, a function can be specified for other distribution of the frailty effect. The specified function should randomly return a positive numeric value. The functions that generate random numbers following a certain distribution from stats package can be directly used. The arguments of the function can be specified through a list named frailty in arguments. For example, if we consider Gamma distribution with mean one as the distribution of frailty effect, we may specify frailty = "rgamma". The shape and scale parameter needs to be specified through a list named frailty in arguments, such as arguments = list(frailty = list(shape = 2, scale = 0.5)).
recurrent	A logical value with default value TRUE indicating whether to generate recurrent event data or survival data.
interarrival	A function object for randomly generating (positive) interarrival time between two successive arrivals/events. The default value is "rexp" (i.e., function stats::rexp) for generating interarrival times following exponential distribution, which leads to a Poisson process. If the assumption of exponential interarrival times cannot be justified, we may consider a renewal process, (a generalization of Poisson process), in which interarrival times between events independently follows an identical distribution. A customized function can be specified in this case. It must have at least one argument named rate for the expected number of arrivals/events in unit time and returns one positive numerical value. If the function contains an argument named n, it is assumed that the function returns n interarrival times in one function call to possibly speed up the random number generation procedure. Other arguments can be specified through a named list inside arguments.
relativeRisk	Relateive risk function for incorporating the covariates and the covariate coefficients into the intensity function. The applicable choices include exponential (the default) for the regular Cox model or Andersen-Gill model, linear for linear model (including an intercept term), excess for excess model, and none for not incorporating the covariates through a relative risk function. A customized function can be specified. The specified function must have at least one argument named z for the covariate vector and another argument named zCoef for covariate coefficient vector. The function should return a numeric value for given z vector and zCoef vector. Other arguments can be specified through a named list inside arguments.
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 foo with two arguments, x (for time) and y, is specified for the time-varying covariates, the value of its second argument, y, can be specified by arguments = list(z = list(y = 1). A partial matching on names is not allowed to avoid possible misspecification. The input arguments will be evaluated within function simEvent, which can be useful for randomly setting function parameters for each process in function simEventData. See examples and vignettes for details.
...	Additional arguements passed from function simEventData to fucntion simEvent. For function simEvent, ... is not used.
nProcess	Number of stochastic processes. If missing, the value will be the number of row of the specified matrix z. Otherwise, a positive number should be speicified.

Details

For each process, a time-invariant or time-varying 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 Andersen-Gill 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 time-varying 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 time-varying 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.

Value

The function simEvent returns a simEvent S4 class object and the function simEventData returns a data.frame.

References

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: 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 of Nonhomogeneous 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 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)


### time-varying covariates and time-varying 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))))


### non-negative time-varying 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 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 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)))

reda documentation built on July 9, 2022, 1:06 a.m.