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
)A list of 2 data.frames: one recording the requisite
longitudinal outcomes data, and one recording the time-to-event data.
the number of subjects to simulate data for.
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.
logical value: if TRUE then the joint model is fitted
with separate association, see Details.
the maximum number of (discrete) time points to simulate repeated longitudinal measurements at.
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.
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.
a vector of specifying the latent association parameter(s) for
the longitudinal outcome. It must be of length 1 if sepassoc =
FALSE.
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.
a numeric value specifying the residual standard error.
parameters controlling the failure rate. See Details.
logical: if TRUE, includes an independent censoring
time.
a scale (\(>0\)) parameter for an exponential distribution
used to simulate random censoring times for when censoring = TRUE.
logical: if TRUE, adds a truncation time for a
maximum event time in the case of model = "int" or model =
"intslope".
a truncation time for use when truncation = TRUE. For
model = "quad", trunctime is required, and defaults to
max(ntms) if not specified.
the step-size for the grid used to simulate event times when
model = "quad". Default is gridstep = 0.01. See Details.
Pete Philipson
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.
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.