mspath {mspath}R Documentation

Multi-state Path-Dependent Models

Description

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.

Usage

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,
      ... )

Arguments

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 0not 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. Specifying

constraint = 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.

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.

Value

Ordinarily a mspath object. See it for details. However, there are the following special cases:

testing non-empty 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.
do.what=-1 Return a matrix with counts for each case.
do.what=10 Return a data frame with a random path as described above.

Simulation

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.

Note

The trace argument will cause tracing information to go into the caller's environment and onto the disk, as described in checkpoint.

Author(s)

Ross D. Boylan ross@biostat.ucsf.edu

References

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.

See Also

master, checkpoint, mspath return object.


[Package mspath version 0.9-9 Index]