*! version 1.4.9 22Jun2022
}help stmt
also see: stmt_postestimation

Title

stmt -- Modelling multiple timescales using flexible parametric survival models on the log hazard scale

Syntax

stmt [varlist] [if] [in] [, time1(sub-options) time2(sub-options) time3(sub-options) options]

Multiple timescale options
time1(sub-options) sub-options for timescale 1

time2(sub-options) sub-options for timescale 2

time3(sub-options) sub-options for timescale 3

timeint(int_options) specifies two-way timescale interactions

timeintknots(int_options) specifies internal knots for timescale interactions

timeintbknots(int_options) specifies boundary knots for timescale interactions

Timescale-specific sub-options
bknots(knots_list) boundary knots for the timescale specified

bknotstvc(knots_list) boundary knots for covariate-interactions with timescale

df(#) degrees of freedom for the timescale specified

dftvc(df_list) degrees of freedom for each interaction with the timescale

indicator(varname) specifies which observations are to have the second or third timescale

knots(knots_list) knot locations for timescale specified

knotstvc(knots_list) knot locations for interactions with the timescale

knscale(scale) scale for user-defined knots (default scale is time)

logtoff smooth the timescale over time (default is log time)

start(varname) the difference between the timescale and the timescale specified stset

tvc(varlist) varlist for covariate-timescale interactions

Estimation options
noconstant suppress constant term

nodes(#) specify the number of nodes used in Gauss-Legendre quadrature numerical
integration of the hazard function (default is 30)

noorthog do not use orthogonal transformation of splines variables

Reporting options
nohr do not report hazard ratios

verbose verbose output details the process of the stmt program

Maximization options
from(matrix) initial coefficient values stored in a matrix

inith(varname) initial hazard values

maximize_options control the maximization process; seldom used

You must stset your data before using stmt; see [ST] stset.

fweights, iweights, and pweights may be specified using stset; [ST] stset.

Description

stmt fits flexible parametric survival models (Royston-Parmar models)
on the log hazard scale where multiple timescales are allowed

See [ST] stpm2 for flexible parametric survival models on the log cumulative hazard scale, and [ST] strcs
for flexible parametric survival models on the log hazard scale.

Options

Multiple timescale options

time1(sub-options) time1() contains sub-options for the first timescale
(see below for a list of sub-options). The first timescale is always specified using the
stset command.

time2(sub-options) time2() contains sub-options for the second timescale
(see below for a list of sub-options). The second and third timescales are functions
of the first timescale; the difference between the second timescale and the first timescale
is specified in the start sub-option.

time3(sub-options) time3() contains sub-options for the third timescale
(see below for a list of sub-options). The second and third timescales are functions
of the first timescale; the difference between the third timescale and the first timescale
is specified in the start sub-option.

timeint(int_list) int_list contains details for two-way timescale interactions.
The following syntax applies:

timeint(t1:t2 2:4)

This will create an interaction between timescale 1 and timescale 2, as specified in the
time1() and time2() options where a restricted cubic spline function of timescale 1 with
2 degrees of freedom is interacted with a restricted cubic spline function of timescale 2
with 4 degrees of freedom. The space separates the specified timescales and their degrees
of freedom.Additional timescale interactions can be added using |
in the following way: timeint(t1:t2 2:4 | t1:t3 2:2).

timeintknots(int_list) int_list contains the internal knots for two-way timescale interactions.
The following syntax applies:

timeint(2 5 : 50 60 70)

This will create restricted cubic splines for the first timescale specified in the timeint()
option with internal knots at 2 and 5 (3 degrees of freedom), and restricted cubic splines for the
second timescale specified in the timeint() option at 50, 60 and 70 (4 degrees of freedom),
and then interact the spline terms together. The spaces separate the internal knot locations, whereas
the colon separates the timescales included in timeint(), i.e. the internal knots for the first
timescale in timeint() are specified prior to the colon, and the internal knots for the second
timescale in timeint() come after the colon. Internal knots for additional timescale interactions
can be added using | in a similar way to the timeint() option, for example:

timeintknots(2 5 : 50 60 70 | 2 3 5 : 2)

Note that knots should be specified on the untransformed
time scale.

timeintbknots(int_list) int_list contains the boundary knots for two-way timescale interactions.
The following syntax applies:

timeint(0 7 : 25 95)

This will create restricted cubic splines for the first timescale specified in timeint()
with boundary knots at 0 and 7, and restricted cubic splines for the second timescale
specified in timeint() with boundaries at 25 and 95, and then interact the spline terms
together. The spaces separate the boundary knot locations, whereas
the colon separates the timescales included in timeint(), i.e. the boundary knots for the first
timescale in timeint() are specified prior to the colon, and the boundary knots for the second
timescale in timeint() come after the colon. Boundary knots for additional timescale interactions can be added using | in
a similar way to the timeint() option, for example:

timeintbknots(0 7 : 25 95 | 0 7 : 1 12)

Note that both the timeint() options and timeintknots() options can be used with the
timeintbknots() option and that knots in the timeintbknots() should be specified
on the untransformed time scale.

Timescale-specific sub-options

bknots(knots_list) knots_list is a two-element numlist giving
the boundary knots. By default these are located at the minimum and maximum
of the uncensored survival times. They are specified on the scale defined
by knscale().

bknotstvc(knots_list) knots_list gives the boundary knots for
any restricted cubic splines created for the timescale, when including an interaction
between a covariate and the timescales. By default these are the same as for the bknots()
option. They are specified on the scale defined by knscale().

For example,

bknotstvc(x1 0.01 10 x2 0.01 8)

df(#) specifies the degrees of freedom for the restricted
cubic spline function used for the baseline function; the number of degrees of
freedom does not include the constant term. # must be between
1 and 10. The knots() option is not applicable if the df() option
is specified. The knots are placed at equally spaced centiles of the distribution
of the uncensored log survival times. For example, for df(5) knots are
placed at the 20th, 40th, 60th and 80th centiles of the distribution of the
uncensored log survival times. Note that these are interior knots and
there are also boundary knots placed at the minimum and maximum of the distribution
of uncensored survival times.

dftvc(df_list) gives the degrees of freedom used when including interactions
between the timescale and covariates in df_list.
With 1 degree of freedom a linear effect is fitted.
If there is more than one interaction and different degrees of freedom
are requested, then the following syntax applies:

dftvc(x1:3 x2:2 1)

This will use 3 degrees of freedom for x1, 2 degrees of freedom for
x2 and 1 degree of freedom for all remaining interactions between covariates
and the timescale.

indicator(varname) specifies an indicator variable which indicates which
observations have more than one timescale. The indicator variable should be
coded 0 for those observations who did not have the second/third timescale,
and 1 for those who did.

knots(# [# ...]) specifies knot locations for the baseline distribution
function, as opposed to the default locations set by df(). Note that
the locations of the knots are placed on the scale defined by knscale.
However, the scale used by the restricted cubic spline function is always
log time unless the logtoff option is specified. Default knot positions are determined by the df() option.

knotstvc(knots_list) defines numlist knots_list as the location
of the interior knots used when recalculating the restricted cubic splines for the
specified timescale-covariate interaction. If different knots
are required for different interactions, the option is
specified, for example, as follows:

knotstvc(x1 1 2 3 x2 1.5 3.5)

knscale(scale) sets the scale on which user-defined knots are specified.
knscale(time) denotes the original time scale, knscale(log) the
log time scale and knscale(centile) specifies that the knots
are taken to be centile positions in the distribution of the uncensored
survival times, or log survival times depending on whether the logtoff option
is specified. The default is knscale(time).

logtoff smooths the timescale over time using restricted cubic splines.
By default, smoothing is over log time.

start(varname) specifies the difference between the timescale specified in stset
and the timescale of interest. For example, if the first timescale of interest (t1) is time
since diagnosis, and the second timescale (t2) is attained age, attained age is equal to
time since diagnosis plus the age at diagnosis; i.e., t2=t1 + age at diagnosis. Thus in this
example, varname would be a variable which contains the age at diagnosis. Not for use when
using time1() since this timescale is specified when using the stset command.

tvc(varlist) gives the name of the variables that are to be included as part
of an interaction with the specified timescale. The degrees of freedom are specified
using the dftvc() option.

Estimation options

noconstant;
see [ST] estimation options.

nodes(#) specifies the number of nodes to be used in Gauss-Legendre quadrature
numerical integration when calculating the estimated cumulative hazard function
from the estimated hazard function. The default number of nodes is 30. Changing the
number of nodes may be useful if there are convergence problems. Too few nodes
may result in a poor approximation involved in the numerical integration. Analyses should
be performed to ensure the results are not sensitive to the number of nodes.

noorthog suppresses orthogonal transformation of spline variables.

Reporting

nohr reports the coefficients instead of hazard ratios.

verbose verbose output details the process of the stmt program.

Maximization options

from(matrix) defines the matrix containing initial coefficient estimates to be used in maximum
likelihood estimation. By default stmt estimates initial hazard estimates by fitting a model
on the log cumulative hazard scale using the stpm2 command.

inith(varname) defines initial hazard estimates to be used in maximum likelihood
estimation. By default stmt estimates initial hazard estimates by fitting a model on the
log cumulative hazard scale using the stpm2 command.

maximize_options; difficult, technique(algorithm_spec),
iterate(#), [no]log, trace, gradient,
showstep, hessian, shownrtolerance, tolerance(#),
ltolerance(#) gtolerance(#), nrtolerance(#),
nonrtolerance, from(init_specs); see [R] maximize. These
options are seldom used, but the difficult option may be useful if there
are convergence problems.

Remarks

Let t denote time. stmt works by first fitting a flexible parametric survival
model on the log cumulative hazard scale using the stpm2 command. The predicted hazard
from this model is used to estimate initial values of the log hazard with specified
covariates. These initial values are used within the maximum likelihood estimation.

Mutliple timescales are included in the model as a function of each other. For example, if one
wishes to model both the effect of time since diagnosis of a disease (timescale 1), and attained age
(timescale 2), timescale 2 is equal to timescale 1 plus the age at diagnosis. Standard models which
use maximum likelihood with analytical forms of the log likelihood function cannot utilize this relationship
between timescales. Here, numerical integration is used to get the log likelihood function meaning that
modelling including multiple timescales in this way is valid. Flexible parametric survival models on the
log hazard scale are used with Gaussian Quadrature.
The log hazard estimate, ln(h_hat(t)), is smoothed on ln(t), or t if logtoff is
specified, using restricted cubic splines with knots placed at certain quantiles of the
distribution of t. The knot positions are chosen automatically if the spline complexity
is specified by the df() option, or manually by way of the knots() option.

The survival function is calculated as

S_hat(t) = exp(-exp H_hat(t)).

With df(1) a Weibull model is fitted.

Estimation is performed by maximum likelihood. Optimisation uses the
default technique (nr, meaning Stata's version of Newton-Raphson
iteration.

Examples

Setup

webuse brcancer

Rename variable for age at diagnosis

rename x1 agediag

Modelling one timescale (time since diagnosis)

stset rectime, failure(censrec = 1) scale(365.24)

stmt hormon, time1(df(4))

Modelling one timescale: time-dependent effects on the log hazard scale

stmt hormon, time1(df(4) tvc(hormon) dftvc(3))

Modelling two timescales (time since diagnosis and attained age)

stmt hormon, time1(df(4)) time2(start(agediag) df(2))

Stored results

stmt stores the following in e():

Scalars

e(rank) rank of e(V)

e(N) number of observations

e(ic) number of iterations

e(k) number of parameters

e(k_eq) number of equations in e(b)

e(k_dv) number of dependent variables

e(converged) 1 if the model converged

e(rc) return code

e(ll) log likelihood

e(k_eq_model) number of equations in overall model test

e(dfbase_t#) degrees of freedom specified for timescale #

e(Ntimeint) number of timescale interactions #

e(df_timeint#_t1) degrees of freedom used for timescale 1 of timescale interaction #

e(df_timeint#_t2) degrees of freedom used for timescale 2 of timescale interaction #

e(df_timeint#_t3) degrees of freedom used for timescale 3 of timescale interaction #

e(minknot) value of minimum knot

e(maxknot) value of maximum knot

e(nodes) number of nodes used in numerical integration

e(dev) deviance

e(AIC) Akaike information criterion

e(BIC) Bayesian information criterion

e(Ntimescales) BNumber of timescales included in model

Macros

e(cmdline) command as typed

e(cmd) stmt

e(predict) stmt_pred

e(varlist) variables included in the model

e(depvar) _d _t

e(orthog) orthog if splines were orthogonalized

e(indicator_t#) orthog if timescale # only applied to a subset of observations defined by the indicator variable

e(rcsterms_base_t#) name of spline terms used to model timescale #

e(bhknots_t#) value of the knots used to model timescale #

e(exp_bhknots_t#) exponential of the knots used to model timescale #

e(tvc_t#) variables to be included in a timescale-covariate interaction #

e(rcsterms_t#_varname) name of spline terms used to model the timescale for the varname-timescale interaction #

e(tvcknots_t#_varname) value of the knots used for the interaction between varname and the timescale #

e(exp_tvcknots_t#_varname) exponential of the knots specified in e(tvcknots_t#_varname) #

e(knots_timeint#_t1) knots for the first timescale as part of timescale interaction #

e(knots_timeint#_t2) knots for the second timescale as part of timescale interaction #varname on timescale #

e(knots_timeint#_t3) knots for the third timescale as part of timescale interaction #varname on timescale #

e(noconstant) noconstant if constant term was supressed

e(opt) type of optimization

e(vce) vcetype specified in vce()

e(user) name of likelihood-evaluator program

e(ml_method) type of ml method

e(technique) maximization technique

e(which) max or min; whether optimizer is to perform maximization or minimization

e(properties) b V

Matrices

e(b) coefficient vector

e(V) variance-covariance matrix of estimators

e(Cns) constraints matrix

e(R_bh_t#) orthogonlization matrix for timescale #

e(ilog) iteration log

e(gradient) gradient vector

e(R_varname) orthogonalization matrix for splines of timescale-covariate interactions varname

e(R_timeint#_t1) orthogonalization matrix for splines of timescale 1 in timescale interaction #

e(R_timeint#_t2) orthogonalization matrix for splines of timescale 2 in timescale interaction #

e(R_timeint#_t3) orthogonalization matrix for splines of timescale 3 in timescale interaction #

Functions

e(sample) marks estimation sample

Author

Hannah Bower, Karolinska Institutet, Stockholm, Sweden. (hannah.bower@ki.se)

Therese M-L Andersson, Karolinska Institutet, Stockholm, Sweden.

Michael Crowther, Karolinska Institutet, Stockholm, Sweden.

Paul Lambert, University of Leicester, UK.

References

H. Bower, M.J. Crowther, P. C. Lambert. strcs: A command for fitting flexible parametric
survival models on the log-hazard scale. The Stata Journal, 2016;16:989-1012.

P. Royston, P.C. Lambert. Flexible parametric survival analysis in Stata:
Beyond the Cox model StataPress, 2011

Also see

Online: [ST] stmt_postestimation;

Online: [ST] strcs;