In a mixture model for competing events, an individual can experience one of a set of different events. We specify a model for the probability that they will experience each event before the others, and a model for the time to the event conditionally on that event occurring first.
flexsurvmix(
formula,
data,
event,
dists,
pformula = NULL,
anc = NULL,
partial_events = NULL,
initp = NULL,
inits = NULL,
fixedpars = NULL,
dfns = NULL,
method = "direct",
em.control = NULL,
optim.control = NULL,
aux = NULL,
sr.control = survreg.control(),
integ.opts,
hess.control = NULL,
...
)List of objects containing information about the fitted model. The
important one is res, a data frame containing the parameter
estimates and associated information.
Survival model formula. The left hand side is a Surv
object specified as in flexsurvreg. This may define various
kinds of censoring, as described in Surv. Any covariates on
the right hand side of this formula will be placed on the location
parameter for every component-specific distribution. Covariates on other
parameters of the component-specific distributions may be supplied using
the anc argument.
Alternatively, formula may be a list of formulae, with one
component for each alternative event. This may be used to specify
different covariates on the location parameter for different components.
A list of formulae may also be used to indicate that for particular individuals, different events may be observed in different ways, with different censoring mechanisms. Each list component specifies the data and censoring scheme for that mixture component.
For example, suppose we are studying people admitted to hospital,and the competing states are death in hospital and discharge from hospital. At time t we know that a particular individual is still alive, but we do not know whether they are still in hospital, or have been discharged. In this case, if the individual were to die in hospital, their death time would be right censored at t. If the individual will be (or has been) discharged before death, their discharge time is completely unknown, thus interval-censored on (0,Inf). Therefore, we need to store different event time and status variables in the data for different alternative events. This is specified here as
formula = list("discharge" = Surv(t1di, t2di, type="interval2"),
"death" = Surv(t1de, status_de))
where for this individual, (t1di, t2di) = (0, Inf) and (t1de,
status_de) = (t, 0).
The "dot" notation commonly used to indicate "all remaining variables" in a
formula is not supported in flexsurvmix.
Data frame containing variables mentioned in formula,
event and anc.
Variable in the data that specifies which of the alternative
events is observed for which individual. If the individual's follow-up is
right-censored, or if the event is otherwise unknown, this variable must
have the value NA.
Ideally this should be a factor, since the mixture components can then be
easily identified in the results with a name instead of a number. If this
is not already a factor, it is coerced to one. Then the levels of the
factor define the required order for the components of the list arguments
dists, anc, inits and dfns. Alternatively, if
the components of the list arguments are named according to the levels of
event, then the components can be arranged in any order.
Vector specifying the parametric distribution to use for each
component. The same distributions are supported as in
flexsurvreg.
Formula describing covariates to include on the component membership proabilities by multinomial logistic regression. The first component is treated as the baseline.
The "dot" notation commonly used to indicate "all remaining variables" in a formula is not supported.
List of component-specific lists, of length equal to the number of components. Each component-specific list is a list of formulae representing covariate effects on parameters of the distribution.
If there are covariates for one component but not others, then a list
containing one null formula on the location parameter should be supplied
for the component with no covariates, e.g list(rate=~1) if the
location parameter is called rate.
Covariates on the location parameter may also be supplied here instead of
in formula. Supplying them in anc allows some components
but not others to have covariates on their location parameter. If a covariate
on the location parameter was provided in formula, and there are
covariates on other parameters, then a null formula should be included
for the location parameter in anc, e.g list(rate=~1)
List specifying the factor levels of event
which indicate knowledge that an individual will not experience particular
events, but may experience others. The names of the list indicate codes
that indicate partial knowledge for some individuals. The list component
is a vector, which must be a subset of levels(event) defining the
events that a person with the corresponding event code may experience.
For example, suppose there are three alternative events called
"disease1","disease2" and "disease3", and for some
individuals we know that they will not experience "disease2", but
they may experience the other two events. In that case we must create a
new factor level, called, for example "disease1or3", and set the
value of event to be "disease1or3" for those individuals.
Then we use the "partial_events" argument to tell
flexsurvmix what the potential events are for individuals with this
new factor level.
partial_events = list("disease1or3" = c("disease1","disease3"))
Initial values for component membership probabilities. By default, these are assumed to be equal for each component.
List of component-specific vectors. Each component-specific
vector contains the initial values for the parameters of the
component-specific model, as would be supplied as the inits argument of
flexsurvreg. By default, a heuristic is used to obtain
initial values, which depends on the parametric distribution being used,
but is usually based on the empirical mean and/or variance of the survival
times.
Indexes of parameters to fix at their initial values and
not optimise. Arranged in the order: baseline mixing probabilities,
covariates on mixing probabilities, time-to-event parameters by mixing
component. Within mixing components, time-to-event parameters are ordered
in the same way as in flexsurvreg.
If fixedpars=TRUE then all parameters will be fixed and the
function simply calculates the log-likelihood at the initial values.
Not currently supported when using the EM algorithm.
List of lists of user-defined distribution functions, one for
each mixture component. Each list component is specified as the
dfns argument of flexsurvreg.
Method for maximising the likelihood. Either "em" for
the EM algorithm, or "direct" for direct maximisation.
List of settings to control EM algorithm fitting. The only options currently available are
trace set to 1 to print the parameter estimates at each iteration
of the EM algorithm
reltol convergence criterion. The algorithm stops if the log
likelihood changes by a relative amount less than reltol. The
default is the same as in optim, that is,
sqrt(.Machine$double.eps).
var.method method to compute the covariance matrix. "louis"
for the method of Louis (1982), or "direct"for direct numerical
calculation of the Hessian of the log likelihood.
optim.p.control A list that is passed as the control
argument to optim in the M step for the component membership
probability parameters. The optimisation in the M step for the
time-to-event parameters can be controlled by the optim.control
argument to flexsurvmix.
For example, em.control = list(trace=1, reltol=1e-12).
List of options to pass as the control argument
to optim, which is used by method="direct" or in the
M step for the time-to-event parameters in method="em". By
default, this uses fnscale=10000 and ndeps=rep(1e-06,p)
where p is the number of parameters being estimated, unless the
user specifies these options explicitly.
A named list of other arguments to pass to custom distribution
functions. This is used, for example, by flexsurvspline to
supply the knot locations and modelling scale (e.g. hazard or odds). This
cannot be used to fix parameters of a distribution --- use
fixedpars for that.
For the models which use survreg to find the
maximum likelihood estimates (Weibull, exponential, log-normal), this list
is passed as the control argument to survreg.
List of named arguments to pass to
integrate, if a custom density or hazard is provided without
its cumulative version. For example,
integ.opts = list(rel.tol=1e-12)
List of options to control covariance matrix computation. Available options are:
numeric. If TRUE then numerical methods are used
to compute the Hessian for models where an analytic Hessian is
available. These models include the Weibull (both versions),
exponential, Gompertz and spline models with hazard or odds
scale. The default is to use the analytic Hessian for these
models. For all other models, numerical methods are always used
to compute the Hessian, whether or not this option is set.
tol.solve. The tolerance used for solve
when inverting the Hessian (default .Machine$double.eps)
tol.evalues The accepted tolerance for negative
eigenvalues in the covariance matrix (default 1e-05).
The Hessian is positive definite, thus invertible, at the maximum
likelihood. If the Hessian computed after optimisation convergence can't
be inverted, this is either because the converged result is not the
maximum likelihood (e.g. it could be a "saddle point"), or because the
numerical methods used to obtain the Hessian were inaccurate. If you
suspect that the Hessian was computed wrongly enough that it is not
invertible, but not wrongly enough that the nearest valid inverse would be
an inaccurate estimate of the covariance matrix, then these tolerance
values can be modified (reducing tol.solve or increasing
tol.evalues) to allow the inverse to be computed.
Optional arguments to the general-purpose optimisation routine
optim. For example, the BFGS optimisation algorithm is the
default in flexsurvreg, but this can be changed, for example
to method="Nelder-Mead" which can be more robust to poor initial
values. If the optimisation fails to converge, consider normalising the
problem using, for example, control=list(fnscale = 2500), for
example, replacing 2500 by a number of the order of magnitude of the
likelihood. If 'false' convergence is reported with a
non-positive-definite Hessian, then consider tightening the tolerance
criteria for convergence. If the optimisation takes a long time,
intermediate steps can be printed using the trace argument of the
control list. See optim for details.
This differs from the more usual "competing risks" models, where we specify "cause-specific hazards" describing the time to each competing event. This time will not be observed for an individual if one of the competing events happens first. The event that happens first is defined by the minimum of the times to the alternative events.
The flexsurvmix function fits a mixture model to data consisting of a
single time to an event for each individual, and an indicator for what type
of event occurs for that individual. The time to event may be observed or
censored, just as in flexsurvreg, and the type of event may be
known or unknown. In a typical application, where we follow up a set of
individuals until they experience an event or a maximum follow-up time is
reached, the event type is known if the time is observed, and the event type
is unknown when follow-up ends and the time is right-censored.
The model is fitted by maximum likelihood, either directly or by using an
expectation-maximisation (EM) algorithm, by wrapping
flexsurvreg to compute the likelihood or to implement the E
and M steps.
Some worked examples are given in the package vignette about multi-state
modelling, which can be viewed by running vignette("multistate", package="flexsurv").
Jackson, C. H. and Tom, B. D. M. and Kirwan, P. D. and Mandal, S. and Seaman, S. R. and Kunzmann, K. and Presanis, A. M. and De Angelis, D. (2022) A comparison of two frameworks for multi-state modelling, applied to outcomes after hospital admissions with COVID-19. Statistical Methods in Medical Research 31(9) 1656-1674.
Larson, M. G., & Dinse, G. E. (1985). A mixture model for the regression analysis of competing risks data. Journal of the Royal Statistical Society: Series C (Applied Statistics), 34(3), 201-211.
Lau, B., Cole, S. R., & Gange, S. J. (2009). Competing risk regression models for epidemiologic data. American Journal of Epidemiology, 170(2), 244-256.