Ties between event times and jump times in the Cox model

Xin, Horrocks and Darlington have a new paper in Statistics in Medicine. This considers approaches for dealing with ties in Cox proportional hazard models, not between event times but between event times and changes to time dependent covariates.

If a change in a time-dependent covariate coincides with a failure time there is ambiguity over which value of the time dependent covariate, z(t+) or z(t-), should be taken for the risk set at time t. By convention, it is usually assumed that z(t-) should be taken, i.e. that the change in the covariate occurs after the failure time. The authors demonstrate that for small sample sizes and/or a large proportion of ties, the estimates can be sensitive to the convention chosen. The authors also only consider cases where z(t) is a binary indicator that jumps from 0 to 1 at some point and cannot make the reverse jump. Obviously this will magnify the potential for bias because the "change after" convention will always underestimate the true risk whereas the "change before" will always overestimate the true risk.

The authors consider some simple adjustments for the problem: compute the "change before" and "change after" estimates and take their average or use random jittering. A problem with the averaging approach is estimating the standard error of the resulting estimator. An upper bound can be obtained by assuming the two estimators have perfect correlation. The jittering estimator obviously has the problem that different random jitters will give different results, though in principle the jittering could be repeated multiple times and combined in a fashion akin to multiple imputation.

It is surprising that the further option of adopting an method akin to the Efron method for ties. Essentially at each failure time there is an associated risk set. It could be argued that every tied covariate jump time had a 50% chance of occurring before or after the failure time. The expected contribution from a particular risk set could then be
It should also be possible to apply this approach using standard software, e.g. coxph() in R. It is simply necessary to replace any (start,stop) interval that ends with a tied "stop" with two intervals (start, stop - 0.00001) and (start, stop + 0.00001) each of which are associated with a weight of 0.5.

Assessing age-at-onset risk factors with incomplete covariate current status data under proportional odds models

Chi-Chung Wen and Yi-Hau Chen have a new paper in Statistics in Medicine. This considers estimation of a proportional odds regression model for current status data in cases where a subset of the covariates may be missing at random for a subset of the patient population.

It is assumed that the probability that a portion of the covariates is missing depends on all the other observable outcomes (the failure status, the survey time and the rest of the covariate vector). The authors propose to fit a logistic regression model, involving all subjects in the dataset, for this probability of missingness. To fit the regression model for the current status data itself, they propose to use what they term a "validation likelihood estimator." This involves only working with the subset of patients with complete data but maximizing a likelihood that conditions on the fact that the whole covariate vector was observed. An advantage of using the proportional odds model over other candidate models (e.g. proportional hazards) is that the resulting likelihood remains of the proportional odds form.

Clearly a disadvantage of this "validation likelihood estimator" is that the data from subjects who have incomplete covariates is not used directly in the regression model. As a result the estimator is likely to be less efficient than approaches that effectively attempt to impute the missing covariate values. The authors argue that the validation likelihood approach will tend to be more robust since it is not necessary to make (parametric) assumptions about the conditional distribution of the missing covariates.

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

Practicable confidence intervals for current status data

Byeong Yeob Choi, Jason P. Fine and M. Alan Brookhart have a new paper in Statistics in Medicine. Essentially the paper clarifies the practical implications of results relating to the asymptotic theory for current status data. In particular, it is known that the nonparametric bootstrap is inconsistent for current status data when the distribution of sampling times is continuous. The authors note that the most reliable method considered previously is a previous study by Ghosh et al concluded that construction of confidence intervals based on inversion of the likelihood ratio statistic (as original proposed in Banerjee and Wellner (2001)) gave the best results particularly for smaller sample sizes. However, they also note that this approach is difficult to implement (e.g. lack of available software). Here they therefore pursue approaches based on using the limiting Chernoff distribution to construct Wald type confidence intervals, but using cloglog or logit transformations to get better coverage, and also look at the performance of simple non-parametric bootstrap.

Perhaps unsurprisingly they find that, when sample sizes are relatively small, the performance of all methods is dependent on the quantile of the failure time distribution at which the confidence interval is computed (e.g. it performs much better when t is close to the median) and the observation density at the time point considered (performance is poorer at times with a lower observation density). Using cloglog or logit transformations is found to improve coverage, but the non-parametric bootstrap tended to outperform this approach for smaller sample sizes, suggesting the admissibility of using the non-parametric bootstrap.

An apparent omission in the paper is any mention of the altered asymptotics in the case where the observation distribution has support at a finite set of time points (or indeed where the rate of increase of points is less than ). This issue is most comprehensively discussed in Tang et al which didn't come out until after the paper was apparently submitted. However, the basic issue of standard asymptotics (and by implication a consistent bootstrap) when there is a finite set of observation points is discussed in Maathuis and Hugdens (2011) which the authors cite. For instance, in the Hoel and Walburg mice dataset used as illustration, the resolution of the data is to the nearest day. It is therefore reasonable to assume that in this case were the sample size to increase, the number of observations would either be bounded by a fixed value (e.g. ~1000) or else the number of unique points would increase at a rate much less than .

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.

Absolute risk regression for competing risks: interpretation, link functions, and prediction

Thomas Gerds, Thomas Scheike and Per Andersen have a new paper in Statistics in Medicine. To a certain extent this is a review paper and considers models for direct regression on the cumulative incidence function for competing risks data. Specifically models of the form where is a known link function and is the cumulative incidence function for event 1 given covariates X. The Fine-Gray model is a special case of this class of models, where a complementary log-log link is adopted. Approaches to estimation based on inverse probability of censoring weights and jackknife based pseudo-observations are considered. Model comparison based on predictive accuracy as measured through Brier score and model diagnostics based on extended models allowing time dependent covariate effects are also discussed.
The discussion gives a clear account of the various pros and cons of direct regression of the cumulative incidence functions. In particular, an obvious, although perhaps not always sufficiently emphasized issue is that if, in a model with two causes, a Fine-Gray (or other direct model) is fitted to the first cause, and another to the second cause, the resulting predictions will not necessarily have the property that an issue that is not problematic if the second cause is essentially a nuisance issue, but obviously problematic if both causes are of interest. In such cases regression of the cause-specific-hazards is preferable even if it makes interpreting the effect on the cumulative intensity functions more difficult.

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.

Nonparametric multistate representations of survival and longitudinal data with measurement error

Bo Hu, Liang Li, Xiaofeng Wang and Tom Greene have a new paper in Statistics in Medicine. This develops approaches for summarizing longitudinal and survival data in terms of marginal prevalences. The authors' use of "multistate" is perhaps not entirely in line with its typical usage. They consider data consisting of right censored competing risks data plus additional continuous longitudinal measurements which persist until a competing event has occurred (or censoring). For the purpose of creating a summary measure, the longitudinal measurement can be partitioned into a set of discrete states. There are thus states corresponding to absorbing competing risks plus a series of transient states corresponding to the longitudinal measurements. The aim of the paper is to develop a nonparametric estimate of the marginal probability of being in a particular state at a particular time.

The approach taken is firstly to use standard non-parametric estimates for competing risks data to get estimates of the probability of being in each of the absorbing states. For the longitudinal part, it is assumed that the "true" longitudinal process is not directly observed but instead observed with measurement error. As a consequence the authors propose to use smoothing splines to get an individual estimate of each subject's true trajectory. The combined state occupancy probability at time t for a longitudinal state then consists of the overall survival probability from the competing risks multiplied by the proportion of subjects still at risk at time t who are estimated (on the basis of their spline smooth) to be within that interval. The probability of being in an absorbing state is computed directly from the competing risks estimates. Overall a stacked probability plot consisting of the stacked CIFs for each of the competing risks plus the (not necessarily monotonic) partition of the longitudinal states.

The use of individual smoothing splines seems to present practical problems. Firstly, it assumes that the true longitudinal process is itself in some way "smooth". In some cases the change in state in a biological process may manifest itself in a rapid collapse of a biomarker. Secondly, it seems to require a relatively large number of longitudinal measurements per person in order to get a reasonable estimate of their "true" process. Presumably the level of the longitudinal measure is likely to have a bearing on the cause-specific-hazards of the competing risks. The occurrence of one of the competing risks is thus informative for the longitudinal process. The authors claim to have got around this by averaging only over people currently in the risk set at time t. However, the longitudinal measurements are intermittent. If they are sparse then someone may be observed at say 1 year, 2 years and then die at 10 years. The method would estimate a smooth spline based on years 1 and 2 and extrapolate up to 10 years not using the fact the subject died at 10 years. Similarly, there might be one or fewer longitudinal observations before a competing event for some patients making estimation of the true trajectory near impossible. Also, the estimator as it stands attempts no weighting to take account of the relative uncertainties about different individuals true trajectories at particular times. Overall as a descriptive tool it may be useful in some circumstances; primarily if subjects has regular longitudinal measurements. In this respect it is similar to the "prevalence counts" (Gentleman et al, 1994) method of obtaining non-parametric prevalence estimates for interval censored multi-state data.

In the appendix, a brief description is given of an approach to allowing the transition probabilities between states to be calculated. They only illustrate the method for a case of going from a longitudinal state to an absorbing state (presumably the procedure for transitions between longitudinal states would be different). Nevertheless, there doesn't seem to be any guarantee that estimated transition probabilities will lie in [0,1].

Use of alternative time scales in Cox proportional hazard models: implications for time-varying environmental exposures

Beth Griffin, Garnet Anderson, Regina Shih and Eric Whitsel have a new paper in Statistics in Medicine. The paper investigates the use of different time scales (e.g. other than either study time or patient age) in cohort studies analysed via Cox proportional hazard models. Of particular focus is the use of calendar time as an alternative time scale, with a motivation in the variation of environmental exposures over time. They perform a simulation study considering two scenarios for the relationship between a time varying environmental exposure variable and calendar year. In the first scenario these are made independent, while the second scenario assumes a linear relationship. As one might expect, when there is no correlation between calendar time and the time dependent environmental exposure, estimates are unbiased regardless of the choice of time scale. When a linear relationship exists then models that account for calendar time, either as the primary time scale or as additional covariates in the model. Again, this isn't necessarily surprising because the model is effectively attempting to include a year effect twice, once in the baseline hazard and again as a large component of the time dependent covariate, e.g. you are fitting a model with "mean environmental exposure in year t" and "environmental exposure" as covariates and expecting the latter to have the correct coefficient. The paper only gives a simulation study, I don't think it would have been that hard to have given some basic theoretical results in addition to the simulations.
The conclusion of the paper, albeit with caveats, is that attempting to adjust for calendar time because you suspect the environmental exposure may be correlated with time is not useful. Clearly if there are other reasons to suspect that calendar year may be important to the hazard in a study then there is an inherent lack of information in the study to establish whether the environmental exposure is directly affecting the hazard or whether it is an indirect effect due to the association with calendar time. Ideally, one would look for other calendar time dependent covariates (e.g. prevailing treatment policy regimes etc.) and perhaps try directly adjusting for them rather than calendar time itself.

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.

Modeling hospital length of stay by Coxian phase-type regression with heterogeneity

Xiaoqin Tang, Zhehui Luo and Joseph Gardiner have a new paper in Statistics in Medicine. This considers modelling right-censored length-of-stay data by using Coxian phase-type distributions, which are distributions defined by the times to absorption of a class of acyclic finite-state time homogeneous Markov processes. Coxian phase-type distributions have been used quite extensively to model right-censored survival data , particularly length-of-stay, see e.g. Marshall and Zenga. The novelty of the current paper is in the estimation of the model. Phase-type distributions suffer from being over parameterized and as a result suffer identifiability issues that in turn cause poor behaviour of optimization procedures. A further issue is the choice of the number of phases of the phase-type distribution. In the current paper a Bayesian reversible jump MCMC approach is taken.

Any acyclic Markov chain can be represented by a Coxian distribution in which the absolute values of the diagonal of the subgenerator matrix are decreasing. If parameters are unrestricted then there are inherent identifiability problems which would hamper the MCMC procedure. To avoid this the authors parameterize based on the first diagonal element and then the ratio (between 0 and 1) of the second element to the first, and so on, with the hazard of absorption from each phase being determined by a proportion (again between 0 and 1) of the diagonal element.

A mixed non-homogeneous hidden Markov model for categorical data, with application to alcohol consumption

Antonello Maruotti and Roberto Rocci have a new paper in Statistics in Medicine. This develops a hidden Markov model for modelling longitudinal data on alcohol consumption in discrete time. The observed data are taken to consist of a three-level ordinal variable denoting whether no drinking (0 drinks), light drinking (1–c drinks), and intense drinking (c+ drinks) occurred in that period of time. The model considered is both time non-homogeneous and mixed, in the sense that there is additional patient level heterogeneity after accounting for covariates. Rather than specifying a continuous distribution for the random effects, the authors adopt the non-parametric mixing distribution approach. Computationally, a finite mixture random effect is much simpler than having a continuous random effect if the random effect is multi-dimensional. However, computation of the full non-parametric maximum likelihood estimate of the mixing distribution is not in itself straightforward. The authors adopt the approach of Aitkin (Statistics and Computing, 1996) which is essentially to work up from a small number of components, performing an EM-algorithm at the fixed level of mixture components. EM based approaches to obtaining the NPMLE of a mixing distribution are known to perform badly and approaches using directional derivatives are preferred (see for instance Wang 2007, JRSS B). The best model, with m components, is assumed to have been reached once taking m+1 components does not produce a better model in terms of AIC or BIC. The main issue with this approach is that the EM algorithm is typically very sensitive to the initial parameter values chosen and prone to fail to find a global maximum. A further danger with these models is to ascribe too great a physical significance to the mixture components estimated.

To reach the final model, choices have to be made regarding: the categorization for the responses (observed level of drinking), the latent Markov states for the HMM (e.g. 2, 3 or 4 latent states), the number of mixture components for the random effect (how many "archetypes" of longitudinal behavior) and the degree of time non-homogeneity of transition probabilities. As a result, while the model is likely to explain the observed data reasonably well, a leap of faith is required to believe the model is an accurate representation of the process of binge drinking/alcoholism.

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.

Estimating net transition probabilities from cross-sectional data with application to risk factors in chronic disease modeling

van de Kassteele, Hoogenveen, Engelfriet, van Baal and Boshuizen have a new paper in Statistics in Medicine. This considers the estimation of the transition probabilities in a non-homogeneous discrete-time Markov model, when the only available information is cross-sectional data, i.e. for each time (or age) we have only a sample of individuals and their state occupancy from which the prevalence at that time can be estimated. Note this type of observation is more extreme than aggregate data, considered for instance by Crowder and Stephens, where we only have prevalences at a series of times but the state occupation counts correspond to the same set of subjects.

The authors take a novel, if slightly quirky approach, to estimation. They firstly use P-splines to smooth the observed prevalences. Having obtained these they then need to translate them into transition probabilities. This is not straightforward since there are more parameters to estimate than degrees of freedom. To get around this problem the authors restrict their estimate to be the values that minimize a transportation problem. Essentially this assigns a "cost" to transitions, penalizing those to further apart states and giving zero cost to remaining in the same state. So gives a solution that aims to maximize the diagonals of the transition probability matrices whilst constraining the prevalences to take their P-spline smoothed values.

What is absent from the paper is formal justification for the approach. Presumably a similar outcome could be achieved by applying a penalized likelihood approach, possibly formulating the problem in continuous time and setting the penalty to be the magnitude of the transition intensities (and possibly their derivatives). However, this would require some calibration to choose the penalty weights and it is not clear how this would be done (the usual approach of cross-validation would not work here).

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.

A proportional hazards regression model for the subdistribution with right-censored and left-truncated competing risks data

Xu Zhang, Mei-Jie Zhang and Jason Fine have a new paper in Statistics in Medicine. This covers the same ground as the paper by Geskus in Biometrics, in developing an approach to fitting the Fine-Gray proportional subdistribution hazard model for competing risks data with left truncated and right censored observations by using inverse probability weights (IPW). Bizarrely, the paper makes no reference at all to the Geskus paper. Presumably this is because the paper was first submitted in 2009 before Geskus's work was published (April 2010). However, it is strange that neither the authors nor the referees became aware of the work in the interim (i.e. acceptance of the paper wasn't until March 2011).

What is interesting is the differences between the approach taken in this paper compared to Geskus. The authors work on the basis that since X = min(T,C) is only observable if X > L, where T is the time of failure, L the time of left truncation and C the time of right censoring, the IPW should be calculated conditional on L < X. Zhang et al use a stabilised weight rather than the IPW to reduce the variability in the original weight. The weights they derive seem quite different to Geskus's as they depend on an estimate of overall survival, which will have to depend on the covariates if the semi-parametric model for the subdistribution hazard is to apply.
The authors suggest using Aalen additive hazard models for the overall survival (thus allowing for time varying covariate effects that can ensure the weights are consistent with the proportional subdistribution hazard model).

Zhang et al start from the general case where the truncation and censoring distributions depend on covariates (but are independent conditional on these covariates), though they only detail non-parametric estimates of the weights. Geskus argued that even if the censoring/truncation distribution depended on covariates that didn't imply it was necessary to include these covariates in the weightings.

Given these discrepancies it would be of interest to contrast and compare the two approaches to the same problem. If both approaches are effective, Geskus's seems preferable because the weights are much easier to calculate.

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.

Estimating and testing for center effects in competing risks

Sandrine Katsahian and Christian Boudreau have a new paper in Statistics in Medicine. This develops methods for including frailty terms within a Fine-Gray competing risks model in order to account for clustering, e.g. effects of different centres.

Since the Fine-Gray model is essentially just a standard Cox proportional hazards regression model with additional time dependent weights, based on the censoring distribution, for individuals who have had a competing event, methods appropriate for standard Cox frailty models can be readily adapted.

Katsahian and Boudreau closely follow the approach taken by Ripatti and Palmgren (Biometrics, 2000). They assume a Gaussian frailty. Computation of the likelihood requires integrating out the frailty terms. Here this is performed using a Laplace approximation. A difficulty with the Laplace approximation is that it still requires the modal value of the frailty distribution conditional on the data and current values of the parameters. The authors therefore take a profile likelihood approach in which they fix the frailty variance and maximize the likelihood term with respect to both the regression parameters and the frailty terms, . Having obtained, and they can then plug into the Laplace approximation to get the profile likelihood for . The procedure gives a local approximation for which can be used to suggest the updated estimate. Thus the process involves alternating between two Newton-Raphson algorithms until convergence.

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.