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.

Nonparametric estimation of the cumulative intensities in an interval censored competing risks model

Halina Frydman and Jun Liu have a new paper in Lifetime Data Analysis. This concerns non-parametric estimation for competing risks models under interval censoring. The problem of estimating the cumulative incidence functions (or sub-distribution functions) under interval censoring has been considered by Hudgens et al (2001) and involves an extension of the NPMLE for standard survival data under interval censoring.

The resulting estimates of the cumulative incidence functions are only defined up to increments on intervals. Moreover, the intervals by which the CIFs are defined are not the same for each competing risk. This causes problems if one wants to convert the CIFs into estimates of the cumulative cause-specific hazards. Frydman and Liu propose estimating the cumulative cause-specific hazards by first constraining the NPMLEs of the CIFs to have the same intervals of support (NB: this is just a sub-set of the set of all NPMLEs involving sub-dividing the intervals) and adopting a convention to distribute the increment within the resulting sub-intervals (they assume an equal distribution across sub-interval).

In addition they show that an ad-hoc estimator of the cumulative hazards based on the convention that the support of each interval of the NPMLE for each CIF occurs at its midpoint leads to biased results. They also show that their estimator has standard N^0.5 convergence when the support of the observation time distribution is discrete and finite.

A semi-Markov model for stroke with piecewise-constant hazards in the presence of left, right and interval censoring

Venediktos Kapetanakis, Fiona Matthews and Ardo van den Hout have a new paper in Statistics in Medicine. This develops a progressive three-state illness-death model for strokes with interval-censored data. The proposed model is a time non-homogeneous Markov model. The main approach to computation is to assume a continuous effect of age (time) on transition intensities but to use a piecewise constant approximation to actually fit it. The intensity to death from the stroke state additionally depends on the time since entry into the state (i.e. age at stroke) and since the exact time is typically unknown, it is necessary to numerically integrate over the possible range of transition times (here using Simpson's rule).

The data include subjects whose time to stroke is left-censored because they have already suffered a stroke before the baseline measurement. The authors state that they cannot integrate out the unknown time of stroke because the left-interval (i.e. the last age at which the subject is known to have been healthy) is unknown. They then proceed to propose a seemingly unnecessary ad-hoc EM-type approach based on estimating the stroke age for these individuals, which requires the arbitrary choice of an age at which it can be assumed the subject was stroke free. However, surely if we can assume a , we can just use as the lower limit in the integral for the likelihood?

The real issue seems to be that all subjects are effectively left-truncated at the time of entry into the study (in the sense that they are only sampled due to not having died before their current age). For subjects who are healthy at baseline this left-truncation is accounted for by just integrating the hazards of transition out of state 1 from their age at baseline rather than age 0. For subjects who have already had a stroke things are more complicated because the fact they have survived provides information on the timing of the stroke (e.g. if stroke increases hazard of death, the fact they have survived implies the stroke occurred sooner than one would assume if no information on survival were known). Essentially the correct likelihood is conditional on survival to time and so the unconditional probability of the observed data needs to be divided through by the unconditional probability of survival to time . For instance, in their notation, a subject in state 2 at baseline censored at time should have likelihood contribution: The authors claim that their convoluted EM-approach has "bypassed the problem of left truncation". In reality, they have explicitly corrected for left-truncation (because the expected transition time is conditional on being in state 2 at baseline) but in a way that is seemingly much more computationally demanding than directly computing the left-truncated likelihood would be.

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.

Estimating survival of dental fillings on the basis of interval-censored data and multi-state models

Pierre Joly, Thomas Gerds, Vibeke Qvist, Daniel Commenges and Niels Keiding have a new paper in Statistics in Medicine. This considers the estimation of survival times of dental fillings from interval censored data. A particular feature of the data is that there is inherent clustering in the form of multiple fillings from the same child.

Bizarrely the authors claim "we are not aware of any paper combining multi-state models for clustered data with interval censoring" implying they (and presumably the referees as well) are unaware of both Cook, Yi, Lee and Gladman (Biometrics, 2004) and
Sutradhar and Cook (JRSS C, 2008), the latter having "clustered", "multistate" and "interval-censored" all in the title!

A progressive four-state model is assumed for each filling, with the states consisting of Treatment, Filling failure, endodontic complication and exfoliation (an absorbing state).
For exfoliation, age of child is taken as the time scale meaning that the time of treatment is taken as a left-truncation time. For transitions from treatment to the other two states, time since treatment is taken as the time scale. Weibull transition intensities are assumed. However, monitoring ended once any filling event (filling failure or endodontic complication) had occurred. Because the exact time of an endodontic complication is known if it occurred during the monitoring period, the transition intensity from endodontic complication to exfoliation is not relevant to the likelihood. The authors also argue that it is necessary to assume that the intensity from filling failure to exfoliation and from treatment to exfoliation is the same, due to never being able to observe a filling failure to exfoliation transition. Strictly speaking, under the assumption of non-informative observation times, there should be some information in the data to estimate something about the separate intensities based on the proportion of cases where a filling was observed compared to the proportion where exfoliation occurred without observing a filling. Indeed in the research report by Frydman et al (2008), using a subset of the data in the current paper and a three-state verison of the model, a discrete-time NPMLE for the intensities was developed.

Random effects are incorporated into the model in a hierarchical way, with a dentist level random effect that affects the intensity to filling failure or endodontic complication and correlated child level random effects determining the correlation to time to failures (thus affecting the transitions to filling failure and endodontic complication) and time to exfoliation (thus affecting only the exfoliation transition intensity). The random effects are taken to be multivariate Normal with a log-additive effect on the intensities. Calculating the likelihood requires numerical integration, which here is achieved via Gauss-Hermite quadrature. 30 quadrature points were used - this seems a rather small number for a multi-dimensional integral. Cook et al (2004) avoided attempting to get a strict approximation to the multivariate Normal by formally restricting the random effects to have a discrete distribution. Sutradhar and Cook (2008) used an MCEM in order to apply a continuous random effects distribution. That approach is likely to be more computationally intensive that Gauss-Hermite quadrature (on 30pts) but more accurate. The recent suggestion by Putter and van Houwelingen to use a simple two-component mixture frailty could also be adaptable for this situation.

Two filling types, amalgam and glass ionomer are compared in terms of probability of surviving without complication, with amalgam performing somewhat better.

A dynamic model for the risk of bladder cancer progression

Núria Porta, M. Luz Calle, Núria Malats and Guadalupe Gómez have a new paper in Statistics in Medicine. This develops a model for progression of bladder cancer with particular emphasis on predicting future risk given events up to a certain point in time.

In many ways the paper is taking a similar approach to Cortese and Andersen in explicitly modelling a time dependent covariate (here recurrence) in order to obtain predictions.
They fit a semi-parametric Cox Markov multi-state model to the data and define a prediction process


where is the time of the second event, is the type of the second event where P denotes progression and represents the history of the process up to time t. Analogously to outcome measures like the cumulative incidence functions, this predictive process is a function of the transition intensities. They also consider time dependent ROC curves to assess the improvement in classification accuracy that can be achieved by taking into account past history in addition to baseline characteristics.

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.

Impact of delayed diagnosis time in estimating progression rates to hepatitis C virus-related cirrhosis and death

Bo Fu, Wenbin Wang and Xin Shi have a new paper in Statistical Methods in Medical Research. The paper is interested in the estimation of rates of progression from hepatitis C infection to cirrhosis and from cirrhosis to death. The primary complication in the model is that the time of cirrhosis development is interval censored, since it is only diagnosed at clinic examination times.

The model consists of a simple unidirectional 3-state model. A parametric Weibull approach is taken, where the time to cirrhosis from infection has a standard Weibull distribution and the time from cirrhosis to death has a Weibull distribution where the scale parameter is a function of the time between infection and cirrhosis.

Unsurprisingly, they show that a full likelihood approach which takes into account the interval censoring is much less biased than an approach that ignores the problem of interval censoring.

A criticism of the model, which the authors acknowledge, is that the possibility of death before cirrhosis is not accommodated. Potentially, even under their current model there might also be cases where death occurs after cirrhosis but before diagnosis of cirrhosis - a scenario which isn't accommodated in the three cases listed in their likelihood development. The rate of increase of the hazard of cirrhosis from infection is observed to increase after 30 years since infection. It is not clear whether this is a real effect, whether it is due to a small number of people at risk or whether it is an artifact of not accommodating the possibility of death without cirrhosis. Given the long follow-up time in the study it might have been more sensible to consider age rather than either time since cirrhosis or time to infection as the primary timescale, at least for time to death.

Finally, whilst the approach might be novel for HPV patients, the 3-state approach is clearly widely used in other contexts, most specifically in HIV/AIDS studies. For instance the penalised likelihood approach of Joly and Commenges would be highly applicable.

Joint modelling of longitudinal outcome and interval-censored competing risk dropout in a schizophrenia clinical trial

Ralitza Gueorguieva, Robert Rosenheck and Haiqun Lin have a new paper in JRSS A. The paper concerns the joint modelling of a longitudinal outcome and an interval censored competing risks outcome that explains drop-out. As is common with these joint longitudinal and survival types of models the two processes are linked via a normally distributed vector of random effects. The novelty of the paper is in the survival part is a competing risks process and the event time is interval censored. The authors adopt a parametric model for the competing risks, using the family of distributions proposed by Sparling et al (Biostatistics, 2006). This makes inference somewhat more straightforward than it would be if a non-parametric baseline cause-specific hazards were used. As recently noted, parametric treatment of competing risks data is surprisingly rare. One problem faced by the authors is that the hazard family of Sparling, while allowing closed form expressions for interval censored univariate survival data, do not result in closed form expressions for interval censored competing risks data (except in special cases). Instead a numerical integral has to be competed. The presence of the overall random effects would mean the likelihood requires nested integration. To avoid this problem the authors adopt an approximation to the true likelihood for competing risks data. If a patient is known to have had a failure of type j in the interval [t0,t1] the authors assume that the patient is censored of all risks except risk j at time t0. It is clear that this approximation will lead to systematic bias as the time at risk from each failure type will be underestimated so the hazards will tend to be overestimated. The amount of bias will depend on the typical length of the intervals [t0,t1].

For the CATIE data example the proposed approximation is probably not an issue. The drop out (competing risks) part of the model is not the primary focus of the inference, and it is really the relative hazards of different types of drop out rather than their absolute values that is important in determining the trajectories of the longitudinal measure without drop out. For instance the estimates for simulated data of a similar type are close to unbiased.
However in extreme cases like current status competing risks data the approximation will do extremely badly.

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.

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.

Parametric Estimation of Cumulative Incidence Functions for Interval Censored Competing Risks Data

Peter Nyangweso[1], Michael Hudgens and Jason Fine have a new paper currently available as a University of North Carolina at Chapel Hill Department of Biostatistics Technical Report. This considers parametric modelling of cumulative incidence functions in the case of interval censored competing risks data. It is somewhat curious that while there has been quite a lot of work on interval-censored competing risks by modelling non-parametrically, little seems to have been done regarding the (presumably easier) problem of parametric modelling (although lots have work has been done for more general multi-state models).

Nyangweso and colleagues consider a parameterisation based on the cumulative incidence functions. This makes sense since the interval-censored likelihood can be directly written in terms of the CIFs. They consider a Gompertz model for each CIF taking the form:


Note that this allows for a proportion who will never experience event k. An obvious complication is that a valid set of CIFs needs to have the property that . It is noted that in practice the unconstrained MLE will be such that the sum of the CIFs will be less than 1 at the greatest time at which someone was right censored as otherwise the likelihood involves taking the log of a non-positive number. Some discussion of constrained optimization, where , so that it is assumed that all patients must eventually experience on of the events (i.e. no cure fraction).

In addition to full likelihood estimation, the parametric analogue of the naive estimator investigated in Jewell et al (2003, Biometrika) is also considered. This just involves modelling each CIF separately as if it were univariate interval censored survival data. Obviously there is an even greater risk of CIFs summing to more than 1 for the naive estimates.

The paper does not get as far as considering models for covariates. The main complication with adding covariates to the Gompertz model seems to be that we would then be almost guaranteed to have covariate patterns where the CIFs sum to more than 1 even at times of interest. There is surely advantage in modelling the cause specific hazards as this way the CIFs are guaranteed to be valid at all times for all covariate values. While for most parametric hazards computation of the CIF requires numerical integration, the extra computation required shouldn't be prohibitive.

Update: The original paper appears to have been superseded by a new version with a slightly different set of authors (Hudgens, Li & Fine).

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.

Estimating the distribution of the window period for recent HIV infections: A comparison of statistical methods

Michael Sweeting, Daniela De Angelis, John Parry and Barbara Suligo have a new paper in Statistics in Medicine. This considers estimating the time to a biomarker crossing a threshold from seroconversion in doubly-censored data (denoted T). Such methods have been used in the past to model time to seroconversion from infection. Sweeting et al state in the introduction that a key motivation is to allow the prevalence of recent infection to be estimated. This is defined as





where S(t) is the time from seroconversion to biomarker crossing. Here we see there is a clear assumption that T is independent of calendar time d. This is consistent with the framework first proposed by De Gruttola and Lagakos (Biometrics, 1989) where the time to first event is assumed to be independent of the time between first and second events, the data are essentially a three-state progressive semi-Markov model where interest lies in estimating the sojourn distribution in state 2. Alternatively, the methods of Frydman (JRSSB, 1992) are based on a Markov assumption. Here the marginal distribution of the sojourn in state 2 could be estimated by integrating over the estimated distribution of times to entry into state 2. Sweeting et al however treat it as standard bivariate survival data (where X denotes the time to serocoversion and Z denotes time to biomarker crossing) using MLEcens for estimation. It's not clear whether this is sensible if the aim is to estimate P(d) as defined above. Moreover, Betensky and Finklestein (Statistics in Medicine, 1999) suggested an alternative algorithm would be required for doubly-censored data because of the inherent ordering (i.e. Z>X). Presumably as long as the intervals for seroconversion and crossing are disjoint the NPMLE for S is guaranteed to have support only at non-negative times.

Sweeting et al make a great deal about the fact that since the NPMLE of (X,Z) is non-unique, e.g. the NPMLE can only ascribe mass to regions of the (X,Z) plane, and that T = Z - X, there will be huge uncertainty over the distribution of S between the extremes of assuming mass only in the lower left corner or only in the upper right corner of the support rectangles. They provide a plot to show this in their data example. The original approach to doubly-censored data of De Gruttola and Lagakos assumes a set of discrete time points of mass for S (the generalization to higher order models was given by Sternberg and Satten). This largely avoids any problems of non-uniqueness (though some problems remain see e.g. Griffin and Lagakos 2010).

For the non-parametric approach, Sweeting et al seem extremely reluctant to make any assumptions whatsoever. In contrast, they go completely to town on assumptions once they get onto their preferred Bayesian method. It should firstly be mentioned that the outcome time is time to crossing of a biomarker and there is considerable auxiliary data available in the form of intermediate measurements. Thus, under any kind of monotonicity assumptions, we can see that extra modelling of the growth process of the biomarker has merit. Sweeting et al use a parametric, mixed effects growth model, to model the growth of the biomarker. The growth is measured in terms of time since seroconversion, which is unknown. They consider two methods: a naive method that assumes seroconversion occurs at the midpoint of the time interval and a uniform prior method that assumes a priori that the seroconversion time is distributed uniformly within the times at which it is interval censored (i.e. last time before seroconversion and the first biomarker measurement time). The first method is essentially like imputing the interval midpoint for X. The second method is like assuming the marginal distribution for X is uniform.

Overall, the paper comes across as a straw man argument against non-parametric methods appended onto a Bayesian analysis that makes multiple (unverified) assumptions.

A semi-competing risks model for data with interval-censoring and informative observation

Jessica Barrett, Fotios Siannis and Vern Farewell have a new paper in Statistics in Medicine. Essentially the paper uses similar methods as in Siannis et al 2006, to investigate informative loss-to-follow-up (LTF) in a study of aging and cognitive function.

LTF, refers to loss-to-follow-up of monitoring cognitive impairment - crucially survival continues to be monitored. LTF (from healthy) is modelled as a separate state in the process with its own transition intensity. Once LTF, the subject experiences different intensities of becoming cognitive impaired or dying. An unidentifiable parameter k determines the relative rate at which people who are lost to follow-up (before becoming cognitively impaired) proceed to the cognitively impaired state rather than the death state, compared to those not lost to follow-up.
k can be varied to see what impact assumptions about those LTF have on overall estimates. In the current study k has quite a large impact on estimates of cumulative incidence of cognitive impairment. This is in contrast to the Whitehall study where these methods were applied to right-censored data, where k had little effect.

A parametric Weibull intensities Markov model is used to model the data. Due to the interval censoring, computation of the likelihood requires numerical integration.

As an informal goodness-of-fit test the authors compare on Cox model of overall survival, with the corresponding survival estimates for the multi-state model with proportional intensity models on each intensity. The authors note that the Cox proportional hazards model for overall survival should be unbiased. Of course, it is only unbiased if the proportional hazards assumption holds on the overall hazard of death. If, however, the covariates are proportional on the individual intensities, as assumed in the multi-state model, the Cox model on overall survival will be biased. This approach is therefore more of a test of robustness to assumptions about the covariates than a goodness-of-fit test because the models aren't nested. The same method was used in Siannis et al 2006 and in Van den Hout et al, 2009.

Multiple imputation for estimating the risk of developing dementia and its impact on survival

Yu, Saczynski and Launer have a new paper in Biometrical Journal. This proposes an approach to fitting Cox-Markov or Cox-semi-Markov models to illness-death models subject to interval censoring via the use of iterative multiple imputation.

In essence this paper is a generalization of the method for interval censored survival data proposed by Pan (2000; Biometrics). The idea is to start by imputing the mid-point of the interval for unknown event times, calculate the parameter estimates treating these times as known. Then based on the new estimates, impute e.g. 20 different complete datasets, estimate the parameters based on these complete datasets, and the re-impute. This process continues until the parameter estimates converge.

The authors consider data regarding mortality and the onset of dementia. They model post-dementia hazard as being proportional to pre-dementia hazard and additionally allow onset age to affect hazard (for the semi-Markov model).

A drawback of the approach, which was mentioned in the paper by Pan, but not seemingly referred to by Yu et al, is that the support points for the baseline intensities will depend entirely on the initial values chosen, i.e. the final baseline intensities will have support at a subset of the original support points. For small datasets, there may be a decay in the number of support points down to 1. Pan suggested using a Link estimate of the baseline hazard (i.e. linear smoothed Breslow estimate) and to impute estimates from this to avoid the problem.

A further issue is determining what exactly constitutes convergence of the parameter estimates. After all, at each iteration a random imputation of M datasets takes place. Hence, we can expect the new parameter estimates to have a certain degree of random variation compared to the previous iteration.

The authors should be commended for having provided code and a subset of the full data in the Supporting Material of the paper. However, there seems to be an error in the code which results in an overestimation of the cumulative disease onset hazard. The code obtains the Breslow estimate of the cumulative baseline hazard. It then estimates the hazard between event times to be But instead of then simulating new event times from a piecewise-exponential distribution, they ascribe h(t) to all event times (e.g. previously imputed onset times and death times) between t_i and t_{i+1}. The coding error seems to also result in considerable bias in the other estimates as corrected code (assuming support for the hazard only at imputed onset times) seems to give quite different results.

The analysis seems to be incorrect in that it treats all patients as left-truncated (i.e. healthy) at age 50 (presumably time 0 in the data is age 50). No deaths occur before age 71.6 (or 21.6 on the time scale used in the data). Even prevalent cases enter the study conditional on survival until their entry age.

A comparison with existing methods can be made by computing the NPMLE under a Markov assumption. For this we use methods developed by Gauzere (2000; PhD Thesis) and Frydman and Szarek (2009; Biometrics). For data removing the prevalent cases, a corrected version of the imputation code gives estimates that are reasonably similar to the NPMLE (see figure), whereas the uncorrected version substantially overestimates the cumulative incidence of dementia.



The results of the analysis of the data in the paper are close to those using the subsample data with the code, suggesting the same code was used for the full data. It is unclear whether the same code was used for the simulation data.

Finally, it should be noted that the method, while undoubtedly more straightforward to implement than either a semi-parametric maximum likelihood estimate or the penalized likelihood approach of Joly and Commenges, is nevertheless still quite slow. Moreover, as with the other methods, the complexity increases rapidly as the number of states in the model increases.

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 novel semiparametric regression method for interval-censored data

Seungbong Han, Adin-Cristian Andrei and Kam-Wah Tsui have a paper currently available as a University of Wisconsin Biostatistics and Medical Informatics Department working paper. This essentially extends the concept of pseudo-observations, which up till now have only concerned right-censored, to interval censored survival data. The idea remains the same except that S(t) is estimated using the NPMLE of the survival function (e.g. via Turnbull self-consistent estimator or iterative convex minorant).

The paper is a little disappointing in only providing a very brief heuristic justification for using pseudo-observations in the interval-censoring case. For right-censored data, an estimate of the baseline survival function can be obtained as well as the regression parameters. No discussion of whether this is possible for the interval censored case is given. However, the baseline estimates are likely to be highly unreliable (e.g. non-monotonic) because particular subjects may have extreme influence because they effect where the mass points of the NPMLE occur. For example, the plot above is based on 1000 subjects with survival generated from an exponential with rate 0.25, subject to independent current status observation (uniformly on (0,10)). Estimating the baseline survival from the pseudo-observations (calcuated at times 1,2,3,...,10) leads to a survivor function which increases at time 5. It seems necessary that there should be more consideration of this issue as well as the choice of how many time points to evaluate the pseudo-observations at.

The authors choose to transform the pseudo-observations before regressing on the covariates, rather than using a link function in a GLM. One problem with this approach is presumably that if the estimate of S(t) is 0 or 1, g(S(t))=-Inf or Inf.

On a practical level the authors use Icens to calculate the NPMLE. As noted previously the MLEcens package seems to perform considerably better than Icens and would presumably speed up computation of the pseudo-observations method.

A natural next step would be to consider pseudo-observations for interval censored multi-state data. The lack of non-parametric methods except in a few simple cases is an obvious bar to development in this direction.

Update: The paper has now been published in Communications in Statistics - Simulation and Computation.