*! version 1.0.0

Title

predictms -- predictions from a multi-state model

Syntax

predictms , transmatrix
statistic
[model_spec
prediction_spec
options]

transition matrix

transmatrix(matname) transition matrix

singleevent shorthand to define a transmatrix() for a single event survival model

cr shorthand to define a transmatrix() for a competing risks model

One of the above options must be specified

statistic

probability calculate transition probabilities

los calculate length of stay in each state

rmst calculate restricted mean survival time

visit probability of ever visiting each state

hazard transition-specific hazard function

survival transition-specific survival function

userfunction(func_name) user-defined Mata function for bespoke predictions; see details

At least one of the above options must be specified

model specification

models(namelist) list of estimates stored for all transition hazards

reset use clock-reset approach

tscale2(numlist) transition models on a second timescale

time2(numlist) time to add to main timescale

prediction specification

from(numlist) starting state(s) for predictions

ltruncated(#) starting time, i.e. time at which the starting state(s) from() are entered

timevar(varname) time points at which to calculate predictions

mint(#) minimum time at which to calculate predictions

maxt(#) maximum time at which to calculate predictions

obs(#) number of time points to calculate predictions at between mint() and maxt()

at#(at_list) calculate predictions (and contrasts) at covariate patterns

standardise calculates standardised/population-averaged predictions; see details

standif(condition) restricts the observations that are standardised over; see details

difference calculate differences between predictions

ratio calculate ratios of predictions

atreference(#) specifies the reference prediction for difference and ratio contrasts

userlink(string) link function used in calculation of confidence intervals for userfunction(); default identity

outsample for out of sample predictions; see below.

confidence intervals

ci calculate confidence intervals of predictions

level(#) calculate confidence intervals at specific level, default is 95

bootstrap calculate confidence intervals using the parametric bootstrap

m(#) number of bootstrap repetitions for calculating confidence intervals

percentile calculate confidence intervals using percentiles of the bootstrap repetitions

novcv(numlist) transition models that should be assumed are not estimated with uncertainty; see details

options

aj use the Aalen-Johansen estimator for transition probabilities; see details

simulate calculate predictions using large-sample simulation; see details

latent use latent times as the method of simulation; see details

n(#) sample size of simulated dataset

chintpoints(#) number of Gauss-Legendre quadrature points; see details

save(name, [replace]) save each simulated dataset used to calculate predictions and confidence intervals; see details

seed(#) set the simulation seed

Description

predictms calculates a variety of predictions from a multi-state survival model, including
transition probabilities, length of stay (restricted mean survival times in each state), the probability of ever visiting
each state, and more. Predictions are made at user-specified covariate patterns. Differences and ratios of predictions across
covariate patterns can also be calculated. Standardised (study population-averaged) predictions can be obtained. Confidence
intervals for all quantities are available. Numerical integration, simulation or the Aalen-Johansen estimator (with
parametric estimates of the cumulative hazard) are used to calculate all quantities. User-defined predictions can also be
calculated by providing a user-written Mata function, to provide complete flexibility.

The transition models must be estimated using the stmerlin or merlin commands, with the current range of
supported models including:

- exponential proportional hazards model

- Weibull proportional hazards model

- Gompertz proportional hazards model

- generalised gamma accelerated failure time model

- log normal accelerated failure time model

- log logistic accelerated failure time model

- piecewise-exponential proportional hazards model

- Royston-Parmar flexible parametric model

- general log hazard scale model

- general log cumulative hazard scale model

- general additive hazard model

- user-defined/custom hazard model

The user (usually) must provide the transition matrix used in the fitted model, through the transmatrix() option.
The transition matrix may be cyclic or acyclic. Default predictions assume all subjects start in state from(1), at time
ltruncated(0).

predictms creates the following variables:

_time times at which all predictions are calculated

If probability is requested, then predictms creates the following variables:

_prob_ati_a_b transition probability for ati() (from state a to state b)

_prob_ati_a_b_lci lower confidence interval of transition probability for ati() (from state a to state b)

_prob_at#_a_b_uci upper confidence interval of transition probability for ati() (from state a to state b)

If los is requested, then predictms also creates the following variables:

_los_ati_a_b length of stay in state b for ati() (given they started from state a)

_los_ati_a_b_lci lower confidence interval of the length of stay in state b for ati() (given they started from state a)

_los_ati_a_b_uci upper confidence interval of the length of stay in state b for ati() (given they started from state a)

If rmst is requested, then predictms also creates the following variables:

_rmst_ati_a total time spent in non-absorbing states for ati() (given they started from state a)

_rmst_ati_a_lci lower confidence interval of the total time spent in non-absorbing states for ati() (given they started from state a)

_rmst_ati_a_uci upper confidence interval of the total time spent in non-absorbing states for ati() (given they started from state a)

If visit is requested, then predictms also creates the following variables:

_visit_ati_a_b probability of ever visiting state b for ati() (given they started from state a)

_visit_ati_a_b_lci lower confidence interval of the probability of ever visiting state b for ati() (given they started from state a)

_visit_ati_a_b_uci upper confidence interval of the probability of ever visiting state b for ati() (given they started from state a)

If hazard is requested, then predictms also creates the following variables:

_hazard_ati_a_b predicted hazard function for the transition going from state a to b for ati()

_hazard_ati_a_b_lci lower confidence interval of the predicted hazard function for the transition going from state a to b for ati()

_hazard_ati_a_b_uci upper confidence interval of the predicted hazard function for the transition going from state a to b for ati()

If a userfunction() is provided, then predictms also creates the following variables:

_user_ati_a_b returned values for column b of the user-function for ati() (given they started from state a)

_user_ati_a_b_lci lower confidence interval of the returned values for column b of the user-function for ati() (given they started from state a)

_user_ati_a_b_uci upper confidence interval of the returned values for column b of the user-function for ati() (given they started from state a)

If a difference is requested, then predictms also creates some or all of the following variables:

_diff_prob_ati_a_b difference in transition probabilities between ati() and that specified in the reference atj()

_diff_los_ati_a_b difference in length of stay between ati() and that specified in the reference atj()

_diff_user_ati_a_b difference in the user-function between ati() and that specified in the reference atj()

*_lci and *_uci are returned as appropriate.

If a ratio is requested, then predictms also creates some or all of the following variables:

_ratio_prob_ati_a_b ratio in transition probabilities between ati() and that specified in the reference atj()

_ratio_los_ati_a_b ratio in length of stay between ati() and that specified in the reference atj()

_ratio_user_ati_a_b	ratio in the user-function between ati() and that specified in the reference atj()

*_lci and *_uci are returned as appropriate.

All appropriate combinations of a and b are returned. predictms first drops any variables it may have calculated
from a previous call (e.g. _prob_*).

predictms is part of the multistate package. Further
details here: mjcrowther.co.uk/software/multistate

Options

Transition matrix

transmatrix(matname) specifies the transition matrix which governs the multi-state model
that was fitted. 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.

singleevent indicates that you have fitted a standard single event survival analysis model. This corresponds to a
transition matrix of (.,1\.,.). transmatrix() does not need to be specified.

cr indicates that you have fitted a competing risks model, and is a useful way of avoiding having to specify a
transmatrix(). For use with models(), as the number of competing risks corresponds to the number of model
objects. For example, for a model with two competing events:

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

Statistic

probability calculate transition probabilities. If aj has been used then I recommend at least obs(500)
is used in your timevar().

los calculate (restricted) length of stay in each state, i.e. restricted mean survival time for each transient state.
This is the integral of the transition probabilities across follow-up time. If aj has been used then length of stay
is calculated using numerical integration of the transition probabilities, and I recommend at least obs(500)
is used in your timevar().

rmst calculate (restricted) mean survival time, i.e. the total time spent in non-absorbing/transient states.
If aj has been used then restricted mean survival time is calculated using numerical integration of the
transition probabilities, and I recommend at least obs(500) is used in your timevar().

visit calculate the probability of ever visiting each state within the time interval defined from mint() and
the prediction time. Can currently only be calculated using the large-sample simulation method.

hazard calculate the predicted hazard function for each transition from the from() state(s).

survival calculate the predicted survival function for each transition from the from() state(s).

userfunction(func_name) defines a Mata function which returns a real matrix, to calculate user-defined quantities.
The function must be of the form:

real matrix func_name(S)

{

...

some code

...

pred = some more code

return(pred)

}

where S is a transmorphic object which should not be changed. It is passed to utility functions, which give access to the
transition probabilites and/or length of stays, for example:

real matrix func_name(S)

{

p1 = ms_user_prob(S,1)

p2 = ms_user_prob(S,2)

p3 = ms_user_prob(S,3)

pred = p1,p2,p3

return(pred)

}

which would give us our transition probabilities for each state. Length of stays are accessed using ms_user_los() in the
same way, but the los option must also be specified.

Model specification

models(namelist) specifies the names of the estimates store objects containing the estimates of the model
fitted for transition 1, 2, 3, ..., for example,

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

. predictms , transmatrix(tmat) models(m1 m2 m3)

If models() is not specified, then predictms assumes you have fitted a single, stacked, multi-state model to
your dataset.

reset tells predictms that the clock-reset approach, i.e. a Markov renewal multi-state model has been fitted. The
default assumes clock-forward. When reset is specified, then the time at which predictions are made, stored in
_time, represent the main time scale, time since entry to the from() state(s), which may be left truncated if
ltruncated() is greater than 0.

tscale2(numlist) specifies any transition models (using the same index as specified in transmatrix()), which
are modelled on a secondary timescale (different to the main one), enabling transition-specific timescales (Weibull et al.
Under review.). This can be used for example when modelling age as the timescale for some transitions, and time since
diagnosis for others.

time2(numlist) specifies the value to be added to the main timescale at entry, for each of the transition models
specified in tscale2(). If age is the second timescale, then this would be age at baseline. Each element of
time2() corresponds to each at#(). If only one time2() is specified then it will be assumed for all
at#()s.

Prediction specification

from(numlist) define the starting state for all predictions (state a in the variable descriptions above), default
is state 1. Multiple starting states can be defined, which will calculate all possible predictions for each starting state.

ltruncated(#) defines the starting time, i.e. the time at which the starting state(s) from() are entered, with a default
of ltruncated(0)

timevar(varname) variable which contains time points at which to calculate predictions, which overrides the default
created by using mint(), maxt() and obs(). Predictions are made at these timepoints, conditional on
starting in state(s) from() at time ltruncated().

mint(#) minimum time at which to calculate predictions, with a default of mint(0). If timevar() is not
specified, then a default time variable will be created, using mint(), maxt() and obs().

maxt(#) maximum time at which to calculate predictions. If timevar() is not specified, then a default time
variable will be created, using mint(), maxt() and obs().

obs(#) the number of time points at which to calculate predictions, equally spaced between mint() and maxt(),
with a default of 100 (500 when using aj).

at#(vn # ...) calculates predictions at specified covariate patterns, e.g. at1(female 1 age 55).
Specifying multiple at#()s means only one call of predictms has to be made to calculate many predictions.
Note that any covariates not specified in at#() are set to 0.

standardise calculates standardised/population-averaged predictions. For each statistic requested, the prediction
will be averaged over all covariate patterns found either in your dataset, or in the observations identified using the
standif() condition. For each observation, any covariates not specified in at#() will take their observed values -
this is done for all observations, and the prediction is then averaged. Predictions can then be interpreted as marginal
predictions with respect to the observed covariate distributions that have been averaged over. A standardised prediction
can be very computationally intenstive to obtain.

standif(condition) restricts the observations that are standardised over, for example standif(_n<100), or
standif(trt==1)

difference calculate the difference between predictions, across different at#() specifications.
See atreference().

ratio calculate the ratio of the predictions, across different at#() specifications.

atreference(#) specifies the reference at#() for calculating prediction contrasts. Default is atref(1), meaning
that each prediction is contrasted to that calculated using at1().

userlink(link_name) specifies the link function used when calculating confidence intervals with the normal approximation
applied to the userfunction(). Options include log, logit, and the default identity. The link function
is applied, and the mean and standard error are then calculated on the transformed scale to calculate confidence intervals,
and then transformed back. This ensures all predictions are within the desired range (e.g. 0 and 1 for probabilities).

outsample specifies predictions are being made out of sample, which suppresses checks that variables specified in
at#()s are in your current dataset. This is of use when transition models() come from different datasets. As
this suppresses some error checks, users should be extra careful when specifying their at#()s.

ci requests confidence intervals of all predictions. The delta method or parametric bootstrap are used, depending on the model
specification and the requested statistics. The most efficient method will be chosen by default. The number of bootstrap draws from
the estimated coefficient vectors, and associated variance-covariance matrix, is specified in the m() option. When bootstrapping,
the default method uses a normal approximation of simulated repetitions to calculate the standard error, see also percentile.

level(#) confidence interval level, default level(95).

bootstrap forces predictms to use the parametric bootstrap to calculate confidence intervals, rather than the
delta method, which is the default for survival, competing risks & illness-death models.

m(#) number of parametric bootstrap repetitions for calculating confidence intervals; default is m(200). I
recommend using an increasing number to check a sufficiently large m() has been used.

percentile calculate confidence intervals based on centiles of the predictions across the m(#) sets, instead of
the default normal based calculation.

novcv(numlist) specifies transitions (using the same index as specified in transmatrix()) that when calculating
confidence intervals on predictions (using ci), are assumed to be estimated free of uncertainty, i.e. draws will not be
made from the multivariate normal centred on the parameter estimates and with VCV the estimated VCV from that transition, but
will use the estimated parameter vector every time.

Options

aj specifies a hybrid Aalen-Johansen estimator be used to calculate predictions, instead of either numerical integration or
large-sample simulation. This method uses the non-parametric formula for the AJ estimator, but uses parametric estimates of the
transition-specific cumulatve hazard functions from the fitted models(), or stacked model, as the ingredients, rather
than the non-parametric Nelson-Aalen estimates. This approach assumes constant hazards between timepoints, and hence more
timepoints should be used when calculating predictions, i.e. at least obs(100). This option is only valid with
Markov models so can't currently be used with a reset model. The aj method is substantially faster than the
large-sample simulation approach.

simulate forces predictms to use large-sample simulation to calculate predictions, rather than numerical integration
which is the default for survival, competing risks, and illness-death settings. Any other transition matrix structure will use
simulation by default.

latent specifies that latent times are used as the method of simulation, rather than the default method of Beyersmann et al. (2009).
The latent times approach conducts a series of competing risks simulations, simulating an event time from each cause-specific hazard
function, and taking the minimum as the observed event, compared to Beyersmann's method which simulated one event time from the
total hazard, and a multinomial draw to determine the event that occured, based on the cause-specific hazard contributions. Both
methods will give essentially identical results, but the Beyersmann method is generally more computationally efficient.

n(#) defines the number of individual simulated trajectories through the multi-state model which are used to calculate predictions
under the large-sample simulation approach. The default is n(100,000) unless ci is specified, then it is n(10,000).
Accuracy increases with a higher sample, reducing Monte Carlo error, but increasing computation time.

chintpoints(#) defines the number of Gauss-Legendre quadrature points used to calculate any analytically intractible integrals.
Numerical integration is required in the default prediction method for survival, competing risks and illness-death settings, and for
simulating transition times from spline and complex time-dependent models. The defaukt is chintpoints(30). It is good
practice to increase this number to ensure the predictions are reliable. Accuracy increases with a higher number, but increasing
computation time.

save(name, [replace]) saves the simulated trajectories dataset(s) generated when using large-sample simulation to calculate predictions and
confidence intervals. The dataset used to calculate point estimates will be saved in name.dta, and each simulated dataset generated in the
mth bootstrap sample will be saved in namei.dta. If the datasets already exist, you will need to add the replace
option to overwrite.

seed(#) sets the simulation seed. When using large-sample simulation to obtain predictions, it is important to set the seed
to ensure reproducibility.

Examples

Example 1: A separate transition-specific illness-death model

This dataset contains information on 2982 patients with breast cancer. Baseline is defined as time of surgery, and patients
can experience relapse, relapse then death, or death with no relapse. Time of relapse is stored in rf, with event
indicator rfi, and time of death is stored in os, with event indicator osi.

Load example dataset:

. use http://fmwww.bc.edu/repec/bocode/m/multistate_example

msset the data (from the multistate package):

. msset, id(pid) states(rfi osi) times(rf os)

Store the transition matrix:

. mat tmat = r(transmatrix)

stset the data using the variables created by msset

. stset _stop, enter(_start) failure(_status=1)

We fit separate Weibull models, so a fully stratified model, also allowing transition specific age effects:

. stmerlin age if _trans1==1, distribution(weibull)
. estimate store m1

. stmerlin age if _trans2==1, distribution(weibull)
. estimate store m2

. stmerlin age if _trans3==1, distribution(weibull)
. estimates store m3

Calculate transition probabilities for a patient with age 50:

. predictms, transmatrix(tmat) models(m1 m2 m3) probability at1(age 50)

Create a stacked transition probability plot:

. graphms

Calculate transition probabilities and length of stay for a patient with age 50:

. predictms, transmatrix(tmat) models(m1 m2 m3) at1(age 50) probability los

Calculate the difference in transition probabilities for a patient with age 60 compared to a patient aged 50:

. predictms, transmatrix(tmat) models(m1 m2 m3) probability at1(age 50) at2(age 60) difference

Calculate the ratio of transition probabilities for a patient with age 60 compared to a patient aged 50:

. predictms, transmatrix(tmat) models(m1 m2 m3) probability at1(age 50) at2(age 60) ratio

Calculate differences and ratios of length of stay and transition probabilities for a patient with age 60 compared to
a patient aged 50, with confidence intervals:

. predictms, transmatrix(tmat) models(m1 m2 m3) probability at1(age 50) at2(age 60) los diff ratio ci

Example 2: A stacked, proportional, illness-death model

Here we replicate the above example, but this time using a single merlin model, and using interactions between covariates
and transition-specific indicators to allow for transition-specific effects.

Load example dataset:

. use http://fmwww.bc.edu/repec/bocode/m/multistate_example

msset the data (from the multistate package):

. msset, id(pid) states(rfi osi) times(rf os)

Store the transition matrix:

. mat tmat = r(transmatrix)

stset the data using the variables created by msset

. stset _stop, enter(_start) failure(_status=1)

We fit proportional Weibull models, with a common age effect:

. merlin (_t age _trans2 _trans3 , family(weibull, failure(_d) ltruncated(_t0)))

We fit proportional Weibull models, with a transition-specific age effects:

. merlin (_t age#_trans1 age#_trans2 age#_trans3 _trans2 _trans3 , family(weibull, failure(_d) ltruncated(_t0)))

All predictions are as above, but we no longer need the models() option. Calculate transition probabilities
for a patient with age 50:

. predictms, transmatrix(tmat) probability at1(age 50)

Create a stacked transition probability plot:

. graphms

Calculate transition probabilities and length of stay for a patient with age 50:

. predictms, transmatrix(tmat) at1(age 50) probability los

Calculate the difference in transition probabilities for a patient with age 60 compared to a patient aged 50:

. predictms, transmatrix(tmat) probability at1(age 50) at2(age 60) difference

Calculate the ratio of transition probabilities for a patient with age 60 compared to a patient aged 50:

. predictms, transmatrix(tmat) probability at1(age 50) at2(age 60) ratio

Calculate differences and ratios of length of stay and transition probabilities for a patient with age 60 compared
to a patient aged 50, with confidence intervals:

. predictms, transmatrix(tmat) probability at1(age 50) at2(age 60) los diff ratio ci

Authors

Michael J. Crowther (1,*), Paul C. Lambert (2,3)

(1) Red Door Analytics, Stockholm, Sweden

(2) Department of Medical Epidemiology and Biostatistics, Karolinska Institutet, Sweden

(3) Biostatistics Research Group, Department of Health Sciences, University of Leicester, UK

(*) michael@reddooranalytics.se

Please report any errors you may find.

References

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

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(29):4719-4742.

Crowther MJ. Extended multivariate generalised linear and non-linear mixed effects models.
https://arxiv.org/abs/1710.02223

Crowther MJ. merlin - a unified framework for data analysis and methods development in Stata.
https://arxiv.org/abs/1806.01615

de Wreede LC, Fiocco M, Putter H. mstate: An R Package for the Analysis of Competing Risks and Multi-State Models.
Journal of Statistical Software 2011;38:1-30.

Putter H, Fiocco M, Geskus RB. Tutorial in biostatistics: competing risks and multi-state models.
Statistics in Medicine 2007;26:2389-2430.

Weibull CE, Lambert PC, Eloranta S, Andersson TM-L, Dickman PW, Crowther MJ. A multi-state model incorporating
estimation of excess hazards and multiple time scales. Statistics in Medicine 2021; (In Press).