*! version 1.0.0 2021-02-28

Title

standsurv -- Standardized survival and related functions

Syntax

standsurv
{stub* | newvar | newvarlist}
[if] [in]
[, [options]

at1()...atn() fix specific covariate values for each cause

atvars() the new variable names (or stub) for each atn() option

atreference() the reference atn() option (default 1)

centile(numlist) centiles of the standardised survival function

centileupper(#) upper starting value when calculating centiles

centvar() the new variable to denote centiles

ci calculates confidence intervals for each atn() option and for contrasts

cif calculates cumulative incidence function for competing risks models

contrast() perform contrast between covariate patterns defined by atn() options

contrastvars() the new variable names (or stub) for each contrast

crmodels(modellist) names of competing risk models

crudeprob calculate crude probabilities

failure calculate standardised failure function (1-S(t))

frame(framename) save predictions to a frame

nogen do not generate at() option variables

genind(stub) output individual predictions

hazard calculate hazard function of standardised survival curve

indweights(varname) variable containing weights (for external standardisation)

lincom(numlist) linear combination of at options

level(#) sets confidence level (default 95)

mestimation use M-estimation for standard errors & confidence intervals

nodes(#) number of nodes for numerical integration (default 30)

ode use ordinary differential equations for some integrations

odeoptions(options) options for ordinary differential equations

over(varlist) estimate effects at unique values of varlist

rmst calculate restricted mean survival time

rmft calculate restricted mean failure time

se calculates standard errors for each atn() option and for contrasts

storeg(name) stores derivatives for each at option

stub2(name) override the default stubnamesn() option for competing risks

timevar(varname) time variable used for predictions (default _t)

toffset(varname) time offset when multiple models

transform() transformation to calculate standard errors when obtaining confidence intervals

userfunction() user defined function

userfunctionvar() the new variable names (or stub) for each user defined function

Description

standsurv is a postestimation command that calculates various standardized
(marginal) measures after fitting a parametric survival model. These include
standardized survival functions and a variety of functions of standardized survival
functions and contrasts of these functions. These functions include
the failure and hazard functions restricted mean survival time and centiles
of the marginal failure function.

When standardizing, specific covariate(s) can be held constant and contrasts between
different groups can be made, for example differences and ratios. Confidence intervals for all quantities are available.
User-defined transformations can be calculated by providing a user-written
Mata function.

The command also allows more than one survival model to be specified in a
competing risks setting. This allows calculation of standardized cause-specific
cumulative incidence functions and other useful measures.

Survival model estimation commands supported include stpm2, stpm3, strcs, and streg.
Note generalized gamma models are not currently implemented for streg models.

Factor variables are only supported for stpm3 models.
For other models you will need to create dummy variables and interactions
yourself when fitting the models.

By default new variables are saved to the original dataset.
However, it is advised that you use the frame option to save results to a new frame.

standsurv creates the following variables by default. However, it is recommended
that you use the atvar(), contrastvar(), lincomvar() and userfunctionvar() options
to give more meaningful variable names.
Note that you can now specify new variables names before the comma, as an alternative to
using the atvars() option.

_ati - standarized estimate for the ith ati() option

If the ci option is specified then the following variables are created

_ati_lci - lower confidence interval for standarized estimate for the ith ati() option

_ati_uci - upper confidence interval for standarized estimate for the ith ati() option

If the contrast() option is specified then the following variables are created

_contrastj_i - contrast between the atj() and ati() options.
By default at1() is the reference, but this can be changed using the atref() option.

_contrastj_i_lci - lower confidence interval for contrast between the atj() and ati() options.

_contrastj_i_uci - upper confidence interval for contrast between the atj() and ati() options.

If the lincom() option is specified then the following variables are created

_lincom - linear combination of standardized estimates.

_lincom_lci - lower confidence interval of linear combination of standardized estimates.

_lincom_uci - upper confidence interval of linear combination of standardized estimates.

If the userfunctionvar() option is specified then the following variables are created

_userfunction - function of standardized estimates specified using userfunction() option.

_userfunction_lci - lower confidence interval for function of standardized estimates specified using userfunction() option.

_userfunction_uci - upper confidence interval for function of standardized estimates specified using userfunction() option.

Options

at1(varname # [varname # ..], suboptions)..atn(varname # [varname # ..], suboptions)
specifies covariates to fix at specific values when averaging predictions.
For example, if x denotes a binary covariate and you want to standardise
over all other variables in the model then using at1(x 0) at2(x 1) will give
two standardised functions, one where x=0 and one where x=1.

Using at1(.) will calculated the standardized quantity with all observations set to their observed values.

It can be sometimes be useful to set certain variables to take the values of a
different covariate. This can be done using at1(x1 = x2) for example.
This can be useful when there are interactions: consider a model with
treat age treat_age as covariates where treat_age is an interaction
between treat and age. When standardising for treat=0 and treat=1,
the at() options should be at1(treat 0 treat_age 0) and
at2(treat 1 treat_age = age).

If you use stpm3 then it is advised you use factor variables when using
interactions. This makes specification of the options much easier,
particularly when there are interactions.

There are some suboptions. These are,

atif(expression) restricts the standardization to that selected by the
expression. For example, if x is an exposure covariate taking 0 for the
unexposed and 1 for the exposed, then using at1(x 0, atif(x==1))
standardizes over the covariate distribution among the exposed.
Note atif() allows different if expressions for each at() option.
Often the same if expression is required for each at() option and so
a standard single if/in statement can be used.

atenter(#) specifies the start time when integrating for RMST or RMFT.

atindweights(varname) Multiplies each individual prediction by the specified
varname. This can be used for external (age) standardization. If the same weights
are being used for all at() options then the main indweights() option can
be used.

attimevar(varname) specifies a different time variables for each ati() option.
Here you should be careful with contrasts.

atvars(stub | newvarnames) gives the names of the new variables to be
created for each ati() option.
This can be specified as a varlist equal to the number of at()
options or a stub where new variables are named stub1 -
stubn. If this option is not specified, the names default to
_at1 to _atn.

atreference(#) the at#() option that defines the reference category.
By default this is at1().

centile(numlist) calculates centiles of the standardised survival curve
for the centiles given in numlist. The centile values are given in a
new variable, _centvar, or that defined using the centvar() option.

centileupper(#) upper starting value when calculating centiles of the
standardised failure curve. The default is four times the maximum survival time
if there is a value of _t, otherwise it is 100.
If you have to set this, and the analysis data is in memory,
it probably means your estimate is based on extrapolating
the survival function way beyond your observed follow-up.

centileupper(newvarname) name of new varaible giving values of centiles.
The default is _centvar.

ci calculates a level(#)% confidence interval for each standardised
function or contrast. The confidence limits are stored using the
suffix _lci and _uci.

cif calculates the cause-specific cumulative incidence functions from
competing risks models. This must be used with crmodels() option to list the
cause-specific models. See XXX for naming rules...

contrast(contrastname) calculates contrasts between standardised measures.
Options are difference and ratio. There will be n-1
new variables created, where n is the number of at() options.
See above for naming rules.

contrastvars(stub | newvarnames) gives the new variables to create when
using the contrast() option. This can be specified as a varlist or a stub,
whereby new variables are named stub1 - stubn-1.
The names default to _contrast1 to _contrastn-1.

crmodels(modellist) gives the names of the cause-specifc models for
competing risks. Each model needs to have been saved using estimates store.

crudeprob calculates crude probabilities. This assumes that a relative
survival model has been fitted and that expected mortality rates have been
given using the expsurv() option.

enter(#) gives the enter time for conditional estimates. This currently only
works with the rmst and rmft options.

expsurv(suboptions) indicates that expected survival should be calculated
and then combined with the model based (relative) survival estimates to give
all cause survival. There are a number of suboptions,

agediag(varname) gives the variable giving the age at diagnosis in years
for subjects in the study population. Note that when using integer age is
is assumed that individuals were diagnosed on their birthday, so avoid this if possible.

datediag(varname) gives the variable giving the date at diagnosis
for subjects in the study population. If you do not have exact dates then you
still need to specify this option. For example, if you had dates in months
and years you could use

. gen datediag = mdy(diagmonth,1,diagyear)

to assume all subjects were diagnosed on the 1st of the month.

genind(varlist) gives the names of new variables to create and store the
individual contributions to the marginal expected survival, failure, rmst or rmft.
Note this only works if the main genind() option is specified.

expsurvvars(stub | newvarnames) gives the names of the new variables to be
created for marginal expected survival (survival option) or marginal expected
life expectency (RMST option). If not specified these variables will
not be stored. The names can be specified as a varlist equal to the number
of at() options or a stub where new variables are named stub1 -
stubn.

pmage(varname) gives the name of the age variable in the population mortality file.

pmmaxage(varname) gives the maximum age in the population mortality file.
When calculating attained age to merge in the expected mortality rates any record
that is over this maximum will be set to this maximum. The default value is 99.

pmmaxyear(varname) gives the maximum year in the population mortality file.
This is potentially useful when extrapolating. When calculating attained year
to merge in the expected mortality rates any record that is over this maximum
will be set to this maximum.

pmother(other) gives the name of other variables in the population
mortality file. For example, this is usually sex, but can also be region,
deprivation etc.

pmrate(rate) gives the rate variables in the population mortality file. Note
that standsurv requires the expected mortality rate and not the expected survival.

pmyear(varname) gives the name of the calendar year variable in the population
mortality file.

split(varname) gives the split times when calculating expected survival (for rmst).
The default is one year.

using(filename) filename of population mortality file.

failure calculates the standardised failure function rather than
the standardised survival function.

frame(framename, options) specifies the name of the frame to store the predictions.
The default name is stpm2_pred. Options are,

merge merge predictions into frame framename. This will extract the
time variable from the framename automatically.

mergecreate is a convenience option when writing predictions in a loop.
If the frame does not exist it will be created.
If the frame does exists then the timevar() option will be ignored
and predictions merged into framename.

replace replace the existing frame framename is it exists.

nogen Does not generate variables for each at() option. This is only
relevent when using the contrast() option and you are not interested
in the predictions for each at option.

genind(stubname) outputs the predictions for each individual, i.e. the
predicted values that feed into the average. Standsurv is concerned with marginal
estimates, but it is sometimes of interest to look at the variation between
individuals. Note that this only will work if timevar is a single value. It
was implemented as a debugging tool when developing the command and may not work for all
options.

hazard calculates the hazard function of the standardised survival
function. Note that this is not the mean of the predicted hazard functions,
but a weighted mean with weights, S(t). The weights are time-dependent.

indweights(varname) Multiplies each individual prediction by the specified
varname. This is used for external (age) standardization.

level(#) sets the confidence level; default is level(95) or as set by set level.

lincom(#...#) calculates a linear combination of atn options. As an example, if
there were two at() options then lincom(1 -1) would calculate the difference in the
standardized estimate. This would be the same as using the contrast(difference) option.
Note that in the competing risks setting the number of values in the lincom() options
should be the number of atn options muliplied by the number of competing risk models.
In such as case the first value corresponds to at1 competing risk model 1 and their
second value to at1 competing risk model 2, i.e. the values correspond to the order in which
the new standardized variables are created. For example with two competing risk models and
two at() options combined with the cif option, lincom(1 1 0 0) would sum the CIFs
for cause 1 and 2 for at1.

lincomvar(newvarname) gives the new variable to create when
using the lincom() option.

mestimation requests that standard errors are obtained using M-estimation
(Stefanski and Boos 2002) rather than the delta-method.

nodes(#) number of nodes when performing numerical integration using
Gaussian quadrature to calculate the restricted mean survival time.

ode use ordinary differential equations (Dormand Prince 45) for numerical integration.
This the default and only option for various predictions and so only relevent
when the default is Gaussian quadrature.

odeoptions(options) Various options for ordinary differential equations.

abstol(#) absolute tolerance - default 1e-6.

error_control(#) error control when reducing step size - default (5/safety)^(1/pgrow).

initialh(#) initial step size - default 1e-8.

maxsteps(#) maxiumum number of steps - default 1000.

pgrow(#) power using when increasing step size.

pshrink(#) power using when decreasing step size.

reltol(#) relative tolerance - default 1e-5.

safety(#) safety factor - default 1.

tstart(#) lower bound of integration.

verbose output details for each step.

over(varlist) specifies that separate standardized estimates are obtained
for the groups defined by varlist. The variables need not be
covariates in your model.

per(#) multiplies predictions by #. For example, when predicting
survival per(100) will express results as a percentage rather than a
proportion, or when predicting hazard functions, per(1000) gives results
per 1000 person years (if your time scale is in years of course).

rmst calculates the restricted mean survival time. These are calculated at
the time points give in variable given in the timevar() option.

rmft calculates the restricted mean failure time. These are calculated at
the time points give in variable given in the timevar() option.
Note that if the rmft option is used together with the cif option
for competing risks models then the area under the CIF is calculated.
This is sometimes referred to as the expected loss in life due to a specific cause.
See Andersen (2013).

se calculates the standard error for each standardised
function or contrast. This is stored using the suffix _se.

storeg(name) store the derivatives of the standsardized estimate w.r.t each
of the model parameters. This can be useful if you need to combine results from
multiple standsurv calls.

stub2 Overide default stub names_se.

timevar(varname) defines the variable used as time in the predictions. The
option is useful for large datasets where, for plotting purposes, predictions
are needed only for (say) 200 observations. Note that predictions are averaged
over the whole sample, not just those where timevar is not missing. It is
recommended that timevar() is used, as otherwise an estimate of the survival
function is obtained at each value of _t for all subjects.
Default varname is _t.

toffset(varlist) defines variables that have an offset for time for use when
predicting cause-specific cumulative incidence functions. For example, if model 1
used time since diagnosis as the time-scale and model 2 uses attained age as
the time scale, then using toffset(. agediag) will make all predictions
on the time since diagnosis time scale, but making appropriate adjustments for
age at diagnosis.

trans(name) transformation to apply when calculating standard errors to
obtain confidence intervals for the standardised curves. The default transformation is
log(). Other possible names are none, loglog, logit.

userfunction(name) gives a Mata function that calculates transformations the
standardized functions. This enables flexibility to calculate
a wide range of potential functions. An example of a Mata function to calculate
a difference between two standardized function is shown below

mata:

function user_eg(at) {

return(at[2] - at[1])

}

end

standsurv, at1(x1 0) at2(x1 1) timevar(tt) ci userfunction(user_eg)

userfunctionvar(newvarname) gives the new variable to create when
using the userfunction() option. The name defaults to _userfunc.

verbose gives details about what is being estimated during the running
of the command, mainly used for developing and debugging.

Example:

For some more detailed examples see https://pclambert.net/software/standsurv/

Load example dataset:

. webuse brcancer

stset the data:

. stset rectime, f(censrec==1) scale(365.24)

Fit stpm2 model:

. stpm2 hormon x5 x1 x3 x6 x7, scale(hazard) df(4) tvc(hormon x5 x3) dftvc(3)

Generate variable that defines timepoints to predict at. The following creates 50 equally spaced time points between 0.05 and 5 years:

. range timevar 0 5 50

Obtain standardised curves for hormon=0 and hormon=1.
In each case the survival curves are the average of the 686
survival curves using the observed covariate values except for hormon.

. standsurv, atvars(S0a S1a) at1(hormon 0) at2(hormon 1) timevar(timevar) ci frame(f1,replace)

Plot standardised curves:

. frame f1: line S0a S1a timevar

Obtain standardised curves for hormon=0 and hormon=1, but apply the covariate distribution amongst those with hormon=1.

. standsurv if hormon==1, atvars(S0b S1b) at1(hormon 0) at2(hormon 1) timevar(timevar) ci frame(f2, replace)

Plot standardised curves:

. frame f2: line S0b S1b timevar

Obtain standardised curves for hormon=0 and hormon=1, and calculate difference in standardised survival curves and 95 confidence interval.

. standsurv, atvars(S0c S1c) at1(hormon 0) at2(hormon 1) timevar(timevar) ci contrast(difference) contrastvar(Sdiffc) frame(f3, replace)

Plot difference in standardised curves and 95% confidence interval:

. frame f3: line Sdiffc* timevar

Authors

Paul C. Lambert

Biostatistics Research Group

Department of Population Health Sciences

University of Leicester

and

Department of Medical Epidemiology and Biostatistics

Karolinska Institutet

E-mail: paul.lambert@le.ac.uk

Michael J. Crowther

Biostatistics Research Group

Department of Health Sciences

University of Leicester

Please report any errors you may find.

References

Andersen PK. Decomposition of number of life years lost according to causes of death.
Statistics in Medicine 2013;32;5278-5285

Crowther MJ, Lambert PC. Parametric multi-state survival models: flexible modelling allowing transition-specific distributions with
application to estimating clinically useful measures of effect differences. Statistics in Medicine 2017;36:4719-4742.

Stefanski L.A. and Boos, DD. The calculus of M-estimation. The American Statistician 2002;56;29-38.