*! version 1.0.0 2012

Title

stmixed -- Multilevel mixed effects survival analysis

Syntax

stmixed [fe_equation] || re_equation [|| re_equation ...]
[, options]

where the syntax of fe_equation is

[varlist] [if] [in]

and the syntax of re_equation for random coefficients and intercept is

levelvar: [varlist]
[, re_options]

levelvar is a variable identifying the group structure for the random
effects at that level.

re_options

Model
noconstant suppress constant term from the random-effects equation

options

Model
distribution(exponential) exponential survival distribution

distribution(gompertz) Gompertz survival distribution

distribution(weibull) Weibull survival distribution

distribution(ggamma) generalised gamma AFT model

distribution(lognormal) log normal AFT model

distribution(loglogistic) log logistic AFT model

distribution(rp) Royston-Parmar survival model

distribution(rcs) restricted cubic splines on the log hazard scale

distribution(pwexponential) piecewise-exponential survival distribution

distribution(user) a user-defined survival model; see user functions

df(#) degrees of freedom for baseline hazard function with d(fpm) or d(rcs)

knots(numlist) knot locations for baseline hazard with d(fpm) or d(rcs)

tvc(varlist) varlist of time-dependent effects

dftvc(df_list) degrees of freedom for each time-dependent effect's spline function

knotstvc(numlist) knot locations for time-dependent effects' spline function

noconstant suppress constant term from the fixed-effects equation

bhazard(varname) expected mortality rate; invokes a relative survival model

covariance(stmixed##vartype:vartype_list) variance-covariance structure of the random effects at each level

Integration
intmethod(intmethod) integration method

intpoints(#) set the number of integration points

adaptopts(adaptopts) options for adaptive quadrature

Estimation
from(matname) specify starting values

restartvalues(sv_list) specify starting values for specific random effect variances

apstartvalues(#) specify the starting value for all ancilary parameters; see details

zeros specify all initial values set to 0; see details

maximize_options control the maximization process; seldom used

user functions

llfunction(fnc_name) Mata function name which returns the observation-level log likelihood contribution

loghfunction(fnc_name) Mata function name which returns the observation-level log hazard function

hfunction(fnc_name) Mata function name which returns the observation-level hazard function

chfunction(fnc_name) Mata function name which returns the observation-level cumulative hazard function

nap(#) number of ancillary parameters to estimate

See merlin user-defined functions.

vartype

independent one variance parameter per random effect,
all covariances zero; the default unless a factor variable is specified

exchangeable equal variances for random effects,
and one common pairwise covariance

identity equal variances for random effects, all
covariances zero; the default for factor variables

unstructured all variances and covariances distinctly
estimated

intmethod

mvaghermite mean-variance adaptive Gauss-Hermite quadrature;
the default

ghermite nonadaptive Gauss-Hermite quadrature

mcarlo Monte-Carlo integration using Halton sequences or anti-thetic sampling; see details

adaptopts

[no]log whether to display the iteration log
for each numerical integral calculation

iterate(#) number of iterations to update integration points; default 1001

Reporting

showmerlin display the merlin syntax used to fit the model

level(#) set confidence level; default is level(95)

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

Weights are not currently supported.

Factor variables are not currently supported.

Description

stmixed fits multilevel survival models including standard parametric distributions,
flexible spline-based approaches (such as Royston-Parmar and the log hazard equivalent)
and allows the user to specify their own hazard function. Simple or complex time-dependent effects can be included,
with the addition of expected mortality for a relative survival model, and allowing for left-truncation/delayed entry.
t-distributed random effects can also be used instead of the default Gaussian.

See stmixed postestimation for a variety of predictions and residuals that can be
calculated following a model fit.

stmixed is part of the merlin family.

Options

Model

noconstant suppresses the constant (intercept) term and may be specified for the fixed effects equation and for the random effects equations.

distribution(string) specifies the survival distribution.

distribution(exponential) fits an exponential survival model.

distribution(weibull) fits a Weibull survival model.

distribution(gompertz) fits a Gompertz survival model.

distribution(ggamma) fits a generalised gamma AFT survival model.

distribution(lognormal) fits a log normal AFT survival model.

distribution(loglogistic) fits a log logistic AFT survival model.

distribution(rp) fits a Royston-Parmar survival model. This is a highly flexible fully parametric
alternative to the Cox model, modelled on the log cumulative hazard scale using restricted cubic splines.

distribution(rcs) fits a log hazard scale flexible parametric survival model. This is a highly flexible
fully parametric alternative to the Cox model, modelled on the log hazard scale using restricted cubic splines.

distribution(pwexponential) fits a piecewise-exponential survival model; requires the knots() option.

distribution(user) specify a user-defined survival model; see options below and
merlin user-defined functions.

df(#) specifies the degrees of freedom for the restricted cubic spline function used for the baseline
function under a rp or rcs survival model. # must be between 1 and 10, but usually a value
between 1 and 5 is sufficient.
The knots() option
is not applicable if the df() option is specified. The knots are placed at the following centiles of
the distribution of the
uncensored log survival times:

df knots Centile positions

1 0 (no knots)
2 1 50
3 2 33 67
4 3 25 50 75
5 4 20 40 60 80
6 5 17 33 50 67 83
7 6 14 29 43 57 71 86
8 7 12.5 25 37.5 50 62.5 75 87.5
9 8 11.1 22.2 33.3 44.4 55.6 66.7 77.8 88.9
10 9 10 20 30 40 50 60 70 80 90

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.

knots(numlist) defines either:

the knot locations for the baseline distribution function under a
rp or rcs survival model, as opposed to the default
locations set by df(). Note that the locations of the knots are placed on the standard time scale.
However, the scale used by the restricted cubic spline
function is always log time. Default knot positions are determined by the df() option.

the knot locations (cut-points) for the piecewise-exponential distribution. knots() should be
specified in increasing order.

tvc(varlist) gives the name of the variables that have time-varying coefficients.
Time-dependent effects are fitted using restricted cubic splines. The degrees of freedom are
specified using the dftvc() option.

dftvc(df_list) gives the degrees of freedom for time-dependent effects
in df_list. The potential degrees of freedom are listed under the
df() option. With 1 degree of freedom a linear effect of log time is fitted.
If there is more than one time-dependent effect and different degress of freedom
are requested for each time-dependent effect then the following syntax applies:

knotstvc(knotslist) defines numlist knotslist as the location
of the interior knots for time-dependent effects.

bhazard(varname) specifies the variable which contains the expected mortality rate, which invokes a
relative survival model.

covariance(vartype_list), where each vartype is

diagonal|exchangeable|identity|unstructured

specifies the structure of the covariance
matrix for the random effects. An diagonal covariance structure allows a distinct
variance for each random effect within a random-effects equation and
assumes that all covariances are zero. exchangeable covariances
have common variances and one common pairwise covariance. identity
is short for "multiple of the identity"; that is, all variances are equal
and all covariances are zero. unstructured allows for
all variances and covariances to be distinct. If an equation consists of
p random-effects terms, the unstructured covariance matrix will have
p(p+1)/2 unique parameters.

covariance(diagonal) is the default.

Integration

intmethod(intmethod),
intpoints(#), and
adaptopts(adaptopts)
affect how integration for the latent variables is numerically
calculated.

intmethod(intmethod) specifies the method and defaults
to intmethod(mvaghermite). The current implementation uses mean-variance adaptive quadrature
at the highest level, and non-adaptive at lower levels. Sometimes it is useful to fall back
on the less computationally intensive and less accurate
intmethod(ghermite) and then perhaps use one of the other more accurate
methods.

intmethod(mcarlo) tells merlin to use Monte-Carlo integration, which either uses Halton
sequences with normally-distributed random effects, or anti-thetic random draws with t-distributed
random effects.

intpoints(#) specifies the number of integration points
to use and defaults to intpoints(7) with intmethod(mvaghermite) or intmethod(ghermite),
and intpoints(150) with intmethod(mcarlo). Increasing the number
increases accuracy but also increases computational time.
Computational time is roughly proportional to the number specified.

adaptopts(adaptopts) affects the adaptive part of
adaptive quadrature (another term for numerical integration) and
thus is relevant only for intmethod(mvaghermite).

adaptopts() defaults to
adaptopts(nolog iterate(1001)).

[no]log
specifies whether iteration logs are shown each
time a numerical integral is calculated.

iterate(#) specifies the number of iterations to update the
integration points, which will include updating prior to iteration 0 in
the maximisation process.

Estimation

from(matname) allows you to specify starting values.

restartvalues(sv_list) allows you to specify starting values for specific random effect variances. See
merlin estimation for further details.

apstartvalues(#) allows you to specify a starting value for all ancillary parameters, i.e those defined by
using the nap() option.

zeros tells merlin to use 0 for all parameters starting values, rather than fit the fixed effect model. Both restartvalues()
and apstartvalues() can be used with zeros.

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.

Reporting

showmerlin displays the merlin syntax used to fit the model.

level(#) specifies the confidence level, as a percentage, for confidence intervals. The default is level(95) or as set by set level.

Example 1

This is a simulated example dataset representing a multi-centre trial scenario, with 100 centres and each centre recuiting 60 patients, resulting in
6000 observations. Two covariates were collected, a binary covariate x1 (coded 0/1), and a continuous covariate, x2, within the range [0,1].

Load dataset:

. use http://fmwww.bc.edu/repec/bocode/s/stmixed_example1

stset the data:

. stset stime, f(event=1)

We fit a mixed effect survival model, with a random intercept and Weibull distribution, adjusting for fixed effects of x1 and x2.

. stmixed x1 x2 || centre: , dist(weibull)

Example 2

This is a simulated example dataset representing an individual patient data meta-analysis, with 15 trials and each trial recuiting 200 patients, resulting in
3000 observations. We are interested in the pooled treatment effect, accounting for heterogeneity between trials.

Load dataset:

. use http://fmwww.bc.edu/repec/bocode/s/stmixed_example2

stset the data:

. stset stime, f(event=1)

Create dummy variables for trial membership:

. tab trial, gen(trialvar)

We fit a flexible parametric model with 3 degress of freedom for the baseline, proportional trial effects with trial = 1 as the reference,
and a random treatment effect.

. stmixed treat trialvar2-trialvar15 || trial: treat, nocons dist(rp) df(3)

Author

Michael J. Crowther

Red Door Analytics AB

Stockholm, Sweden

E-mail: michael@reddooranalytics.se.

Please report any errors you may find.

References

Crowther MJ, Riley RD, Staessen JA, Wang J, Gueyffier F, Lambert PC. Individual patient data meta-analysis of survival data using Poisson regression models. BMC Med Res Methodol 2012;12:34.

Crowther MJ, Look MP, Riley RD. Multilevel mixed effects parametric survival models using adaptive Gauss-Hermite quadrature with application to recurrent events and IPD meta-analysis. Statistics in Medicine 2014;(In Press).

Crowther MJ. Extended multivariate generalised linear and non-linear mixed effect models. 2017; Under review.