Estimating parametric semi-Markov models from panel data using phase-type approximations

Andrew Titman has a new paper in Statistics and Computing. This extends previous work on fitting semi-Markov models to panel data using phase-type distributions. Here, rather than assume a model in which each state is assumed to have a Coxian phase-type sojourn distribution, the model assumes a standard parametric sojourn distribution (e.g. Weibull or Gamma). The computational tractability of phase-type distributions is exploited by approximating the parametric semi-Markov model by a model in which each state has a 5-phase phase-type distribution. In order to achieve this, a family of approximations to the parametric distribution, with scale parameter 1, is developed by solving a relatively large one-off optimization assuming the optimal phase-type parameters for given shape parameters evolve as B-spline functions. The resulting phase-type approximations can then be scaled to give an approximation for any combination of scale and shape parameter and then embedded into the overall semi-Markov process. The resulting approximate likelihood appears to be very close to the exact likelihood, both in terms of shape and magnitude.

In general a 2-phase Coxian phase-type model is likely to give similar results to a Weibull or Gamma model. The only advantages of the Weibull or Gamma model are that the parameters are identifiable under the null (Markov) model so standard likelihood ratio tests can be used to compare models (unlike for the 2-phase model). Also the Weibull or Gamma model requires fewer parameters so may be useful for smaller sample sizes.

Effect of vitamin A deficiency on respiratory infection: Causal inference for a discretely observed continuous time non-stationary Markov process

Mingyuan Zhang and Dylan Small have a paper to appear in The Canadian Journal of Statistics, currently available here. The paper uses a multi-state model approach to obtain estimates of the causal effects of vitamin D deficiency on respiratory infection.

The observed data consist of observations of respiratory infection status, vitamin D deficiency status and whether the child is stunted at time t. Each of these is a binary variable, leading to 8 possible observation patterns.

The data are assumed to be generated from an underlying non-homogeneous Markov chain on 32 states, consisting of a latent 4-level definition of vitamin D deficiency, the stunting variable, the observed respiratory infection status and additionally a counter-factual respiratory infection status defined as the status at time t hat would have occurred had the child maintained the lowest level of vitamin D deficiency from time 0 to time t.

Let represent the observed infection status, the underlying vitamin deficiency and the counter-factual infection status, the assumed relationship between them is given by , and for j>0. then measures the additional risk of having a respiratory infection at a particular time, given a current vitamin D deficiency at level j>0.

The overall model is a pretty innovative use of a hidden Markov model structure to obtain those causal estimates. The true process is assumed to occur in continuous time. However, it is desired that the underlying transition intensities are not time constant. As a result, the authors choose to approximate the process by one in discrete time (with some similarities to the approach of Bacchetti et al 2010).

In practical terms, the weakness of the model seems to be the assumption that the relative effect of a current vitamin deficiency compared to a perfect record of vitamin D levels, is both constant in time and does not depend on the past history of vitamin D deficiency. The latter assumption is essentially the Markov assumption and would be quite difficult to relax. The observed infection status is effectively a misclassified version of the counter-factual infection status. As a result the former assumption could be relaxed by letting the depend on time, perhaps in a piecewise constant fashion.

The definition of the process as initially being continuous is a little artificial and ill-specified in places. For instance, the transition intensities are defined in terms of a logit transformation from the outset. Also, it is stated that the underlying 32 state process (including both of ) is a continuous-time Markov process. However if is defined only through the misclassification equations, there is no limiting intensity for as . Once the process is in discrete time this is not a problem, but it would be more sensible to define the underlying process in continuous time to be 16 states not including and then specify that the observed are observations in a hidden Markov model.

Mixture distributions in multi-state modelling: Some considerations in a study of psoriatic arthritis

Aidan O'Keeffe, Brian Tom and Vern Farewell have a new paper in Statistics in Medicine. This considers random effects models for clustered multi-state models, specifically considering the psoriatic arthritis example considered in their previous paper. The particular emphasis in the current paper is comparing models using a continuous gamma frailty term with an extended model that additionally allows a "stayer" component. The latter can be thought of as a joining of the continuous random effects model (e.g. Cook et al 2004) and the "mover-stayer" model (e.g. Cook et al 2002). In a discussion on random effects, particularly where there is a question of finite mass points, it is strange there is no mention of the non-parametric mixing distribution (Laird, JASA 1978). While this hasn't really been considered for continuous time processes (the nearest example is Frydman's model) it wouldn't be particularly hard to implement at least the (not entirely reliable) EM algorithm type approach used in discrete-time by Maruotti and Rocci. The apparent presence or absence of a "stayer" component is likely to be heavily dependent on the parametric assumptions made about the rest of the mixing distribution. A very small random effect is indistinguishable from a zero random effect. The authors do emphasize the need to consider several possible mover-stayer models and consider the biological plausibility of them. It is also worth mentioning that all these issues will also hinge on the appropriateness of other assumes e.g. conditionally time homogeneous Markov processes.

Fitting and Interpreting Continuous-Time Latent Markov Models for Panel Data

Jane Lange and Vladimir Minin have a paper currently available as a University of Washington Working paper. This considers the computational aspects of fitting a latent continuous time Markov model to multi-state interval censored or panel data. The latent continuous time Markov model proposed is essentially a generalization of the phase-type semi-Markov models considered by Titman and Sharples. It is assumed that for each observable state r, there exist latent states: Observation in state r can correspond to occupation into any one of these latent states. This induces a type of hidden Markov model and can allow the observed process to be semi-Markov. Lange and Minin propose an EM algorithm to fit the model and compare this with other possible optimization techniques, concluding that their EM algorithm generally performs better, both in terms of reliability and speed of convergence. On this issue, it might have been worth considering the performance of, say, L-BFGS-B (or other constrained optimization methods) using the natural [0,Inf) range for the transition intensities, because many convergence problems for BFGS are due to cases where the true value is near (or at) zero and the exponential parameterization means it never gets there. **Update: This paper is now published in Statistics in Medicine.

Optimal designs for epidemiologic longitudinal studies with binary outcomes

Juha Mehtala ,Kari Auranen and Sangita Kulathina have a new paper in Statistical Methods in Medical Research. This considers optimal design of panel studies for binary stochastic processes which are assumed to be time-homogeneous Markov. Optimality is assessed based on the trace of the expected Fisher information. The authors also consider a two-phase design where the time spacing is improved upon after the first stage.

This paper, and also the paper by Hwang and Brookmeyer, restrict attention to equally spaced designs. Potentially, there may be efficiency gains from having unequally spaced observations particularly if equilibrium is not assumed at time 0. Moreover, an adaptive "doctor's care" type design where the gap to the next clinic visit depends on the current state is also likely to give some efficiency gains. There has been relatively little work on design for multi-state models. Clearly the efficiency gains achieved in practice depend on the accuracy of the initial assumptions about the true parameter values, although if a two-phase approach is possible it can militate against this. In complicated studies a simple relationship between expected parameter values and optimal design may not exist but may nevertheless be calculable by bespoke numerical optimization. Perhaps the main problem is that studies are rarely designed with an eventual multi-state model analysis in mind.

Intermittent observation of time-dependent explanatory variables: a multistate modelling approach

Brian Tom and Vern Farewell have a new paper in Statistics in Medicine. This considers the problem of estimating the effect of a time dependent covariate on a multi-state process when both the disease process of interest and the time dependent covariate are only intermittently observed. The most common existing approach to dealing with this problem is to assume that the time dependent covariate is constant between observations, taking the last observed value. The authors instead jointly model the two processes as an expanded multi-state model, if the disease process had n states and the covariate process m states the resulting process will have states. An additional assumption, that movements in the covariate process are not directly affected by the state of the disease process, is also made.

A simulation study is performed which shows that the approach of assuming the time dependent covariate is constant leads to biased estimates, particularly when there is a bias in the trend of the covariate process (e.g. much more likely to decrease than increase in value).

The overall approach taken by the authors is to model their way out of difficulty. They assume that both the disease process of interest and the covariate process are jointly time homogeneous Markov and the validity of the results will depend on these assumptions being correct. As noted by the authors, if the covariate can take more than a small number of values the approach becomes unattractive because of the large number of nuisance parameters required. A point not really emphasized, but related to the analogous approach taken by Cortese and Andersen for continuously observed competing risks data (bizarrely not referenced in this paper despite massive relevance!), is that having modelled the time dependent covariate, the model can then be used to make overall predictions.

One could argue that the convention of following forward the covariate value observed from the previous period is a way of allowing a prediction to be made about the trajectory in the next period. A fairer comparison in some cases might therefore be to look at the bias in estimating the transition probabilities to time given a covariate variate value of at time . While we would expect these estimates still to be biased, the amount of bias is likely to be less than found by looking at the regression coefficients directly.

An open problem seems to be the development of methods that do not require strong assumptions, or else are robust to misspecification, to deal with intermittently observed time dependent covariates.

A case-study in the clinical epidemiology of psoriatic arthritis: multistate models and causal arguments

Aidan O'Keeffe, Brian Tom and Vern Farewell have a new paper in Applied Statistics. They model data on the status of 28 individual hand joints in psoriatic arthitis patients. The damage to joints for a particular patient is expected to be correlated. Particular questions of interest are whether the damage process has 'symmetry', meaning if a specific joint in one hand is damaged does this increase the hazard of the corresponding joint in the other hand becoming damaged, and whether activity, meaning whether a joint is inflamed, causes the joint damage.

Two approaches to analysing the data are taken. In the first each pair of joints (from the left and right hands) is modelled by a 4 state model where the initial state is no damage to either joint, the terminal state is damage to both joints and progession to the terminal state may be through either a state representing damage only to the left joint or through a state representing damage only to the right joint. The correlation between joints within a patient is incorporated by allowing a common Gamma frailty which affects all transition intensities for all 14 pairs of joints. Symmetry can then be assessed by comparing the hazard of damage to the one joint with or without damage having occurred to the other joint. Under this approach, activity is incorporated only as a time dependent covariate. There are two limitations to this. Firstly, panel observation means that some assumption about the status of activity between clinic visits has to be made (e.g. it keeps the status observed at the last clinic visit until the current visit) and, from a causal inference perspective, activity is an internal time dependent covariate so treating it as fixed makes it harder to infer causality.

The second approach seeks to address these problems by jointly modelling activity and joint damage as linked stochastic processes. In this model each of the 28 joints is represented by a three-state model where the first two states represent presence and absence of activity for an undamaged joint, whilst state three represents joint damage. For this model, rather than explicitly fitting a random effects model, a GEE type approach is taken where the parameter point estimates are based on maximizing the likelihood assuming independence of the 28 joint processes, but standard errors are calculated using a sandwich estimator based on patient level clustering. This approach is valid provided the point estimates are consistent, which will be the case if the marginal processes are Markov. For instance if the transition intensities of type {ij} are linked via a random effect u such that but crucially if other transition intensities also have associated random effects, these must be independent of .

The basic model considered is (conditionally) Markov. The authors attempt to relax this assumption by allowing the transition intensities to joint damage from the non-active state to additionally depend on a time dependent covariate indicating whether activity has ever been observed. Ideally the intensity should depend on whether the patient has ever been active rather than whether they have been observed to be active. This could be incorporated very straightforwardly by modelling the process as a 4 state model where state 1 represents never active, no damage, state 2 represents active, no damage, state 3 represents not currently active (but have been in the past) and state 4 represents damage. Clearly, we cannot always determine whether the current occupied state is 1 or 3 so the observed data would come from a simple hidden Markov model.

On a practical level, the authors fit the likelihood for the gamma random effects model by using the integrate function in R, for each patient's likelihood contribution for each likelihood evaluation in the (derivative-free) optimization. Presumably this is very slow. Using a simpler more explicit quadrature formulation would likely improve speed (at a very modest cost to accuracy) because the contributions for all patients for each quadrature point could be calculated at the same time and the overall likelihood contributions could then be calculated in the same operation. Alternatively, the likelihood for a single shared Gamma frailty is in fact available in closed form. This follows from arguments extending the results in the tracking model of Satten (Biometrics, 1999). Essentially, we can write every term in the conditional likelihood as some weighted sum of exponentials:

the product of these terms is then still in this form:

Calculating the marginal likelihood then just involves computing a series of Laplace transforms.
Provided the models are progressive, keeping track of the separate terms is feasible, although the presence of time dependent covariates makes this approach less attractive.

Parametric inference for time-to-failure in multi-state semi-Markov models: A comparison of marginal and process approaches

Yang Yang and Vijayan Nair have a new paper in Canadian Journal of Statistics as part of the Kalbfleisch and Lawless special issue. The paper considers inference for progressive semi-Markov processes for complete, right-censored and interval-censored data for parametric models with Gamma or inverse Gamma sojourn distributions. The advantage of considering Gamma or inverse Gamma distributions is that their convolutions have a closed form meaning the overall failure (absorption) time distribution is available in closed form. Inference using only the time-to-failure data is therefore relatively straightforward. The main focus of the paper is investigating the loss in efficiency of only using time-to-failure data when additional data on intermediate states (either through continuous observation or panel data) are available. The authors show that the loss in efficiency can be rather substantial, particularly if interest lies in making predictions about time to failure given an existing process history rather than just mean or median survival. It is therefore suggested that information on intermediate states should be incorporated where possible despite this presenting computationally difficulties in the case of panel data.

Discrete-time semi-Markov modeling of human papillomavirus persistence

Mitchell, Hudgens, King, Cu-Uvin, Lo, Rompalo, Sobel and Smith have a new paper in Statistics in Medicine. This considers a non-parametric estimator for a 2-state discrete-time semi-Markov process. Following Kang and Lagakos who assume one of the two states (e.g. state 0) is Markov, the process can be characterised by

where is the probability of making a transition from 1 to 0 given i time units spent in state 1,
is the maximum length of state sequence observable in the data and
Extensions to the model, allowing to depend on time in state and to allow an additional disease-free state corresponding to disease-free with no past disease, are also proposed. Missing observations can be dealt with by summing over all possible observed states at the missing times. Estimation of is by maximum likelihood. An acknowledged limitation is the inability to cope with the case of unknown initiation times if either the sojourn distribution of state 0 is non-geometric or observation can start in the disease state.

Of particular interest to the authors is an estimate of 'persistence' of the disease state. This is defined as spending j time units in the disease state, counting single disease free (negative) observations surrounded by positive observations as time spent in the disease. The probability of persistence is just a function of the transition probabilities and so readily estimable.

The authors claim that their discrete time model doesn't not require a "guarantee time" unlike Kang and Lagakos. This is obviously ridiculous, the discrete time model requires a guarantee time of 1 time unit for all transitions! While adopting a discrete time model simplifies the problem of inference to something quite trivial, one has to question how realistic it is to model something that is clearly a continuous time process as discrete time. Bachetti et al's more general approach is along similar lines. Similarly, while the estimation is nominally non-parametric, the discrete time assumption is in many respects more severe than, say, constraining sojourn distributions to be Weibull distributed.

The clinical definition of persistence which makes the assumption that a negative observation between two positive observations counts as a positive is easily accommodated for via the discrete time model. However, a more satisfactory approach would be to adopt a more formal definition, based in continuous time, e.g. persistence if disease free period is less than say 6 months. This would have parallels with the approach taken by Mandel (2010) in defining a hitting time in terms of having a sojourn of more than some length in the disease state. Farewell and Su also dealt with a similar problem but their approach seems to be best avoided.

Non-homogeneous Markov process models with informative observations with an application to Alzheimer's disease

Baojiang Chen and Xiao-Hua Zhou have a new paper in Biometrical Journal. The methodological development of the paper is to extend the methods of Chen, Yi and Cook (Stat Med, 2010) to the case of a non-homogeneous Markov model using the time transformation model of Hubbard et al (Biometrics, 2008) rather than the piecewise constant intensities used in Chen, Yi and Cook. The authors claim that the time transformation model is more appealing than piecewise constant intensities because it requires fewer parameters. However, this parsimony is at the cost of flexibility, as the time transformation model assumes the same temporal trend for all intensities.

The assumption of non-informative examination times is an ever present spectre for multi-state models from panel data. Chen et al's method provides some methods when a complete set of planned examination times is known and it is simply the case that some examination times are missed, meaning the problem can be dealt within the Rubin framework of MAR/MNAR. A more general situation would be where the multi-state process and the process that generates the examination times are dependent. Here the only option seems to be to jointly model the two processes explicitly. A starting model might be one where the intensities of the counting process generating the examination times and the multi-state model are linked through a joint frailty, e.g. something analogous to models for joint modelling of longitudinal and (informative) drop-out (survival) processes.

Progression of liver cirrhosis to HCC: an application of hidden Markov model.

Nicola Bartolomeo, Paolo Trerotoli and Gabriella Serio have a new paper in BMC Medical Research Methodology. This applies a three state hidden Markov model to data on the progression of liver cirrhosis to Hepatocellular carcinoma. A time homogeneous continuous time progressive model is fitted with death as the absorbing state. Covariate effects are included via a proportional intensities model.
Schoenfeld residuals, which are appropriate for right censored data, are applied here as a test of proportionality. It isn't made clear precisely how this is done here. If time of death is known exactly then a Schoenfeld type residual could be defined for the times of death replacing the standard formulation



with



where If the times of death are interval censored then this approach is inappropriate.

On a more trivial level the matrix of misclassification probabilities is missing a 1 in the third row corresponding to the absorbing state.

Flexible Nonhomogeneous Markov models for panel observed data

Andrew Titman has a new paper in Biometrics. This develops an approach to fitting non-homogeneous Markov models to panel data by use of direct numerical solution of the Kolmogorov forward equations. Existing methods to non-homogeneity have concentrated on special cases where the forward equations have matrix analytic solutions (i.e. piecewise constant intensities or time transformation models), although numerical solutions have been used in Bayesian analyses mainly to accommodate use in WinBUGS (see e.g. Welton and Ades or Pan and Chen). The approach is clearly somewhat more computationally intensive than matrix analytic methods. However, a couple of computational tricks are used to improve the situation. In particular a Fisher scoring algorithm is maintained by solving an extended system of ODEs incorporating the first derivatives of the transition probabilities with respect to the model parameters. The real point at which the method struggles is when there are continuous covariates because a separate ODE must be integrated for each covariate value in the data. For large datasets an exact approach becomes untenable. An approximate method is proposed in these situations where a clustering algorithm is used to reduce the number of unique covariate patterns and then each patient is assumed to have covariate pattern equal to the mean value within their cluster. This approach gives pretty close results to the exact method for 10 clusters and is even better for 50 or 100 covariate values where the method is often still practical.

B-spline functions for the transition intensities based on a known set of knot points are proposed. This gives a model which can viewed both as a generalization of the flexible time transformation approach of Hubbard et al, and also as a smooth alternative to piecewise constant intensities. One downside of the great flexibility is that it is quickly possible to run into models with identifiability problems and singular Fisher information matrices. Titman proposes to limit the spline to a maximum time and assume constant intensities beyond this range. Similarly, he suggests there often wont be enough data to allow inhomogeneity on all transition intensities and the method is thus most useful when one or two intensities are of most interest. For the CAV data analyzed, disease onset is of most importance and the method performs better at picking up the increasing hazard than a time transformation model that requires the inhomogeneity to be proportional between intensities.

As part of the supplementary materials some fairly general R code, working in conjunction with the R package deSolve for the ODE solver, is provided. This is flexible in allowing user defined forms for the generator matrix, but there is no easy interface for someone to specify they want a 5-state model with a certain set of allowable intensities and inhomogeneity in a particular set of states, in the same way as for piecewise constant intensities in msm say.

A comparison of non-homogeneous Markov regression models with application to Alzheimer's disease progression

Rebecca Hubbard and Andrew Zhou have a new paper in Journal of Applied Statistics. This considers panel data relating to the progression of Alzheimer's disease using non-homogeneous Markov models. The time transformation method proposed by Hubbard et al (Biometrics, 2008) is extended to allow for (fixed) covariates. In the time transformation model the generator matrix is for some increasing function , which in Hubbard's method can be estimated non-parametrically (or at least flexibly using kernel weights centered on pre-specified time points). Covariates are incorporated as proportional intensities on the intensities in the homogeneous time generator .

A simulation study is presented which aims to compare the robustness of piecewise constant intensity models and time transformation models for estimating covariate effects on intensities. It is no secret that in proportional hazard type models, estimates are reasonably robust to misspecification of the shape of the baseline hazard (in contrast to relative lack of robustness to misspecification of the model for the effect of the covariate on the hazard). In general therefore, each model performs fairly well when the other is true. The simulation puts the piecewise intensities model at a bit of a disadvantage since it requires Q to be estimated separately for each time period (four extra parameters), whereas the time transformation model only requires 1 extra parameter over the homogeneous model. It might have been a fairer comparison if a time transformation model of similar complexity was used. The conclusion that time transformation models are more robust for small sample sizes is therefore not particularly convincing.

The argument against piecewise constant intensities in general is a bit weak as it revolves around the notion that one must estimate a separate generator matrix for each time period leading to many extra parameters. The obvious argument for piecewise constant intensities is that we can look at a particular intensity of interest and let that be time varying while leaving the others constant. In contrast the time transformation method requires the non-homogeneity to be the same for all intensities. Obviously the smoothness of the intensities in the time transformation model is an advantage though.

In the Alzheimer's example a 4-state model is fitted where patients can be in states Normal, Mildly Cognitively Impaired, Alzheimer's or Death. Since the observed data include patients who make transitions back from Alzheimer's to MCI and from MCI to Normal, the Markov models assume backward transitions are possible. In the conclusion it is noted that there is the possibility of misclassification. However, the possible remedy of a hidden Markov model is not mentioned.

Multi-State Models for Panel Data: The msm Package for R

The final paper in the special issue in Journal of Statistical Software is by Chris Jackson and is about his package msm. msm has been around for many (over 8) years and has steadily built up functionality over time. As noted by Hein Putter in his introduction, msm is different from the other packages featured in the issue in that it deals with panel observed/interval censored data and concentrates on parametric models. In particular, hidden Markov models with both discrete and continuous responses may be fitted.

The paper covers much old ground, the early part repeating similar themes from Jackson & Sharples 2002 (Statistics in Medicine) and Jackson et al 2003 (JRSS D), the section on model diagnostics closely follows Titman and Sharples 2010 (SMMR). One error in msm is in the implementation of prevalence counts/plots. It is made clear by Gentleman et al (1994, Stat Med) and Titman and Sharples that the denominator for the counts at time t is the number "under observation". In particular, subjects who reach the absorbing state should stay under observation only until the time at which they would have been censored. msm assumes subjects who reach the absorbing state remain under observation indefinitely which leads to overestimation of the empirical prevalence in the absorbing state(s). This can be seen clearly in the bottom-right panel of figure 4 of the paper that is suggesting spurious lack of fit. Patients either need to be removed from observation at a known administrative censoring time, or else the censoring distribution needs to be empirically estimated to allow time-dependent weighting of dead patients.

In Section 6 Jackson discusses model extensions which are generally not available in the package, but seems to suggest that to maintain generality these extensions, such as a wider range of time inhomogeneous models, random effects models and semi-Markov models, will not be incorporated into msm.

In general msm is a very good package and lots of effort (e.g. C programming, use of first derivatives, direct coding of transition probabilities for more basic model structures) has been put in to provide a fast performance. However, some improvement in computation speed is no doubt possible since msm uses the BFGS algorithm in optim to fit models rather than applying the well-known Fisher scoring algorithm.

Parameterization of treatment effects for meta-analysis in multi-state Markov models

Malcolm Price, Nicky Welton and Tony Ades have a new paper in Statistics in Medicine. This considers how to parameterize Markov models of clinical trials, so that meta-analysis can be performed. They consider data that is panel observed (at regular time intervals) but aggregated into the number of transitions of each type among all patients between two consecutive observations. As a result, it is necessary to assume a time homogeneous Markov model, since there is no information to consider time inhomogeneity or patient inhomogeneity. These types of cost-effectiveness trials are usually analyzed using a discrete time framework (Markov transition model), with the assumption that only one transition is allowed between observations. Price et al advocate using a continuous time method. The main advantage of this is that fewer parameters are required e.g. rare jumps to distant states can be explained by the presence of multiple jumps rather than having to include the transition probability as an extra parameter in the model.

Price et al consider asthma trial data where they are interested in combining data from 5 distinct two-arm RCTs, there is some overlap between the treatments tested so evidence networks comparing treatments can be constructed.

The main weakness of the approach is the reliance on the DIC. Different models for the set of non-zero transition intensities are considered using the DIC, with the models allowing different intensities for each trial treatment arm. While later in the paper there are benefits of adopting the Bayesian paradigm, here it doesn't seem useful. For these fixed effects models and given the flat priors chosen, DIC should (asymptotically) be the same as AIC. However, AIC is dubious here also because of the aspect of testing on the boundary of the parameter space and is likely to put too much favour on simpler models in this context. Likelihood ratio tests based on mixtures of Chi-squared distributions (Self and Liang, 1987) could be applied as in Gentleman et al 1994. However, judging from the example data given for treatments A and D, the real issue is lack of data, e.g. there are very few transitions to and from state X.

maxLik: A package for maximum likelihood estimation

Arne Henningsen and Ott Toomet have a new paper in Computational Statistics. The paper isn't directly related to multi-state modelling, but rather is on their package, maxLik, for general maximum likelihood estimation. Primarily their package is a wrapper for existing optimization packages in R such as optim and nlm. However, they do in addition provide an implementation of the BHHH (Berndt, Hall, Hall and Hausman) algorithm. This is a quasi-Newton algorithm in a similar vein to Fisher scoring, but rather than use the Expected Fisher information, it uses the mean of the outer product of the scores of each observation. Like the BFGS algorithm, a line search is performed to find the step length at each iteration.

For panel observed Markov multi-state models the Fisher scoring algorithm proposed by Kalbfleisch and Lawless (1985, JASA) and generalised by Gentleman et al (1994, Stat Med) is superior to BHHH. However, for data with a mixture of panel observed observations and exact times of absorption (e.g. death), the Fisher scoring algorithm cannot be applied. Here BHHH seems to perform significantly better than the BFGS algorithm supplying first derivatives.

A multistate model for events defined by prolonged observation

Vern Farewell and Li Su have a new paper in Biostatistics. This models remission in psoriatic arthritis, based on panel data assessing joints at clinic visits. Existing models use a two-state model. However, a spell in remission should last a discernible amount of time. Rather than specify an artificial length of time (e.g. 6 months) e.g. a guarantee time in a state, Farewell and Su adopt a model with two states referring to remission: they refer to these as "early stage remission" and "established remission". It is assumed an individual must progress through both early and established remission before returning to active disease. A subject is observed to be in early stage remission if they have no active joints at a visit not preceded by at least 2 other zero count visits. They are in established remission if there is a zero count and at least 2 previous zero counts. State misclassification is allowed in the model through misclassification of the active disease count, i.e. patients may have 0 active joints without being in remission. It is assumed that misclassification to the early stage remission is possible but not to the established remission stage. In the example the misclassification probability is also allowed to depend on whether the previous observed count was zero or not.

The basic problem with the method is that having the states defined by the pattern of previous observations means that it is essentially impossible for the observed data (in terms of the three-state model) to come from the claimed Markov model: In the observed data, an established remission stage must be preceded by two early stage remission observations. Yet the actual Markov model allows the passage time from active disease to established remission to be arbitrarily close to zero (e.g. just the sum of two independent exponential - or perhaps piecewise exponential - distributions). As a result it is not clear how to interpret the resulting transition intensity estimates since the estimated process will not reproduce the original data.

Misclassification is effectively dealt with twice in the model. Firstly in an ad hoc way through rules on what early and established remission are. Then by allowing these observed states to have classification error over some true states. But in fitting a hidden Markov model the misclassification is assumed independent conditional on the underlying state. There is then an inherent contradiction because on the one hand the model says P(Observed zero | Active disease) >0, but at the same time P(Observed zero and two previous zero | Active disease) =0.

An approach using a guarantee time (e.g. Kang and Lagakos (2007)) or perhaps an Erlang distribution through latent states would be far more satisfactory even if it might require "special software". Potentially the guarantee time could be dependent on covariates.

Analysis of interval-censored disease progression data via multi-state models under a nonignorable inspection process

Chen, Yi and Cook have a new paper in Statistics in Medicine. This covers similar ground to the paper of Sweeting et al as it is also concerned with a nonignorable inspection process for a panel observed progressive Markov multi-state model. Again there are a set of known planned observation times, but unlike Sweeting et al there is no auxiliary data. A MNAR model where the missingness probability depends on the true state at the potential observation time is used. This is similar to the MNAR model in Sweeting et al, except that here piecewise constant transition intensities are used rather than a time homogeneous model. The model is fitted using an EM algorithm.

Both this study and Sweeting et al's depend on there being a known set of planned examination times. In the majority of cases, only the actual examination times are known and the observation process is not fully understood. As Chen et al note, informative observation times could be dealt with by determining a stochastic model for the observation process, but would require assumptions about the process.

Multistate Markov models for disease progression in the presence of informative examination times

Sweeting, Farewell and De Angelis have a new paper in Statistics in Medicine. This deals with the problem of a panel observed disease process where the examination times are generated by a non-ignorable mechanism. This is in general a very difficult problem. The authors consider a special case where, while the disease process is only observed at informative examination times, an auxiliary variable is able to be observed at a full set of planned (or ignorable) times. They consider a model where the disease process is missing at random (MAR) conditional on the values of the auxiliary variable. The likelihood then becomes of the form of a (partially) hidden Markov model. An alternative missing not at random model (MNAR) with missingness dependent on the current state is also considered. However the comparison between the MAR and MNAR models is somewhat unfair. In the simulations and the application, the disease process is only observed at a minority of planned examination times. The auxiliary variable has a fairly strong correlation with the disease process, meaning significant information about the process can be obtained through using the auxiliary variable in a hidden Markov model. Indeed, even without the suspicion of informative missingness the HMM would be worth using (if the strong assumptions about the relationship with the auxiliary variable could be reliably assumed) for improved efficiency. However, the MNAR model does not use the auxiliary variable at all. A better model to compare would be a hybrid of the two where a partially HMM is used in conjunction with the logistic model for the MNAR (although this wouldn't necessarily be consistent under AD-MAR unless the logistic model was correctly specified). This does not seem like a difficult extension.

A further shortcoming of the paper is the lack of any suggestions for model checking. Measurements from the auxiliary variable are assumed to be Normally distributed and independent conditional on the underlying disease states at the examination times. This seems fairly unrealistic and should at least be supported by the data.

Non-Markov Multistate Modeling Using Time-Varying Covariates

Bacchetti et al have a new paper in The International Journal of Biostatistics. This considers modelling progression of liver fibrosis due to hepatitis C following liver transplant using a 5-state progressive multi-state model. The data are panel observed at irregular time points, so it would be most natural to model the data in continuous time. In addition the observed states are subject to classification error. To avoid having to making either Markov or time homogeneity assumptions, the authors adopt a discrete time assumption: assuming 4 time periods per year. They then model the transition probabilities as linear on the log-odds scale, depending on covariates such as medical center, donor age, year of transplant as well as log time since entry into the current state (relaxing the Markov assumption). Potentially, time since transplant could also be included in the model (for non-homogeneity).

The main challenge in fitting the models is to enumerate all possible complete "paths" of true states at the discrete time points that could result in the observed data. Obviously an iterative algorithm is needed for this. The computational complexity will depend on the complexity of the state transition matrix, the misclassification probability matrix and the degree of discretization.

The main drawback of approximating a continuous time process by a discrete-time process is the restriction that only 1 transition may occur between time points. While this can be acceptable for progressive models such as the one considered here, it may be more problematic when backward transitions are allowable.

The authors have developed an R package called mspath, designed to complement the existing package for continuous-time Markov and hidden Markov model msm, to allow non-Markov models through the discrete-time approximation.