Effect of an event occurring over time and confounded by health status: estimation and interpretation. A study based on survival data simulations with application on breast cancer

Alexia Savignoni, David Hajage, Pascale Tubert-Bitter and Yann De Ryckea have a new paper in Statistics in Medicine. This considers developing illness-death type models to investigate the effect of pregnancy on the risk of recurrence of cancer amongst breast cancer patients. The authors give a fairly clear account of different potential models with particular reference to the hazard ratio The simplest model to consider is a Cox model with a single time dependent covariate representing pregnancy, here . This can be extended by assuming non-proportional hazards which effectively makes the effect time dependent i.e. . Alternatively, an unrestricted Cox-Markov model could be fitted with separate covariate effects and non-parametric hazards from each pregnancy state, yielding: This model can be restricted by allowing a shared baseline hazard for giving either under a Cox model with a fixed effect or for a time dependent effect.

If we were only interested in and any of these models seems feasible, there doesn't actual seem that much point in formulating the model as an illness-death model. Note that the transition rate does not feature in any of the above equations but would be estimated in the illness-death model. The above models can be fitted by a Cox model with a time dependent covariate (representing pregnancy) that has an interaction with the time fixed covariates. The real power of a multi-state model approach would only become apparent if we were interested in the overall survival for different covariates, treating pregnancy as a random event.

The time dependent effects are represented simply via a piecewise constant time indicator in the model. The authors do acknowledge that a spline model would have been better. The other issue that could have been considered is whether the effect of pregnancy depends on time since initiation of pregnancy (i.e. a semi-Markov effect). An issue in their data example is that pregnancy is only determined via a successful birth meaning there may be some truncation in the sample (through births prevented due to relapse/death).

Bayesian analysis of a disability model for lung cancer survival

Armero, Cabras, Castellanos, Perra, Quirós, Oruezábal and Sánchez-Rubio have a new paper in Statistical Methods in Medical Research. This develops a Bayesian three-state illness-death type model to the progression and survival of lung cancer patients. The data considered are assumed to be complete up to right censoring (in reality there may be some interval censoring but the authors argue the patients can be considered `quasi-continuously' followed. A Weibull semi-Markov model is assumed for the transition intensities and covariates are accommodated via an accelerated failure time model (it's worth noting that for a Weibull distribution the proportional hazard and accelerated failure time models are equivalent up to reparameterization). A feature of the dataset used is the rather small sample size (35 patients) which is perhaps the strongest reason for taking a parametric and Bayesian approach in this case.

Investigating hospital heterogeneity with a multi-state frailty model

Benoit Liquet, Jean Francois Timsit and Virginie Rondeau have a new paper in BMC Medical Research Methodology. They consider data on the evolution of patient's status for patients in intensive care units. The data originate from 16 distinct ICUs and there is therefore inherent clustering in the data. A progressive four-state model, with admission as state 0, ventilator-associated pneumonia infection (VAP) as state 1, death as state 2 and discharge as state 3 is considered. To account for the clustering, shared Gamma frailty terms are included in the intensities for particular transitions. To maintain simplicity of the method the authors consider cases where either each transition intensity has a ICU related frailty that is assumed independent of frailties for other intensities, or where there are individual level frailties that act in common across intensities. In each, there is only one level of clustering (i.e. either independent frailties on each intensity common to each centre or a subject specific intensity affecting multiple intensities). In the first case, a non-homogeneous Markov or semi-Markov allows separate models to be fitted for each transition using methods applicable to univariate survival analysis. In the second case, multiple intensities can be fitted in a single survival model by specifying separate "strata" for each intensity. In either case the presence of only one level of clustering makes estimation, at least if a Gamma frailty is assumed, relatively straightforward. The authors use the R package frailtypack which Rondeau maintains. This allows fully parametric models (via either Weibull or piecewise constant intensities) or semi-parametric models using spline intensities fitted using penalized likelihood. The emphasis on simple(r) approaches, for which existing software exist, is a good one. But some mention of the existing capacity of the more general survival package to fit Cox models with gamma frailties by using frailty() in the formula.

The wider issue of what to do in the case of nested frailties, i.e. where there are multiple layers of clustering, is an interesting one, but is not discussed which is surprising given that frailtypack has some capability for fitting such models. It is also questionable whether the assumption of independent frailties across intensities is a realistic one. To some extent this doesn't matter because if the data are Markov or semi-Markov conditional on the frailties, then even if the frailties are dependent assuming independence should produce unbiased estimates of the marginal frailty distribution (provided it is Gamma). Similarly estimates of the individual frailties should be reasonably unbiased (empirical Bayes estimates aren't fully unbiased anyway) but will be inefficient if frailty terms across intensities are actually correlated. It would be relatively straightforward to look at the correlation in the empirical Bayes estimates of the frailties associated with the intensities under the independence model as a semi-informal diagnostic. The difficult part would be fitting the multivariate frailty model if the independence model appeared inadequate.

Estimating Discrete Markov Models From Various Incomplete Data Schemes

Alberto Pasanisi, Shuai Fu and Nicolas Bousquet have a new paper in Computational Statistics & Data Analysis. This considers approaches to Bayesian inference for time-homogeneous discrete-time Markov models under incomplete observation. Firstly, they consider where there are missing observations in a sequence of states (considering different missingness assumptions). Secondly, they consider the case of aggregate data where all that is known is the number of subjects in each state at each time. A Bayesian approach is adopted throughout, which the authors claim is the most convenient in this situation.

The part involving incomplete data follows similar ground to Deltour et al (Biometrics, 1999). The problem only becomes non-trivial if the missingness mechanism is non-ignorable.

The treatment of aggregate data is incorrect as the authors state the likelihood as being the product of independent multinomial random variables with probabilities corresponding to the probability of being in a state at time t given the initial state distribution at time 0. As a result they claim the likelihood is proportional to the case of current status data where each subject or unit is only observed once. The reason that likelihood-based inference for aggregate data is so difficult is that we observe all units multiple times but don't know number or nature of the transitions that occurred. Hence, the full likelihood would require summing over all possible transitions consistent with the aggregate counts. Kalbfleisch and Lawless (Canadian Journal of Statistics, 1984) derived the mean and covariance of the aggregated counts across times to establish a least-squares estimation procedure. Pasanisi et al's procedure is only relevant when the data consist of a series of independent cross-sectional surveys at different time points all assumed to come from different units. An MCMC or simulation based approach would be necessary to compute the exact likelihood or posterior distribution in the true aggregate data case, which the authors did not pursue. However, the gain in efficiency compared to the least-squares approach is probably not worth the trouble except for very small counts. Crowder and Stephens (2011) pursued an approach based on matching the coefficients of the probability generating function of the aggregate counts.

Bayesian analysis of multistate event history data: beta-Dirichlet process prior

Yongdai Kim, Lancelot James and Rafael Weissbach have a new paper in Biometrika. This develops a conjugate prior process suitable for non-parametric and semi-parametric Bayesian modelling of right-censored multi-state Markov data. The model is parametrised in terms of the sum of the intensities out of each state

and instantaneous transition probabilities

A possible choice for a prior process is a Dirichlet distribution but this is not independent in the limit of a continuous time process. Instead the authors propose a new beta-Dirichlet process consisting of a beta distributed part which determines the increment in (between 0 and 1) and a Dirichlet part determining the instantaneous transition probabilities for each particular transition. The authors prove this prior process is conjugate in the continuous limit.

A semi-parametric regression model is proposed, which the authors term as a semi-proportional intensities model. This consists of a proportional intensities model for the all-cause hazard of exiting state h and a multinomial type model for the instantaneous transition probabilities out of state h and bears some resemblance to the vertical modeling parametrization for competing risks regression.

In an aside the authors claim that interval censoring can easily be dealt with by treating the unknown transition time as missing data that can be accounted for in the Gibbs sampling. This only works under the assumption that only one transition can have occurred between examination times. While other authors have made this assumption (e.g. Foucher et al 2007) it is dubious to say the least and likely to result in biased estimates. Similarly, the authors claim right-censoring can be dealt with by treating a censoring event as an additional state. While this will obviously allow the observed process to be modelled, it is not clear how this approach would allow the underlying process (without censoring) is estimated?

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.

Relative survival multistate Markov model

Ella Huszti, Michal Abrahamowicz, Ahmadou Alioum, Christine Binquet and Catherine Quantin have a new paper in Statistics in Medicine. This develops a illness-death Markov model with piecewise constant intensities in order to fit a relative survival model. Such models seek to compare the mortality rates from disease states with those in the general population, so that the hazard from state r to absorbing state D are given by

The population hazard is generally assumed known and taken from external sources like life tables. Transition intensities between transient states in the Markov model are not subject to any such restrictions. One can think this as being a model where there is an unobservable cause of death state "death from natural causes" which has the same transition intensity regardless of the current transient state.

The paper has quite a lot of simulation results, most of which seem unnecessary. They simulate data from the relative survival model and show, unsurprisingly, that a misspecified model that assumes proportionality with respect to the whole hazard (rather than the excess hazard) is biased. They also compare the results with Cox regression models (on each transition intensity) and Lunn-McNeil competing risks model (i.e. assuming a Cox model assuming common baseline for the competing risks).
The data are a mixture of clinic visits that yield interval censored times of events such as recurrence, but times of death are known exactly. Presumably for the Cox models a further approximation of interval-censored to right-censored data is made.

Combined survival analysis of cardiac patients by a Cox PH model and a Markov chain

Michal Shauly, Gad Rabinowitz, Harel Gilutz and Yisrael Parmet have a new paper in Lifetime Data Analysis. This considers methods for modelling the effect of a mixture of time dependent and time constant covariates on overall survival, with the complication that the time dependent covariates are only observed at a discrete set of time points. They propose to firstly fit a Cox proportional hazard model assuming all covariates are fixed at their baseline values. The main modelling approach is to assume a discrete-time homogeneous Markov model with states corresponding to the combinations of the time dependent covariates (which are categorical or need to be categorized) and death. Transitions between all covariates states are assumed to be possible between each time point. The Cox model is used to determine which of the constant covariates should be considered in the Markov model. For this the authors propose to again categorize the covariates and consider a separate Markov model for each level of the covariates. Having obtained the estimates from the Markov models, it is then possible to calculate expected survival times for patients conditional on their baseline characteristics.

In general the approach proposed is reasonably sensible. However, there is panel data available for the time dependent covariates. It therefore seems possible to fit a time continuous Markov model to the data using methods appropriate for panel data (e.g. Kalbfleisch and Lawless, 1985). This approach has the advantage that the exact time of the death events can still be used.

The authors rely on categorization throughout. While this seems necessary for the time dependent covariates, there seems scope for using multinomial logit models for other covariates. Similarly, by allowing a different mortality probability for each combination of covariates they are effectively fitting covariate models with interactions (i.e. the effect of being in covariate level 2 compared to 1, is different depending on which level(s) of the other time dependent covariate(s) a subject is in). While such interactions may be necessary, it might be better to allow simpler models where only the evolution of the time dependent covariates is kept general. This is another advantage of a continuous time model since covariate effects could remain on the hazard (transition intensity) scale as in the Cox PH model.

Finally, the authors give a partial justification of the use of a time homogeneous Markov model through the Cox PH model having an approximately constant baseline hazard. It should be noted that a time homogeneous Markov model does not imply a constant absorption hazard (unless the model begins in the quasi-stationary distribution). Conversely, while a constant hazard might suggest homogeneity more than non-homogeneity, it is nevertheless possible to construct non-homogeneous processes with constant (or near constant) marginal absorption hazards. The authors do however report a statistic which gives a better justification of homogeneity.

Comparison of prediction models for competing risks with time-dependent covariates

Giuliana Cortese, Thomas Gerds and Per Kragh Andersen have a new paper available as a University of Copenhagen Department of Biostatistics technical report. The paper concerns the development of models for prediction for competing risks in the presence of internal time dependent covariates. This is a follow-up to Cortese and Andersen's 2010 Biometrical Journal paper. Like the previous paper, the authors compare a multi-state modelling approach that explicitly models the progression of the (categorical) time dependent covariate and its effect on the cause specific hazards, and a landmarking approach that sets (arbitrarily chosen) time points and performs separate regressions to estimate the hazards for conditional on the value of the time dependent covariate at time . The authors consider two modelling approaches under landmarking, one based on Cox regression of cause-specific hazards and the other based on Fine-Gray subdistribution hazards.

They compare the predictive ability of the models to predict the outcome by landmark given data up to . This is assessed by using a time dependent Brier score (Gerds & Schumacher, 2006). Rather than use inverse probability weighting, the authors instead use a pseudo-value to estimate outcomes when a subject is lost to follow-up between and . The authors perform the comparison using a bone marrow transplant study where the competing events are relapse and death, and the internal time dependent covariate is the development of Graft versus Host Disease (GvHD). The predictive abilities are estimated via cross-validation involving randomly choosing 2/3 of patients as training data and using the remainder as test data, repeating the process 100 times. For the data considered the three methods performed equally well in terms of prediction error. As might be expected, there was significantly improved predictive ability of these models compared to one that ignored GvHD (i.e. only considered baseline time constant covariates).

As the authors note, there are advantages and disadvantages to both approaches. The multi-state modelling approach requires modelling of the covariate process (e.g. Markov or semi-Markov assumptions and proportional hazard assumptions on the effect of baseline covariates on transition rates through covariate states) and requires a categorical covariate. Landmarking can accommodate continuous covariates but relies on an arbitrary set of landmark times and requires fitting regressions at each landmark. The extra modelling required for the multi-state approach may either be a blessing, in terms of having the potential to give a greater insight into the whole process, or a curse (questions of robustness to incorrect modelling assumptions).

Lie Markov Models

Jeremy Sumner, Jesus Fernandez-Sanchez and Peter Jarvis have a paper recently made available at thehttp://www.blogger.com/img/blank.gif Arxiv. The paper is theoretical in nature but gives an interesting application of group theory to Markov models.

The practical problem addressed is determining the conditions under which a non homogeneous continuous time Markov model can be represented by a "time averaged" homogeneous Markov model, i.e. what constraints are required to ensure a rate matrix exists such that

for rate matrices . This has application for phylogenetic methods, where typically a single rate matrix is fitted to an evolutionary history. However, rates may in fact change over time. The question is then under what conditions could the single rate still be in some sense valid in summarising the time averaged process.

The authors show that the model for Q must be a Lie algebra. They also give the possible model forms for three and four state models under symmetry constraints. Update: This paper has now been published in the Journal of Theoretical Biology.

On inference from Markov chain macro-data using transforms

Martin Crowder and David Stephens have a new paper in Journal of Statistical Planning and Inference. This considers estimation of a discrete-time homogeneous Markov chain from aggregate data (which they term macro-data), ie. where only the overall state counts for N patients are known at time j, while the transition counts are unknown. The likelihood for such data is intractable because it involves summing over the vast number of possible transitions that are consistent with the observed aggregate counts. As a result, inference methods have focused on moment based estimation (see e.g. Kalbfleisch, Lawless and Vollmer, Biometrics 1983) which are reasonably effective in practice.

Crowder and Stephens note that the probability generating function for the observed aggregate counts has a fairly simple form. This motivates an estimation procedure based on trying to match the quantities with their expectations (i.e. the pgfs). An obvious practical issue is the choice of vectors to use to compare the closeness between the pgf and sample quantities. The authors propose choices that avoid computational problems for large sample sizes.

Through a series of simulations, the authors demonstrate that an improvement in efficiency compared to methods based on second moments is possible for small sample sizes (e.g. n <= 25). A comparison is made with the efficiency for micro-data (i.e. where transition counts are known). However, for such small samples computation of the full likelihood for the aggregate data must become a viable option. The practical issues are whether the pgf based approach gives any real improvement in efficiency compared to second moment approaches for say n=100 (the second moment approaches are asymptotically efficient with appropriate weights) or whether the pgf method outperforms (or matches) the full likelihood for small samples sizes (i.e. n<=25) where the full likelihood is calculable. These issues aren't really addressed in the paper.

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.

msSurv: Nonparametric Estimation for Multistate Models

Nicole Ferguson, Guy Brock and Somnath Datta have written an R package for non-parametric estimation in multi-state models. To some extent the package covers similar ground to mstate and etm, the focus being on data continuously observed up to right censoring. Unlike mstate there is no possibility of semi-parametric modelling. The main area of new functionality in msSurv is the ability to estimate state entry and exit time distributions, and the ability to cope with state dependent censoring mechanisms using the methodology of Datta and Satten (Biometrics, 2002). As with mstate, all computations appear to be performed within R itself. Thus if a standard Aalen-Johansen type estimate is required, etm is still the best package to use. For instance, the example simulated right censored data provided in the package takes over 3 minutes to fit using msSurv, compared to just 1.2 seconds in etm. Since, for the more bespoke parts of the package, e.g. robust estimates of state occupancy for non-Markov models or state dependent censoring, bootstrapping is required for confidence intervals, the lack of speed of msSurv is a little disappointing. Update: A paper on the msSurv package has now been published in the Journal of Statistical Software.

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.

Application of multistate models in hospital epidemiology: Advances and challenges

Jan Beyersmannm, Martin Wolkewitz, Arthur Allignol, Nadine Grambauer and Martin Schumacher have a new paper in Biometrical Journal. The Freiburg group have been very active in promoting multi-state modelling to a wider audience in recent years. This paper looks at the advantages adopting a multi-state model approach can bring to modelling processes related to hospital epidemiology.

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.

mstate: An R Package for the Analysis of Competing Risks and Multi-State Models

The sixth paper in the special issue of Journal of Statistical Software is by Liesbeth de Wreede, Marta Fiocco and Hein Putter and is about the R package mstate. A journal article on the package already exists in Computer Methods and Programs in Biomedicine. However, while that paper primarily dealt with theoretical aspects, the current paper is largely a case-study example based on a 6 state model for leukemia patients after bone marrow transplantation.

mstate uses the existing survival package to fit the Cox proportional hazards models. Much of the infrastructure of mstate is in functions for preparing data to be in the correct form. To use mstate for a Cox-Markov or Cox-semi-Markov model, a single coxph() object is required. This results in somewhat messy commands being required because separate strata are required and the same covariate needs to appear multiple times to allow it to have different effects for each transition intensity. For instance the call to coxph in the paper requires 16 lines. While there is perhaps some pedagogical advantage to this complication, in ensuring the user really understands what they are fitting, there is surely scope to allow some automation to this process so that the user only need specify which covariates are required for each transition intensity and this could be passed to coxph behind the scenes.

Finally, as has been noted elsewhere, mstate currently does virtually all calculations within R itself. As a result computation times are sometimes disappointing, especially for models on large datasets.

Using Lexis Objects for Multi-State Models in R

The fourth and fifth papers in the special issue of Journal of Statistical Software concern Lexis objects in the R package Epi and are by Martin Plummer and Bendix Carstensen. The first paper is mainly a general introduction to Lexis objects whereas the second paper looks specifically at their use for multi-state models. Lexis diagrams, which have a long history, provide a way of illustrating data with multiple time scales and have been promoted by Niels Keiding among others.

The Lexis function makes initial exploration of multi-state data easier, by allowing plotting of events on different timescales so that, for instance, the correct time scale(s) for analysis can be chosen. Lexis objects can also be passed to mstate.

The authors consider the problem is determining whether two baseline transition intensities are proportional. A model allowing proportionality is straightforward to implement as it just involves putting the two transition intensities in the same strata of a Cox model and adding an extra indicator variable to denote which transition intensity an event derives from. However, if one is interested in formally testing an assumption of proportionality there is a difficulty because the Cox model factorizes the baseline intensities out. The authors therefore propose full likelihood modelling where the baseline intensities are assumed to be piecewise constant within short time intervals but follow a smooth spline function. Testing proportionality of intensities can then be performed through a standard likelihood ratio test.

Empirical Transition Matrix of Multi-State Models: The etm Package

The third paper in the special issue of Journal of Statistical Software is by Arthir Allignol, Martin Schumacher and Jan Beyersmann and concerns the etm package. etm has been mentioned here before. It's a fairly straightforward package that computes the Aalen-Johansen estimator of the transition probabilities (and their standard errors) for general right-censored and left-truncated Markov models. mstate also has the same functionality (but also allows Cox regression). The niche of etm is that if only the Aalen-Johansen estimator is required it provides a faster implementation than mstate since etm incorporates C code, whereas mstate does most of its computation in R only.