Simulate data from a joint model
This function simulates longitudinal and time-to-event data from a joint model.
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 )
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
|
sepassoc |
logical value: if |
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 |
gamma |
a vector of specifying the latent association parameter(s) for
the longitudinal outcome. It must be of length 1 if |
sigu |
a positive-definite matrix specifying the variance-covariance
matrix. If |
vare |
a numeric value specifying the residual standard error. |
theta0, theta1 |
parameters controlling the failure rate. See Details. |
censoring |
logical: if |
censlam |
a scale (>0) parameter for an exponential distribution
used to simulate random censoring times for when |
truncation |
logical: if |
trunctime |
a truncation time for use when |
gridstep |
the step-size for the grid used to simulate event times when
|
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
λ_0(t) = \exp{θ_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
λ_0(t) = \exp{θ_0 + θ_1 t}.
In the usual representation of the Gompertz distribution, θ_1 is the shape parameter, and the scale parameter is equivalent to \exp(θ_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} λ(t) dt = \lim_{dt \rightarrow 0} P[t ≤ T ≤ t + dt | T ≥ 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.
A list of 2 data.frame
s: one recording the requisite
longitudinal outcomes data, and one recording the time-to-event data.
Pete Philipson (pete.philipson@northumbria.ac.uk)
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.
simjoint(10, sepassoc = TRUE)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.