*! version 0.1.0 ?????2017

Title

merlin -- Command syntax for model specification

Syntax

merlin (model1) [(model2)] [...] [, covariance()]

where the syntax of a model is

[depvar] [component1] [component2] [...] [if] [in] [, model_options]

where the syntax of a component is

element1[#element2][#element3][...][@real]

and each elementn can take one of the forms described in element,
with @real described in elements details.

model_options

family(family) distributional family

timevar(varname) time variable

noconstant omit the constant term

nap(#) number of ancillary parameters to estimate

family

exponential [, survival] exponential

bernoulli Bernoulli

beta beta

cox [, survival] Cox model

gamma gamma distribution

gaussian Gaussian (normal)

ggamma [, survival] generalised gamma

gompertz [, survival] Gompertz

hazard [, survival] additive hazard model

lquantile [, quantile(#)] linear quantile model with asymmetric Laplace distribution

logchazard [, survival] general log cumulative hazard model

loghazard [, survival] general log hazard model

loglogistic [, survival] log logistic

lognormal [, survival] log normal

nbinomial negative binomial with mean dispersion

null does not contribute to the log-likelihood; see details

ologit ordinal response with logit link

oprobit ordinal response with probit link

poisson Poisson

pwexponential [, survival pweopts] piecewise-exponential

rp [, survival rpopts] Royston-Parmar model on the log cumulative hazard scale

user [, user [survival]] user-defined

weibull [, survival] Weibull

survival

failure(varname) indicator for failure event

linterval(varname) lower interval time for interval censored observations

ltruncated(varname) entry time for left-truncated/delayed-entry model

bhazard(varname) expected mortality rate at event times, invokes a relative survival model

pweopts - family(pwexponential)

knots(knots_list) knot locations (cut-points), for the piecewise-exponential distribution. Should be in increasing order.

rpopts - family(rp)

df(#) degrees of freedom for the baseline log cumulative hazard function (does not include the intercept)

knots(knots_list) knot locations for the baseline log cumulative hazard function - includes
boundary knots, should be in increasing order.

noorthog turns off the default orthogonalisation of the spline terms

user

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

See merlin user-defined functions.

element

varname a variable in the dataset

M#[levelvar...] a random effect at the specified level; see details

fp(varname, fp_opts) a fractional polynomial function of varname

pc(varname, pc_opts) a piecewise-constant spline function of varname

rcs(varname, rcs_opts) a restricted cubic spline function of varname

bs(varname, bs_opts) a B-spline function of varname

mf(function_name) a user-defined Mata function; see details

EV[depvar|#] expected value of another outcome; see details

dEV[depvar|#] d/dt of expected value of another outcome; see details

d2EV[depvar|#] d2/dt2 of expected value of another outcome; see details

iEV[depvar|#] integral w.r.t time of the expected value of another outcome; see details

XB[depvar|#] expected value of the complex predictor of another outcome; see details

dXB[depvar|#] d/dt of the expected value of the complex predictor of another outcome; see details

d2XB[depvar|#] d2/dt2 of the expected value of the complex predictor of another outcome; see details

iXB[depvar|#] integral w.r.t time of the expected value of the complex predictor of another outcome; see details

More detailed descriptions can be found in elements details

fp_opts

powers(#) powers of the fractional polynomial; up to 2 degress can be used

offset(varname) to add before the fractional polynomial is calculated

moffset(varname) adds the negative of varname before the fractional polynomial is calculated

pc_opts

knots(knots_list) knot locations (cut-points) for the piecewise-constant spline function

noreference include an indicator variable for the first interval

offset(varname) to add before the spline function is calculated

moffset(varname) adds the negative of varname before the spline function is calculated

rcs_opts

df(#) degrees of freedom for the spline function; see details

knots(knots_list) knot locations for the spline function, including boundary knots

orthog use Gram-Schmidt orthogonalisation on the spline variables

event when using df(), calculate internal knot locations based on centiles of the observations that had an event (i.e. for survival models)

log calculate splines of the log of varname, rather than varname

offset(varname) to add before the spline function is calculated

moffset(varname) adds the negative of varname before the spline function is calculated

bs_opts

order(#) order of the B-spline function; see details

degree(#) degree of the B-spline function; see details

df(#) degrees of freedom for the spline function; see details

knots(knots_list) internal knot locations for the spline function

bknots(knots_list) boundary knot locations for the spline function

event when using df(), calculate knot locations based on centiles of the observations that had an event (i.e. for survival models)

log calculate splines of the log of varname, rather than varname

intercept include the intercept basis function

offset(varname) to add before the spline function is calculated

moffset(varname) adds the negative of varname before the spline function is calculated

Description

The command syntax of merlin is fully specified by the models and the covariance() option.

For full details and many tutorials, take a look at the accompanying website:
reddooranalytics.se/products/merlin

Options

Model

family(family) distributional family

timevar(varname) specifies the variable which represents time, and which will be used in conjunction with
elements. If you are modelling time-dependent effects in a survival analysis, for example, then timevar() should
match the response variable.

noconstant suppresses the constant (intercept) term in the extended linear predictor.

nap(#) specifies any ancillary parameters to be estimated. To be used in conjunction with a family(user) model,
or an mf() element, to be able the user to estimate any additional parameters required.

Elements

Elements are the fundamental flexibility of merlin. Each component can have any number of elements.
Elements within the same component are split using a #. Each element can be one of the following:

varname a variable in the dataset. Any of Stata's . notation is not currently supported.

M#[levelvar...] a random effect at the specified level.

M# defines a random effect, which must begin with M, and be followed by a positive integer

levelvar defines the cluster variables, for example, M1[level1] defines a random effect called M1
at level1, and M2[level1>level2] defines a random effect called M2 at level2 which is nested within level1.

fp(varname, fp_opts) specifies a fractional polynomial function of varname. fp_opts includes:

powers(numlist) can be any within the following
{-2, -1, -0.5, 0, 0.5, 1, 2, 3}. Fractional polynomials of degree 1 and 2 are allowed.

offset(varname) defines an offset to be added to the fp() variable prior to the fractional polynomial being derived.

moffset(varname) defines a negative offset to be taken away from the fp() variable prior to the fractional polynomial being derived.

pc(varname, pc_opts) a piecewise-constant spline function of varname, where pc_opts are:

knots(numlist) specifies the knot locations (cut-points), which must be in ascending order. Binary indicator variables are
created for each interval, except the first which is the reference group, and included in the linear predictor.

noreference specifies an indicator variable is also included for the first interval.

offset(varname) defines an offset to be added to the pc() variable prior to the spline function being derived.

moffset(varname) defines a negative offset to be taken away from the pc() variable prior to the spline function being derived.

rcs(varname, rcs_opts) a restricted cubic spline function of varname, where rcs_opts are:

df(#) specifies the degrees of freedom for the spline function. Boundary knots are assumed to be the min and max of
the varname. Internal knots are placed at equally spaced centiles.

knots(numlist) specifies the knot locations, which includes the boundary knots. Must be in ascending order.

log use splines of log varname rather than varname, the default.

orthog orthogonalise the spline terms.

event can be used in combination with df() and specifies that the internal knots are calculated
based on centiles of event times. Only valid when the family() is a survival model.

offset(varname) defines an offset to be added to the rcs() variable prior to the spline function being derived.

moffset(varname) defines a negative offset to be taken away from the rcs() variable prior to the spline function being derived.

bs(varname, bs_opts) a B-spline function of varname, where bs_opts are:

order(#) specifies the order of the spline function.

degree(#) specifies the degree of the spline function.

df(#) specifies the degrees of freedom (not strictly speaking) for the spline function, which allows you to specify
internal knots at equally spaced centiles, instead of using knots. df() is consistent with rcs() elements in how internal knots
are chosen.

knots(numlist) specifies the internal knot locations. Must be in ascending order.

bknots(numlist) specifies the lower and upper boundary knot locations.
Must be in ascending order. Default is the minimum and maximum of varname.

log use splines of log varname rather than varname, the default.

intercept includes the intercept basis function, which by default is not included.

event can be used in combination with df() and specifies that the internal knots are calculated
based on centiles of event times. Only valid when the family() is a survival model.

offset(varname) defines an offset to be added to the bs() variable prior to the spline function being derived.

moffset(varname) defines a negative offset to be taken away from the bs() variable prior to the spline function being derived.

mf(function_name) a user-defined Mata function. It should be defined as follows:

real matrix function_name(gml , real colvector t)

{

some code

}

Any of the utility functions can be used, which gives the user complete flexibility.
See merlin user-defined functions for further details on writing a Mata function
compatible with merlin.

EV[depvar|#] expected value of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar().

dEV[depvar|#] d/dt of expected value of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar(). The derivative is calculated using numerical differentiation.

d2EV[depvar|#] d2/dt2 of expected value of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar(). The second derivative is calculated using numerical differentiation.

iEV[depvar|#] integral w.r.t time of the expected value of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar(). The integral is calculated using numerical integration.

XB[depvar|#] expected value of the main complex predictor of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar().

dXB[depvar|#] d/dt of expected value of the main complex predictor of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar(). The derivative is calculated using numerical differentiation.

d2XB[depvar|#] d2/dt2 of expected value of the main complex predictor of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar(). The second derivative is calculated using numerical differentiation.

iXB[depvar|#] integral w.r.t time of the expected value of the main complex predictor of a model (other than the current one). Use the appropriate depvar to refer to
the desired model, or the model number # can be used (which may be needed if the model being referenced has family(null).
May also require timevar(). The integral is calculated using numerical integration.

@real specifies a constraint to be applied to all elements of the component. This is generally used in combination with random
effects, i.e to specify a random intercept, one would type M1[id]@1, which would constrain the coefficient of the random effect
M1 to be 1. Without the @1, merlin would attempt to estimate a coefficient.

Familys

family(family) specifies the distributional family. The inbuilt distributions listed above are self explanatory.
Some special cases include

family(lquantile [, quantile(#)]) specifies a linear quantile model utilising the asymmetric Laplace distribution for the particular quantile(#),
where # must be between 0 and 1. The default is quantile(0.5). In my experience, these models are very difficult to fit reliably.
I recommend you explore different integration techniques through intmethod(), along with increasing intpoints(),
and most importantly vary your starting values with from().

family(cox) fits a Cox model using maximum partial likelihood, assuming the Breslow method for handling ties.

family(user, user) which can take the following forms

family(user, llfunction(fnc_name)) defines a Mata function which returns the observation level log
likelihood contribution of a user-defined distributional model.

family(user, hfunction(fnc_name)) specifies a user-defined survival model, by defining a Mata
function which returns the observation-level hazard function. The cumulative hazard is calculated using numerical
quadrature, which means our cumulative hazard need not have a closed-form solution, such as when using splines of fractional
polynomials to model the baseline log hazard function.

family(user, hfunction(fnc_name) chfunction(fnc_name)) specifies a user-defined
survival model, by defining Mata functions which return the observation-level hazard and cumulative hazard functions.

See merlin user-defined functions for further details on writing a Mata function compatible with merlin.

family(null) tells merlin that theres is no depvar in the extended_lin_pred, and that
this model does not contribute to the log likelihood. It is a convenient way of specifying a further extended linear
predictor for use within one (or more) distributional model(s).

Survival

failure(varname) specifies the censoring/failure event. Should be coded 0 for right-censored observations,
1 for exactly observed events, or 2 for interval-censored observations.

linterval(varname) specifies the lower interval for interval censored observations. The upper interval
should be specified in the response variable for the model. Observations that are right censored or events should
be coded as missing.

ltruncated(varname) specifies the time at which observations become at risk of the event. Allows
fitting a delayed-entry survival model. If there are random effects in the associated survival model, then
the likelihood is calculated by dividing through by the marginal survival function at the entry times, which results
in a second set of numerical integration.

bhazard(varname) invokes a relative survival (excess hazard) model, by specifying the expected mortality (event) rate in the reference population at the observed event times.

Piecewise-exponential options

knots(knots_list) defines the knot locations (cut-points) for the piecewise-constant baseline function.
Knots should be specified in increasing order.

Royston-Parmar options

df(#) degrees of freedom for the baseline log cumulative hazard function, i.e. number of restricted cubic
spline terms. Internal knots are placed at centiles of the event times. Boundary knots are placed at the minimum and maximum event times.

knots(knots_list) defines the knot locations for the spline functions used to model the baseline
log cumulative hazard function. Must include boundary knots. Knots should be specified in increasing order.

Examples

For detailed examples, see reddooranalytics.se/products/merlin.

Author

Michael J. Crowther

Red Door Analytics

Stockholm, Sweden

michael@reddooranalytics.se

References

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. Stata Journal 2020;20(4):763-784.

Crowther MJ. Multilevel mixed effects parametric survival analysis: Estimation, simulation and application. Stata Journal 2019;19(4):931-949.

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.

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;40(9):2139-2154.