*! 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 user-defined distribution:

survsim newvarname1 newvarname2 , maxtime(#|varname) [options]

See help survsim for more on simulating survival times in other settings.

Options

loghazard(string) user-defined log baseline hazard function; see details

hazard(string) user-defined baseline hazard function; see details

logchazard(string) user-defined log cumulative baseline hazard function; see details

chazard(string) user-defined baseline cumulative hazard function; 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

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

ltruncated(#|varname) left truncation time(s) (delayed entry); either a common number or a varname

nodes(#) number of Gauss-Legendre quadrature nodes, default 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 Markov multi-state model.
This help file centres on the user-defined distribution setting.

Survival times can be simulated from bespoke, user-defined [log] [cumulative] hazard functions. The function must be specified
in Mata code (using colon operators), with survival times generated using a combination of numerical integration and root
finding techniques. Baseline covariates can be included, with specified associated log hazard ratios. Time-dependent effects
can also be specified with a user-defined function of time. See Crowther and Lambert (2013) for more
details.

newvarname1 specifies the new variable name to contain the generated survival times. newvarname2 specifies the
new variable name to contain the generated event indicator.

Options

loghazard(string) is the user-defined log hazard function. The function can include:

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

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 hazard function.

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

See examples below.

hazard(string) is the user-defined baseline hazard function. See loghazard() for more details, and
examples below.

logchazard(string) is the user-defined log cumulative baseline hazard function. See loghazard() for
more details, and examples below.

chazard(string) is the user-defined baseline cumulative hazard function. See loghazard() for more
details, and examples below.

covariates(varname # [# ...] ...) defines baseline covariates to be included in the linear predictor of the
survival model, 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. If chazard() or logchazard() are used, then covariates() effects are additive
on the log cumulative hazard scale.

tde(varname # [# ...] ...) creates non-proportional hazards by interacting covariates with a function of time,
defined by tdefunction(), on the appropriate log hazard or log cumulative 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 tdefunction(). To circumvent this, you can directly specify them in your user function.

tdefunction(string) defines the function of time to which covariates specified in tde() are interacted
with, to create time-dependent effects. 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 hazard function.

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

maxtime(#|varname) specifies the 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.

ltruncated(#|varname) specifies the left truncated/delayed entry time(s). Either a common time # can be
specified for all observations, or observation specific left truncation times can be specified by using a varname.

nodes(#) defines the number of Gauss-Legendre quadrature points used to evaluate the cumulative hazard function
when loghazard() or hazard() is specified. 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 from a user-defined loghazard() or hazard() function, numerical integration is used to
evaluate the 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 integration points, through use of the nodes() option. Any survsim call that requires Brent's method
is executed under a tolerance of 0, to ensure accurracy of the simulated survival times.

Always set seed, to ensure reproducibility.

Examples

Simulate 1000 observations and generate a binary treatment group:

. set obs 1000

. gen trt = rbinomial(1,0.5)

Generate times from user-defined log hazard function:

. survsim stime1 died1, loghazard(-1 :+ 0.02:* :- 0.03:*:^2 :+ 0.005:*:^3) maxtime(1.5)

Generate times from user-defined log hazard function with diminishing treatment effect:

. survsim stime2 died2, loghazard(-1 :+ 0.02:* :- 0.03:*:^2 :+ 0.005:*:^3) covariates(trt -0.5) tde(trt 0.03) maxtime(1.5)

Generate survival times from a joint longitudinal-survival model:

. clear

. set obs 1000

. gen trt = rbinomial(1,0.5)

. gen age = rnormal(65,12)

Define the association between the biomarker and survival

. local alpha = 0.25

Generate the random intercept and random slopes for the longitudinal submodel

. gen b0 = rnormal(0,1)

. gen b1 = rnormal(1,0.5)

Generate survival times from an exponential baseline hazard

. survsim st1 event, loghazard(`=log(0.1)' :+ `alpha' :* (b0 :+ b1 :* )) maxtime(5) covariates(trt -0.5 age 0.02)

Generate observed biomarker values at times 0, 1, 2, 3 , 4 years

. gen id = _n

. expand 5

. bys id: gen meastime = _n-1

Remove observations after event or censoring time

. bys id: drop if meastime >= st1

Generate observed biomarker values incorporating measurement error

. gen response = b0 + b1*meastime + rnormal(0,0.5)

For more examples please see Crowther and Lambert (2013).

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.

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.