mspath {mspath} | R Documentation |
Ordinarily, fit a path-dependent (non-Markov) multi-state model by maximum likelihood. Observations of the process can be made at arbitrary times. Covariates can be fitted to the transition intensities. When the true state is observed with error, a model can be fitted which simultaneously estimates the state transition intensities and misclassification probabilities, with optional covariates on both processes.
This function will also perform some auxiliary operations, including
counting paths and generating data according to the model. See
do.what
.
Based on the Christopher Jackson's msm package v 0.3.3. The key
differences are that
this package can handle non-Markov models, and it does so using a
discrete-time approximation. Note that while some of the inputs to
msm were transition rates or misclassification probabilities, the
corresponding terms here are the intercepts of multinomial logits.
Finally, fromto
and related options do not exist; they are not
appropriate in a path-dependent setting and have also been dropped
from later versions of msm.
The interface to msm changed with more recent releases to a style that combined specification of allowed transitions or observation errors and their values. This would not be appropriate here, since 0 is a legitimate value for a coefficient. msm also has considerably richer options for handling the different kinds of information conveyed by the times of different observations. There is no theoretical problem with handling those options here; there just hasn't been a need to implement those options.
mspath ( formula, qmatrix, misc = FALSE, ematrix, inits, subject, covariates = NULL, constraint = NULL, misccovariates = NULL, miscconstraint = NULL, qconstraint=NULL, econstraint=NULL, pathvars = NULL, pathoffset = 0, pathconstraint = NULL, initprobs = NULL, data = list(), isexact = FALSE, fixedpars = NULL, stepnumerator = 1, stepdenominator = 1, do.what = 1, testing, comm, profile = FALSE, calcFactory = mspathCalculator, seed, trace, ... )
formula |
A formula giving the vectors containing
states and the corresponding observation times. For example,
states ~ times
Observed states should be in the set 1, ..., n , where
n is the number of states. Use state 0 —not
missing—when the state is not observed but the covariates are.
Note that analysis will exclude any observation with missing
covariates. |
qmatrix |
Matrix of indicators for the allowed transitions.
If a transition is allowed from state r to state s,
then qmatrix should have (r,s) entry 1, otherwise
it should have (r,s) entry 0. The diagonal of qmatrix
is ignored. For example,
rbind(
c( 0, 1, 1 ),
c( 1, 0, 1 ),
c( 0, 0, 0 )
)
represents a 'health - disease - death' model, with transitions allowed from health to disease, health to death, disease to health, and disease to death. |
misc |
Set misc = TRUE if misclassification between
observed and underlying states is to be modelled, FALSE if
observation is assumed accurate, or SIMPLE (no quotes)
if the values given in inits , in conjunction with
econstraint , specify the probabilities of measurement
error.
In the last case, those probabilities are fixed constants and you should omit misccovariates and
miscconstraint . All probabilities should be between 0 and 1,
as should 1 - their sum, which the routine computes as the
probability of accurate observation. For SIMPLE the values
are interpreted directly, rather than via a multinomial logit
transformation. You do not need to list them in fixedpars ;
the routine always considers them fixed.
|
ematrix |
(required when misc == TRUE ) Matrix of indicators for the allowed
misclassifications.
The rows represent underlying states, and the columns represent
observed states.
If an observation of state s is possible when the subject
occupies underlying state r, then ematrix should have
(r,s) entry 1, otherwise
it should have (r,s) entry 0. The diagonal of ematrix
is ignored. For example,
rbind(
c( 0, 1, 0 ),
c( 1, 0, 1 ),
c( 0, 1, 0 )
)
represents a model in which misclassifications are only permitted between adjacent states. |
inits |
(required) Vector of initial parameter estimates for the
optimisation. The “Details” section below gives full
interpretation of the parameters. Parameters are in this order:
- transition intercepts (reading across first rows of intensity matrix, then second row ... ) - covariate effects on transitions - history-dependent path variables effects on transitions - misclassification intercepts (reading across first row of misclassification matrix, then second row ...) - covariate effects on misclassification Covariate effects are given in the following order, - effects of first covariate on transition/misclassification matrix elements (reading across first row of matrix, then second row ...) - effects of second covariate ... |
subject |
Vector of subject identification numbers, when the data
are specified by formula . If missing, then all observations
are assumed to be on the same subject. These must be sorted so that
all observations on the same subject are adjacent. |
covariates |
Formula representing the covariates on the
transition intensities, for example,
~ age + sex + treatment
|
constraint |
A list of one vector for each named covariate. The
vector indicates which covariate effects are
constrained to be equal. Take, for example, a model with five
transition intensities and two covariates. Specifyingconstraint = list (age = c(1,1,1,2,2), treatment = c(1,2,3,4,5)) constrains the effect of age to be equal for the first three intensities, and equal for the fourth and fifth. The effect of treatment is assumed to be different for each intensity. Any vector of increasing numbers can be used as indicators. The intensity parameters are assumed to be ordered by reading across the rows of the transition matrix, starting at the first row. For categorical covariates, defined using factor(covname) ,
specify constraints as follows:list(..., covnameVALUE1 = c(...), covnameVALUE2 = c(...), ...) where VALUE1 , VALUE2 , ... are the levels of the factor.
Make sure the contrasts option is set appropriately, for
example, the default options(contrasts=c(contr.treatment,
contr.poly)) sets the first (baseline) level of unordered factors to
zero.
To assume no covariate effect on a certain transition, set its initial value to zero and use the fixedpars argument to fix
it during the optimisation.
|
misccovariates |
A formula representing the covariates on the
misclassification probabilities, analogously to covariates .
|
miscconstraint |
A list of one vector for each named covariate on
misclassification probabilities. The vector indicates which
covariate effects on misclassification probabilities are
constrained to be equal, analogously to constraint .
|
qconstraint |
A vector of indicators specifying which
transition intercepts are equal. For example,
qconstraint = c(1,2,3,3)
constrains the third and fourth transition intercepts to be equal, in a model with four allowed instantaneous transitions. |
econstraint |
A similar vector of indicators specifying which misclassification intercepts are constrained to be equal. |
pathvars |
A vector of the names of path-dependent variables to
use as covariates on transition intensities. Allowed values
currently include TIS , Time in State; TSO , Time Since
Origin, TIP , Time in Previous states, and LN(TIS) ,
LN(TIP) or LN(TSO) for natural logs of the same. |
pathoffset |
This value will be added to the origin time of path-dependent variables. Use it pick where in the interval to use for your time values, or to avoid calculations with 0. |
pathconstraint |
Constraints on path variables' effects on transition rates. |
initprobs |
Currently only used in misclassification
models. Vector of assumed underlying state occupancy
probabilities at each individual's first observation. Defaults to
c(1, rep(0, nstates-1)) , that is, in state 1 with a
probability of 1. The current implementation only supports
a single initial state.
|
data |
Optional data frame in which to interpret subject as
well as the time, state, and covariates defined in formula ,
covariates , and misccovariates formulae. |
isexact |
By default, the transitions of the Markov process
are assumed to take place at unknown occasions in between the
observation times. However, isexact = TRUE implies observation
times for entry to absorbing states are exact. Only paths with
transitions at the observed time will be counted. Typically, use this
for death. Note this option has no effect on non-absorbing states.
It is an error to use this option when observation of the absorbing
state is imprecise. |
fixedpars |
Vector of indices of parameters whose values will be
fixed at their initial values during the optimisation. These
correspond to indices of the inits vector, whose order is
specified above. This can be useful for building complex models
stage by stage. |
stepnumerator |
The calculation will use a time grid for each case
with steps every stepnumerator/stepdenominator . This allows
exact representation of rationals. The numerator and denominator
must both be integers. |
stepdenominator |
See stepnumerator just above. |
do.what |
1, the default, calculates a maximumlikelihood. To evaluate a
single likelihood, set all parameters to fixed.
0 counts the number of paths and related statistics without evaluating a likelihood. -1 get detailed counts, but not likelihoods, associated with each case. The return value is a matrix. 10 uses the model to simulate a random patch for each case, returning a data.frame with simulated observed states and times
and all other data as observed.
|
testing |
This argument is only for use by developers. Set it to a list with the expected low-level arguments to test if that's what is actually generated. The result will either be a stop or true; the underlying C code is not called. See the code for the exact arguments. |
comm |
An optional MPI communicator. If you provide it, perform a
distributed calculation. Ordinarily, comm = 0 . Ordinarily a
call using this argument will be invoked from the commands supplied
via master 's channel argument. See that
routine for more information on distributed calculation and the
necessary setup.
|
profile |
Profile the time taken in the distributed calculation. |
calcFactory |
Optional, advanced argument. This should be a
function that, when called with the arguments to the standard
mspathCalculator , produces a calculator for internal
use. The comm and profile arguments will be ignored
in this case. |
seed |
Optional seed to set for simulation (do.what=10 ).
Use an integer for full precision; it will be cast to
unsigned. |
trace |
Optional list of arguments to checkpoint .
If present, the objective function will be wrapped by
checkpoint and every function evaluation's parameters
and value will be recorded in the caller's environment and in files
on disk. See checkpoint for details. |
... |
Optional arguments to the general-purpose R
optimization routine optim . Useful options include
method="BFGS" for using a quasi-Newton optimisation
algorithm, which can often be faster than the default Nelder-Mead.
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 and the standard
errors cannot be calculated due to 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. |
For each case, the model generates all possible paths consistent with
the observed data, measurement model, and transition model. It
computes and sums the likelihood of each path, including the initial
probability specified in initprobs
.
It is quite possible to specify a model which makes the observed data impossible; you will get an error in that case.
Both the transition and misclassification models use the multinomial logit. Self transition has assumed coefficients of 0; the other coefficients are multiplied by the covariates, summed and exponentioned. All coefficients, including the intercept or baseline terms, are so transformed. With J possible outcomes and p(j) the probability of the j'th outcome, the formula is
p(j) = exp[X*b(j)]/sum{exp[X*b(k)], k=0 to J}.
b(j) is the vector of coefficients for outcome j and X are the covariates. Calculations use the covariates at the endpoint of the transition to set the transition rate.
TIP
, Time in Previous states, is the start time of the current
state. Note it is time in all previous state{s}, not just the
previous state. It will remain constant as TIS
, Time in State,
changes. If pathOffset
=0, TSO = TIP + TIS
. TIP
is a constant for the first state; we recommend constraining the
corresponding coefficient to 0 if you allow effects of TIP
to
vary with stage.
There must be enough information in the data on each state to estimate each transition rate, otherwise the likelihood will be flat and the maximum will not be found. It may be appropriate to reduce the number of states in the model, or reduce the number of covariate effects, to ensure convergence.
Choosing an appropriate set of initial values for the optimisation can also be important. For flat likelihoods, 'informative' initial values will often be required.
If a single discrete time point includes multiple observations, one is selected. In the initial interval, the first point is used to get the correct start time. For later times, selection favors observations with (non-0) state, and picks the last possible observation. If any observations include state, the last among those is the one picked.
Ordinarily a mspath
object. See it for
details. However, there are the following special cases:
|
Either returns TRUE or stops
with an error if the calculated variables do not those in the
testing list. Specifying testing supersedes any
actions given in do.what . |
|
Return a matrix with counts for each case. |
|
Return a data frame with a random path as described above. |
do.what=10
requests generation of a simulated data.frame
based on the specified model and data.
The first step of the simulation drops rows as in the likelihood
estimation. Rows with missing data may be dropped depending upon the
options you have requested globally or in data
. Then rows in
the same discrete time interval are dropped.
Second, all the data are collected: data
and subject
and
any variables used in formula
, covariates
and
misccovariates
. Data not present in data
will be pulled
in from the appropriate environment and added to the output
data.frame
. If this happens with subject
, the name of
the column will be "subject"
; in these circumstances
subject
may not appear as a term in any formula. Note that
variables in data
that are not referenced in formulae
will be retained for later analysis.
The variable names will be those used in the formulae. If the caller
specifies subject
with a name (e.g., mspath(subject=id,
...)
) that name will be the name of the subject variable;
otherwise the name will be "subject"
. In either case, that
name should not be used elsewhere.
Third, for each case a true path and corresponding observed states
will be simulated using the model. Ordinarily the simulated data will
have the same observation times as the real data. However, if the
simulated path enters an absorbing state and isexact == TRUE
then the simulated data will show an observation at that time, using
the then current covariates. Furthermore, the simulated data will
never have more than one observation in an absorbing state, since
otherwise mspath
will refuse to estimate for the simulated
data.
The trace
argument will cause tracing information to go
into the caller's environment and onto the disk, as described in
checkpoint
.
Ross D. Boylan ross@biostat.ucsf.edu
Boylan, R.D. 2010 Multi-state Path-Dependent Models. Working paper. University of California, San Francisco. Available as doc/mspath-program.* within this package.
Bacchetti, P., Boylan, R.D., Terrault, N.A., Monto, A. and Berenguer, Marina. Forthcoming Non-Markov multistate modeling using time-varying covariates, with application to progression of liver fibrosis due to hepatitis C following liver transplant. The International Journal of Biostatistics.
Jackson, C.H., Sharples, L.D., Thompson, S.G. and Duffy, S.W. and Couto, E. 2003 Multi-state Markov models for disease progression with classification error. The Statistician, 52(2), 193–209.
Jackson, C.H. and Sharples, L.D. 2002 Hidden Markov models for the onset and progresison of bronchiolitis obliterans syndrome in lung transplant recipients Statistics in Medicine, 21(1), 113–128.
Kay, R. 1986 A Markov model for analysing cancer markers and disease states in survival studies. Biometrics 42, 855–865.
Gentleman, R.C., Lawless, J.F., Lindsey, J.C. and Yan, P. 1994 Multi-state Markov models for analysing incomplete disease history data with illustrations for HIV disease. Statistics in Medicine 13(3), 805–821.
Satten, G.A. and Longini, I.M. 1996 Markov chains with measurement error: estimating the 'true' course of a marker of the progression of human immunodeficiency virus disease (with discussion) Applied Statistics 45(3), 275-309.
master
, checkpoint
,
mspath
return object.