msm: Multi-state Markov and hidden Markov models in continuous time

Description

Fit a continuous-time Markov or hidden Markov multi-state model by maximum likelihood. Observations of the process can be made at arbitrary times, or the exact times of transition between states can be known. Covariates can be fitted to the Markov chain transition intensities or to the hidden Markov observation process.

Usage

msm ( formula, subject=NULL, data = list(), qmatrix, gen.inits = FALSE,
      ematrix=NULL, hmodel=NULL, obstype=NULL,
      covariates = NULL, covinits = NULL, constraint = NULL,
      misccovariates = NULL, misccovinits = NULL, miscconstraint = NULL,
      hcovariates = NULL, hcovinits = NULL, hconstraint = NULL,
      qconstraint=NULL, econstraint=NULL, initprobs = NULL,
      death = FALSE, exacttimes = FALSE, censor=NULL,
      censor.states=NULL, cl = 0.95, fixedpars = NULL, center=TRUE,
      hessian=TRUE, ... )

Arguments

formula

A formula giving the vectors containing the observed states and the corresponding observation times. For example, state ~ time

Observed states should be in the set 1, ..., n, where n is the numb

subject

Vector of subject identification numbers for the data 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.

data

Optional data frame in which to interpret the variables supplied in formula, subject, covariates, misccovariates, hcovariates and obstype.

qmatrix

Initial transition intensity matrix of the Markov chain. If an instantaneous transition is not allowed from state $r$ to state $s$, then qmatrix should have $(r,s)$ entry 0, otherwise it should be non-zero. Any diagonal ent

gen.inits

If TRUE, then initial values for the transition intensities are estimated by assuming that the data represent the exact transition times of the process. The non-zero entries of the supplied qmatrix are assumed to indi

ematrix

If misclassification between states is to be modelled, this should be a matrix of initial values for the misclassification probabilities. The rows represent underlying states, and the columns represent observed states. If an observation of

hmodel

Specification of the hidden Markov model. This should be a list of return values from the constructor functions described in the hmm-dists help page. Each element of the list corresponds

obstype

A vector specifying the observation scheme for each row of the data. This can be included in the data frame data along with the state, time, subject IDs and covariates. Its elements should be either 1, 2 or 3, meaning as follows

covariates

Formula representing the covariates on the transition intensities via a log-linear model. For example,

~ age + sex + treatment

covinits

Initial values for log-linear effects of covariates on the transition intensities. This should be a named list with each element corresponding to a covariate. A single element contains the initial values for that covariate on each transition

constraint

A list of one vector for each named covariate. The vector indicates which covariate effects on intensities are constrained to be equal. Take, for example, a model with five transition intensities and two covariates. Specifying constraint = lis

misccovariates

A formula representing the covariates on the misclassification probabilities, analogously to covariates. Only used if the model is specified using ematrix, rather than hmodel.

misccovinits

Initial values for the covariates on the misclassification probabilities, defined in the same way as covinits. Only used if the model is specified using ematrix.

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. Only use

hcovariates

List of formulae the same length as hmodel, defining any covariates governing the hidden Markov outcome models. The covariates operate on a suitably link-transformed linear scale, for example, log scale for a Poisson outcome mode

hcovinits

Initial values for the hidden Markov model covariates. A list of the same length as hcovariates. Each element is a vector with initial values for each covariate on that state. For example, the above hcovariates can b

hconstraint

A named list. Each element is a vector of constraints on the named hidden Markov model parameter. The vector has length equal to the number of times that class of parameter appears in the whole model.

For example consider the three-state hid

qconstraint

A vector of indicators specifying which baseline transition intensities are equal. For example, qconstraint = c(1,2,3,3)

constrains the third and fourth intensities to be equal, in a model with four allowed instantaneous

econstraint

A similar vector of indicators specifying which baseline misclassification probabilities are constrained to be equal. Only used if the model is specified using ematrix, rather than hmodel.

initprobs

Currently only used in hidden Markov 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.

death

Vector of indices of the death states. A death state is an absorbing state whose time of entry is known exactly, but the individual is assumed to be in an unknown transient state ("alive") at the previous instant. This is the usual situation

censor

A state, or vector of states, which indicates censoring. Censoring means that the observed state is known only to be one of a particular set of states. For example, censor=999 indicates that all observations of 999 i

censor.states

Specifies the underlying states which censored observations can represent. If censor is a single number (the default) this can be a vector, or a list with one element. If censor is a vector with more than one element

exacttimes

By default, the transitions of the Markov process are assumed to take place at unknown occasions in between the observation times. If exacttimes is set to TRUE, then all observation times are assumed to represent

Width of symmetric confidence intervals for maximum likelihood estimates, by default 0.95.

fixedpars

Vector of indices of parameters whose values will be fixed at their initial values during the optimisation. These are given in the order: transition intensities (reading across rows of the transition matrix), covariates on intensities (ordered

center

If TRUE (the default) then covariates are centered at their means during the maximum likelihood estimation. This usually improves convergence.

hessian

If TRUE (the default) then the Hessian matrix is computed at the maximul likelihood estimates, to obtain standard errors.

...

Optional arguments to the general-purpose Roptimization routine optim. Useful options include method="BFGS" for using a quasi-Newton optimisation algorithm, which can often be faster than

Value

A list of class msm, with components:
callThe original call to msm.
QmatricesA list of matrices. The first component, labelled logbaseline, is a matrix containing the estimated transition intensities on the log scale with any covariates fixed at their means in the data. Each remaining component is a matrix giving the linear effects of the labelled covariate on the matrix of log intensities. To extract an estimated intensity matrix on the natural scale, at an arbitrary combination of covariate values, use the function qmatrix.msm.
QmatricesSEThe standard error matrices corresponding to Qmatrices.
QmatricesL,QmatricesUCorresponding lower and upper symmetric confidence limits, of width 0.95 unless specified otherwise by the cl argument.
EmatricesA list of matrices. The first component, labelled logitbaseline, is the estimated misclassification probability matrix with any covariates fixed at their means in the data. Each remaining component is a matrix giving the linear effects of the labelled covariate on the matrix of logit misclassification probabilities. To extract an estimated misclassification probability matrix on the natural scale, at an arbitrary combination of covariate values, use the function ematrix.msm.
EmatricesSEThe standard error matrices corresponding to Ematrices.
EmatricesL,EmatricesUCorresponding lower and upper symmetric confidence limits, of width 0.95 unless specified otherwise by the cl argument.
sojournA list with components:
mean = estimated mean sojourn times in the transient states, with covariates fixed at their means.
se = corresponding standard errors.
minus2loglikMinus twice the maximised log-likelihood.
estimatesVector of untransformed maximum likelihood estimates returned from optim. Transition intensities are on the log scale and misclassification probabilities are the logit scale.
estimates.tVector of transformed maximum likelihood estimates with intensities and probabilities on their natural scales.
fixedparsIndices of estimates which were fixed during the maximum likelihood estimation.
covmatCovariance matrix corresponding to estimates.
ciMatrix of confidence intervals corresponding to estimates
foundseLogical value indicating whether the Hessian was positive-definite at the supposed maximum of the likelihood. If not, the covariance matrix of the parameters is unavailable. In these cases the optimisation has probably not converged to a maximum.
dataA list of constants and vectors giving the data, for use in post-processing.
qmodelA list of objects specifying the model for transition intensities, for use in post-processing.
emodelA list of objects specifying the model for misclassification.
qcmodelA list of objects specifying the model for covariates on the transition intensities.
ecmodelA list of objects specifying the model for covariates on misclassification probabilities.
hmodelA list of class "hmodel", containing objects specifying the hidden Markov model. Estimates of "baseline" location parameters are presented with any covariates fixed to their means in the data.
cmodelA list of objects specifying any model for censoring.
Printing a msm object by typing the object's name at the command line implicitly invokes print.msm. This formats and prints the important information in the model fit. This includes the fitted transition intensity matrix, matrices containing covariate effects on intensities, and mean sojourn times from a fitted msm model. When there is a hidden Markov model, the chief information in the hmodel component is also formatted and printed. This includes estimates and confidence intervals for each parameter.
To extract summary information from the fitted model, it is recommended to use the more flexible extractor functions, such as qmatrix.msm, pmatrix.msm, sojourn.msm, instead of directly reading from list components of msm objects.

Details

For full details about the methodology behind the msm package, refer to the PDF manual msm-manual.pdf in the doc subdirectory of the package. This includes a tutorial in the typical use of msm.

Users upgrading from versions of msm less than 0.5 will need to change some of their model fitting syntax. In particular, initial values are now specified in the qmatrix and covinits arguments instead of inits, and qmatrix is no longer a matrix of 0/1 indicators. See the appendix to the PDF manual or the NEWS file in the top-level installation directory for a full list of changes.

For simple multi-state Markov models, the likelihood is calculated in terms of the transition intensity matrix $Q$. When the data consist of observations of the Markov process at arbitrary times, the exact transition times are not known. Then the likelihood is calculated using the transition probability matrix $P(t) = \exp(tQ)$, where $\exp$ is the matrix exponential. If state $i$ is observed at time $t$ and state $j$ is observed at time $u$, then the contribution to the likelihood from this pair of observations is the $i,j$ element of $P(u - t)$. See, for example, Kalbfleisch and Lawless (1985), Kay (1986), or Gentleman et al. (1994).

For hidden Markov models, the likelihood for an individual with $k$ observations is calculated directly by summing over the unknown state at each time, producing a product of $k$ matrices. The calculation is a generalisation of the method described by Satten and Longini (1996), and also by Jackson and Sharples (2002), and Jackson et al. (2003).

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. Hidden Markov models are particularly susceptible to non-identifiability, especially when combined with a complex transition matrix.

Choosing an appropriate set of initial values for the optimisation can also be important. For flat likelihoods, 'informative' initial values will often be required.

References

Kalbfleisch, J., Lawless, J.F., The analysis of panel data under a Markov assumption Journal of the Americal Statistical Association (1985) 80(392): 863--871. Kay, R. A Markov model for analysing cancer markers and disease states in survival studies. Biometrics (1986) 42: 855--865. Gentleman, R.C., Lawless, J.F., Lindsey, J.C. and Yan, P. Multi-state Markov models for analysing incomplete disease history data with illustrations for HIV disease. Statistics in Medicine (1994) 13(3): 805--821. Satten, G.A. and Longini, I.M. 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 (1996) Jackson, C.H. and Sharples, L.D. Hidden Markov models for the onset and progression of bronchiolitis obliterans syndrome in lung transplant recipients Statistics in Medicine, 21(1): 113--128 (2002).

Jackson, C.H., Sharples, L.D., Thompson, S.G. and Duffy, S.W. and Couto, E. Multi-state Markov models for disease progression with classification error. The Statistician, 52(2): 193--209 (2003)

Examples

Run this code

### Heart transplant data
### For further details and background to this example, see
### the PDF manual in the doc directory. 
data(heart)
print(heart[1:10,])
twoway4.q <- rbind(c(-0.5, 0.25, 0, 0.25), c(0.166, -0.498, 0.166, 0.166),
c(0, 0.25, -0.5, 0.25), c(0, 0, 0, 0))
statetable.msm(state, PTNUM, data=heart)
crudeinits.msm(state ~ years, PTNUM, data=heart, qmatrix=twoway4.q)
heart.msm <- msm( state ~ years, subject=PTNUM, data = heart, 
                 qmatrix = twoway4.q, death = 4,
                 control = list ( trace = 2, REPORT = 1 )  )
heart.msm
qmatrix.msm(heart.msm)
pmatrix.msm(heart.msm, t=10)
sojourn.msm(heart.msm)

Run the code above in your browser using DataLab