simjoint {joineR}

Simulate data from a joint model

Description

This function simulates longitudinal and time-to-event data from a joint model.

Usage

simjoint(
  n = 500,
  model = c("intslope", "int", "quad"),
  sepassoc = FALSE,
  ntms = 5,
  b1 = c(1, 1, 1, 1),
  b2 = c(1, 1),
  gamma = c(1, 0.1),
  sigu,
  vare = 0.01,
  theta0 = -3,
  theta1 = 1,
  censoring = TRUE,
  censlam = exp(-3),
  truncation = FALSE,
  trunctime = max(ntms),
  gridstep = 0.01
)

Arguments

n

the number of subjects to simulate data for.

model

a character string specifying the type of latent association. This defaults to the intercept and slope version as seen in Wulfsohn and Tsiatis (1997). For association via the random intercept only, choose model = "int", whereas for a quadratic association, use model = "quad". Computing times are commensurate with the type of association structure chosen.

sepassoc

logical value: if TRUE then the joint model is fitted with separate association, see Details.

ntms

the maximum number of (discrete) time points to simulate repeated longitudinal measurements at.

b1

a vector specifying the coefficients of the fixed effects in the longitudinal sub-model. The order in each row is intercept, a continuous covariate, covariate, a binary covariate, the measurement time.

b2

a vector of length = 2 specifying the coefficients for the time-to-event baseline covariates, in the order of a continuous covariate and a binary covariate.

gamma

a vector of specifying the latent association parameter(s) for the longitudinal outcome. It must be of length 1 if sepassoc = FALSE.

sigu

a positive-definite matrix specifying the variance-covariance matrix. If model = "int", the matrix has dimension dim = c(1, 1); if model = "intslope", the matrix has dimension dim = c(2, 2); else if model = "quad", the matrix has dimension dim = c(3, 3). If D = NULL (default), an identity matrix is assumed.

vare

a numeric value specifying the residual standard error.

theta0, theta1

parameters controlling the failure rate. See Details.

censoring

logical: if TRUE, includes an independent censoring time.

censlam

a scale (>0) parameter for an exponential distribution used to simulate random censoring times for when censoring = TRUE.

truncation

logical: if TRUE, adds a truncation time for a maximum event time in the case of model = "int" or model = "intslope".

trunctime

a truncation time for use when truncation = TRUE. For model = "quad", trunctime is required, and defaults to max(ntms) if not specified.

gridstep

the step-size for the grid used to simulate event times when model = "quad". Default is gridstep = 0.01. See Details.

Details

The function simjoint simulates data from a joint model, similar to that performed in Henderson et al. (2000). It works by first simulating longitudinal data for all possible follow-up times using random draws for the multivariate Gaussian random effects and residual error terms. Data can be simulated assuming either random-intercepts only (model = "int") in each of the longitudinal sub-models; random-intercepts and random-slopes (model = "intslope"); or quadratic random effects structures (model = "quad"). The failure times are simulated from proportional hazards time-to-event models, using the following methodologies:

model = "int"

The baseline hazard function is specified to be an exponential distribution with

\lambda_0(t) = \exp{\theta_0}.

Simulation is conditional on known time-independent effects, and the methodology of Bender et al. (2005) is used to simulate the failure time.

model = "intslope"

The baseline hazard function is specified to be a Gompertz distribution with

\lambda_0(t) = \exp{\theta_0 + \theta_1 t}.

In the usual representation of the Gompertz distribution, \theta_1 is the shape parameter, and the scale parameter is equivalent to \exp(\theta_0). Simulation is conditional on on a predictable (linear) time-varying process, and the methodology of Austin (2012) is used to simulate the failure time.

model="quad"

The baseline hazard function is specified as per model="intslope". The integration technique used for the above two cases is complex in quadratic (and higher order) models, therefore we use a different approach. We note that hazard function can be written as

\lim_{dt \rightarrow 0} \lambda(t) dt = \lim_{dt \rightarrow 0} P[t \le T \le t + dt | T \ge t].

In the simulation routine the parameter gridstep acts as dt in that we choose a nominally small value, which multiplies the hazard and this scaled hazard is equivalent to the probability of having an event in the interval (t, t + dt), or equivalently (t, t + gridstep). A vector of possible times is set up for each individual, ranging from 0 to trunctime in increments of dt (or gridstep). The failure probability at each time is compared to an independent U(0, 1) draw, and if the probability does not exceed the random draw then the survival time is set as trunctime, otherwise it is the generated time from the vector of candidate times. The minimum of these candidate times (i.e. when we deem the event to have first happened) is taken as the survival time.

Value

A list of 2 data.frames: one recording the requisite longitudinal outcomes data, and one recording the time-to-event data.

Author(s)

Pete Philipson

References

Austin PC. Generating survival times to simulate Cox proportional hazards models with time-varying covariates. Stat Med. 2012; 31(29): 3946-3958.

Bender R, Augustin T, Blettner M. Generating survival times to simulate Cox proportional hazards models. Stat Med. 2005; 24: 1713-1723.

Henderson R, Diggle PJ, Dobson A. Joint modelling of longitudinal measurements and event time data. Biostatistics. 2000; 1(4): 465-480.

Examples

simjoint(10, sepassoc = TRUE)