*! version 1.0.0

Title

survsim -- Simulate survival data from a parametric distribution,
a user-defined distribution, from a fitted merlin model, from a cause-specific
hazards competing risks model, or from a general multi-state model

Syntax

Syntax for simulating survival times from a multi-state model:

survsim timestub statestub eventstub , hazard1(haz_options)
hazard2(haz_options) maxtime(#|varname)
[transmatrix(name) hazard3(haz_options)
... options]

Each hazard#(haz_options) represents a transition-specific hazard function, where the number # is used to index the transition in the transition matrix, defined in transmatrix().

Transition-specific hazard#() options

distribution(exponential) exponential survival distribution

distribution(gompertz) Gompertz survival distribution

distribution(weibull) Weibull survival distribution

lambda(#) scale parameter for exponential, Weibull or Gompertz

gamma(#) shape parameter for Weibull or Gompertz

user(function) user-defined hazard function, written in Mata code; see details

covariates(varname # [# ...] ...) baseline covariates

tde(varname # [# ...] ...) time-dependent effects

tdefunction(string) function of time to interact with covariates specified in tde(); see details

reset clock-reset timescale, i.e. the timescale is reset to zero on state entry

Options

transmatrix(name) a transition matrix defining possible transitions between states

maxtime(#|varname) right censoring time(s); either a common number or a varname

startstate(#|varname) starting state for observations; either a common number or a varname

ltruncated(#|varname) starting time for observations, i.e. left truncation time(s); either a common number or a varname

nodes(#) number of Gauss-Legendre quadrature nodes, default nodes(30)

Description

survsim simulates survival data from a parametric distribution, a user-defined distribution, from a fitted
merlin model, from a cause-specific hazards competing risks model, or from a multi-state model.
This help file centres on simulating from a competing risks or multi-state model, specified using transition-specific
hazard functions.

Each transition-specific hazard function can be specified as a standard parametric model, such as the exponential, Weibull or
Gompertz distributions, or a user-defined bespoke hazard function written in Mata code. Baseline covariates can be included
in each transition-specific hazard, with specified associated log hazard ratios. Non-proportional hazards can also be included.
Transition-specific hazard models can be as similar, or as different, as required. Up to 50 transition-specific hazard
functions may be specified.

From each potential starting state, event times are simulated from the overall, total hazard function, using the method of
Beyersmann et al. (2009), i.e. the sum of all the transition-specific hazards out of the starting
state. This is implemented using a combination of numerical integration and root finding techniques using the algorithm
developed in Crowther and Lambert (2013).

timestub specifies the new variable name stub to contain the generated survival times for each transition.
statestub specifies the new variable name stub to contain the generated state identifiers, i.e which state each observation
has transitioned to/is in, at the associated times. eventstub specifies the new variable name stub to contain the
event indicators, i.e. whether an event or right-censoring occurred at the associated times. When observations have
reached an absorbing state, or reached maxtime(), they will have missing observations in any further generated variables.
When all observations have reached an absorbing state, or are right-censored, then the simulation will cease and the generated
variables are returned.

Transition-specific hazard#() options

distribution(string) specifies the parametric survival distribution to use, including exponential,
gompertz or weibull.

lambda(#) defines the scale parameter for a exponential/Weibull/Gompertz distribution(s).

gamma(#) defines the shape parameter for a Weibull/Gompertz parametric distribution(s).

user(function) defines a custom hazard function. The function can include:

which denotes the main timescale, measured on the time since starting state, startstate(),
timescale (which may be ltruncated())

which denotes the time of entry to the current state of the associated transition hazard, measured on the
time since initial starting state timescale

varname which denotes a variable in your dataset

+-*/^ standard Mata mathematical operators, using colon notation i.e. 2 :* , see
help [M-2] op_colon. Colon operators must be used as behind the scenes, gets replaced by an
_N x nodes() matrix when numerically integrating the transition hazard function.

mata_function any Mata function, e.g. log() and exp()

For example, dist(weibull) lambda(0.1) gamma(1.2) is equivalent to user(0.1:*1.2:*:^(1.2:-1))

covariates(varname # ...) defines baseline covariates to be included in the linear predictor of the
transition-specific hazard function, along with the value of the corresponding coefficient. For example, a treatent variable coded 0/1 can be
included, with a log hazard ratio of 0.5, by covariates(treat 0.5). Variable treat must be in the dataset before
survsim is called.

tde(varname # ...) creates non-proportional hazards by interacting covariates with a function of time.
Covariates are interacted with tdefunction(), on the log hazard scale. Values should be entered as tde(trt 0.5),
for example. Multiple time-dependent effects can be specified, but they will all be interacted with the same function of time.

tdefunction(string) defines the function of time to which covariates specified in tde() are interacted
with, to create time-dependent effects in the transition-specific hazard function. The default is , i.e. linear time.
The function can include:

which denotes the main timescale, measured on the time since starting state, startstate(),
timescale (which may be ltruncated())

+-*/^ standard Mata mathematical operators, using colon notation i.e. 2 :* , see
help [M-2] op_colon. Colon operators must be used as behind the scenes, gets replaced by an
_N x nodes() matrix when numerically integrating the transition hazard function.

mata_function any Mata function, e.g. log() and exp()

reset specifies that this transition model is on a clock-reset timescale. The timescale is reset to 0 on entry,
i.e. the timescale for this transition is measured on a time since state entry timescale, rather than the default
clock-forward. If you specify a user() function with reset, then survsim will replace any occurences of
with -, including those specified in tdefunction().

Options

transmatrix(matname) specifies the transition matrix which governs the multi-state model. Transitions must be
numbered as an increasing sequence of integers from 1,...,K, from left to right, top to bottom of the matrix.
Reversible transitions are allowed. If transmatrix() is not specified, a competing risks model is assumed.

maxtime(#|varname) specifies right censoring time(s). Either a common maximum follow-up time # can be
specified for all observations, or observation specific censoring times can be specified by using a varname.

startstate(#|varname) specifies the state(s) in which observations begin. Either a common state # can
be specified for all observations, or observation specific starting states can be specified by using a varname.
Default is startstate(1).

ltruncated(#|varname) specifies left truncated/delayed entry time(s), which is the time(s) at which observations
start in the initial starting state(s). Either a common time # can be specified for all observations, or observation
specific left truncation times can be specified by using a varname. Default is ltruncated(0).

nodes(#) defines the number of Gauss-Legendre quadrature points used to evaluate the total cumulative hazard
function for each potential next transition. To simulate survival times from such a function, a combination of numerical
integration and root-finding is used. The default is nodes(30).

Remarks

When simulating multi-state event data using survsim, numerical quadrature is used to calculate the total
cumulative hazard function, within iterations of Brent's univariate root finder. As with all model frameworks
which use numerical integration, it is important to assess the stability of the simulated survival times with an increasing
number of quadrature nodes. Any survsim call that requires Brent's method is executed under a tolerance of 0, to
ensure accurracy of the simulated survival times.

Note that even if reset transitions are used in the data-generating model, the returned simulated event times
are still on the main timescale, time since initial startstate() (which could be ltruncated()).

Always set seed, to ensure reproducibility.

Examples

Example 1: Simulate from a competing risks model

We'll simulate 1000 observations, and generate a binary treatment group indicator, remembering to set seed first.

. set obs 1000
. set seed 9865
. gen trt = runiform()>0.5

Simulate from a competing risk model with 2 competing events. The first cause-specific hazard has a Weibull distribution,
with no covariates. The second cause-specific hazard model has an exponential distribution, with a beneficial treatment effect.
Right censoring applied at 10 years.

. survsim time state event , hazard1(dist(weibull) lambda(0.1) gamma(1.2)) ///
hazard2(dist(exponential) lambda(0.01) covariates(trt -0.5)) ///
maxtime(10)
variables time0 to time1 created
variables state0 to state1 created
variables event1 to event1 created

Example 2: Simulate from an illness-death model

We first define the transition matrix for an illness-death model. It has three states:

State 1 - A "healthy" state. Observations can move from state 1 to state 2 or 3.

State 2 - An intermediate "illness" state. Observations can come from state 1, and move on to state 3.

State 3 - An absorbing "death" state. Observations can come from state 1 or 2, but not leave.

This gives us three potential transitions between states:

Transition 1 - State 1 -> State 2

Transition 2 - State 1 -> State 3

Transition 3 - State 2 -> State 3

which is defined by the following matrix:

. matrix tmat = (.,1,2\.,.,3\.,.,.)

The key is to think of the column/row numbers as the states, and the elements of the matrix as the transition numbers.
Any transitions indexed with a missing value . means that the transition between the row state and the column state
is not possible. Let's make it obvious, sticking with our "healthy", "ill" and "dead" names for the states:

. mat colnames tmat = "healthy" "ill" "dead"
. mat rownames tmat = "healthy" "ill" "dead"
. mat list tmat

tmat[3,3]
healthy ill dead
healthy . 1 2
ill . . 3
dead . . .

Now we've defined the transition matrix, we can use survsim to simulate some data. We'll simulate 1000 observations, and
generate a binary treatment group indicator, remembering to set seed first.

. set obs 1000
. set seed 9865
. gen trt = runiform()>0.5

The first transition-specific hazard has a user defined baseline hazard function, with a harmful treatment effect.
The second transition-specific hazard model has a Weibull distribition, with a beneficial treatment effect. The
third transition-specific hazard has a user-defined baseline hazard function, with an initially beneficial treatment
effect that reduces linearly with respect to log time. Right censoring is applied at 3 years.

. survsim time state event, transmatrix(tmat) ///
hazard1(user(exp(-2 :+ 0.2:* log() :+ 0.1:*)) covariates(trt 0.1)) ///
hazard2(dist(weibull) lambda(0.01) gamma(1.3) covariates(trt -0.5)) ///
hazard3(user(0.1 :* :^ 1.5) covariates(trt -0.5) tde(trt 0.1) ///
tdefunction(log())) ///
maxtime(3)
variables time0 to time2 created
variables state0 to state2 created
variables event1 to event2 created

survsim creates variables storing the times at which states were entered, with the associated state number, and whether
an event or right-censoring occured. It begins by creating the 0 variables, which represents the time at which
observatations entered the inital state, time0, and the associated state number, state0. As ltruncated()
and startstate() were not specified, all observations are assumed to start in state 1 at time 0. Subsequent transitions
are simulated until all observations have either entered an absorbing state, or are right-censored at their maxtime().
For simplicity, I will assume time is measured in years. We can see what survsim has created:

. list if inlist(_n,1,4,16,112)

+----------------------------------------------------------------------------------+
| trt time0 state0 time1 state1 event1 time2 state2 event2 |
|----------------------------------------------------------------------------------|
1. | 0 0 1 3 1 0 . . . |
4. | 1 0 1 .95636156 2 1 3 2 0 |
16. | 0 0 1 1.0755764 2 1 2.4401409 3 1 |
112. | 1 0 1 2.3290322 3 1 . . . |
+----------------------------------------------------------------------------------+

All observations start initially in state 1 at time 0, which are stored in state0 and time0, respectively. Then,

Observation 1 is right-censored at 3 years, remaining in state 1

Observation 4 moves to state 2 at 0.956 years, and is subsequently right-censored at 3 years, still in
state 2

Observation 16 moves to state 2 at 1.076 years, and then moves to state 3 at 2.440 years. Since state 3
is an absorbing state, there are no further transitions.

Observation 112 moves to state 3 at 2.329 years. Again, since state 3 is absorbing, there are no further
transitions

Author

Michael J. Crowther

Red Door Analytics

E-mail: michael@reddooranalytics.se

Please report any errors you may find.

References

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

Beyersmann J, Latouche A, Buchholz A and Schumacher M. Simulating competing risks data in survival analysis. Statistics in Medicine 2009;28:956-971.

Crowther MJ and Lambert PC. Simulating complex survival data. The Stata Journal 2012;12(4):674-687.

Crowther MJ and Lambert PC. Simulating biologically plausible complex survival data. Statistics in Medicine 2013;32(23):4118-4134.

Jann, B. 2005. moremata: Stata module (Mata) to provide various functions. Available from http://ideas.repec.org/c/boc/bocode/s455001.html.