Tutorial in biostatistics: competing risks and multi-state models

Standard survival data measure the time span from some time origin until the occurrence of one type of event. If several types of events occur, a model describing progression to each of these competing risks is needed. Multi-state models generalize competing risks models by also describing transitions to intermediate events. Methods to analyze such models have been developed over the last two decades. Fortunately, most of the analyzes can be performed within the standard statistical packages, but may require some extra effort with respect to data preparation and programming. This tutorial aims to review statistical methods for the analysis of competing risks and multi-state models. Although some conceptual issues are covered, the emphasis is on practical issues like data preparation, estimation of the effect of covariates, and estimation of cumulative incidence functions and state and transition probabilities. Examples of analysis with standard software are shown.     

Modeling disease‐state transition heterogeneity through Bayesian variable selection

In many diseases, Markov transition models are useful in describing transitions between discrete disease states. Often the probability of transitioning from one state to another varies widely across subjects. This heterogeneity is driven, in part, by a possibly unknown number of previous disease states and by potentially complex relationships between clinical data and these states. We propose use of Bayesian variable selection in Markov transition models to allow estimation of subject‐specific transition probabilities. Our approach simultaneously estimates the order of the Markov process and the transition‐specific covariate effects. The methods are assessed using simulation studies and applied to model disease‐state transition on the expanded disability status scale (EDSS) in multiple sclerosis (MS) patients from the Partners MS Center in Boston, MA. The proposed methodology is shown to accurately identify complex covariate–transition relationships in simulations and identifies a clinically significant interaction between relapse history and EDSS history in MS patients. Copyright © 2009 John Wiley & Sons, Ltd.     

Structured fusion lasso penalized multi‐state models

Multi‐state models generalize survival or duration time analysis to the estimation of transition‐specific hazard rate functions for multiple transitions. When each of the transition‐specific risk functions is parametrized with several distinct covariate effect coefficients, this leads to a model of potentially high dimension. To decrease the parameter space dimensionality and to work out a clear image of the underlying multi‐state model structure, one can either aim at setting some coefficients to zero or to make coefficients for the same covariate but two different transitions equal. The first issue can be approached by penalizing the absolute values of the covariate coefficients as in lasso regularization. If, instead, absolute differences between coefficients of the same covariate on different transitions are penalized, this leads to sparse competing risk relations within a multi‐state model, that is, equality of covariate effect coefficients. In this paper, a new estimation approach providing sparse multi‐state modelling by the aforementioned principles is established, based on the estimation of multi‐state models and a simultaneous penalization of the L₁‐norm of covariate coefficients and their differences in a structured way. The new multi‐state modelling approach is illustrated on peritoneal dialysis study data and implemented in the R package penMSM . Copyright © 2016 John Wiley & Sons, Ltd.     

Parametric multistate survival models: Flexible modelling allowing transition‐specific distributions with application to estimating clinically useful measures of effect differences

Multistate models are increasingly being used to model complex disease profiles. By modelling transitions between disease states, accounting for competing events at each transition, we can gain a much richer understanding of patient trajectories and how risk factors impact over the entire disease pathway. In this article, we concentrate on parametric multistate models, both Markov and semi‐Markov, and develop a flexible framework where each transition can be specified by a variety of parametric models including exponential, Weibull, Gompertz, Royston‐Parmar proportional hazards models or log‐logistic, log‐normal, generalised gamma accelerated failure time models, possibly sharing parameters across transitions. We also extend the framework to allow time‐dependent effects. We then use an efficient and generalisable simulation method to calculate transition probabilities from any fitted multistate model, and show how it facilitates the simple calculation of clinically useful measures, such as expected length of stay in each state, and differences and ratios of proportion within each state as a function of time, for specific covariate patterns. We illustrate our methods using a dataset of patients with primary breast cancer. User‐friendly Stata software is provided.     

Flexible multistate models for interval‐censored data: Specification,estimation, and an application to ageing research

Continuous‐time multistate survival models can be used to describe health‐related processes over time. In the presence of interval‐censored times for transitions between the living states, the likelihood is constructed using transition probabilities. Models can be specified using parametric or semiparametric shapes for the hazards. Semiparametric hazards can be fitted using P‐splines and penalised maximum likelihood estimation. This paper presents a method to estimate flexible multistate models that allow for parametric and semiparametric hazard specifications. The estimation is based on a scoring algorithm. The method is illustrated with data from the English Longitudinal Study of Ageing.     

Nonparametric estimation of the cumulative incidence function under outcome misclassification using external validation data

Misclassification of outcomes or event types is common in health sciences research and can lead to serious bias when estimating the cumulative incidence functions in settings with competing risks. Recent work has shown how to estimate nonparametric cumulative incidence functions in the presence of nondifferential outcome misclassification when the misclassification probabilities are known. Here, we extend this approach to account for misclassification that is differential with respect to important predictors of the outcome using misclassification probabilities estimated from external validation data. Moreover, we propose a bootstrap approach in which the observations from both the main study data and the external validation study are resampled to allow the uncertainty in the misclassification probabilities to propagate through the analysis into the final confidence intervals, ensuring appropriate confidence interval coverage probabilities. The proposed estimator is shown to be uniformly consistent and simulation studies indicate that both the estimator and the standard error estimation approach perform well in finite samples. The methodology is applied to estimate the cumulative incidence of death and disengagement from HIV care in a large cohort of HIV infected individuals in sub-Saharan Africa, where a significant death underreporting issue leads to outcome misclassification. This analysis uses external validation data from a separate study conducted in the same country.     

A marginal structural model for multiple‐outcome survival data:assessing the impact of injection drug use on several causes of death in the Canadian Co‐infection Cohort

It is often the case that interest lies in the effect of an exposure on each of several distinct event types. For example, we are motivated to investigate in the impact of recent injection drug use on deaths due to each of cancer, end‐stage liver disease, and overdose in the Canadian Co‐infection Cohort (CCC). We develop a marginal structural model that permits estimation of cause‐specific hazards in situations where more than one cause of death is of interest. Marginal structural models allow for the causal effect of treatment on outcome to be estimated using inverse‐probability weighting under the assumption of no unmeasured confounding; these models are particularly useful in the presence of time‐varying confounding variables, which may also mediate the effect of exposures. An asymptotic variance estimator is derived, and a cumulative incidence function estimator is given. We compare the performance of the proposed marginal structural model for multiple‐outcome data to that of conventional competing risks models in simulated data and demonstrate the use of the proposed approach in the CCC. Copyright © 2013 John Wiley & Sons, Ltd.     

Comparison of variance estimation approaches in a two-state Markov model for longitudinal data with misclassification

We examine the behaviour of the variance-covariance parameter estimates in an alternating binary Markov model with misclassification. Transition probabilities specify the state transitions for a process that is not directly observable. The state of an observable process, which may not correctly classify the state of the unobservable process, is obtained at discrete time points. Misclassification probabilities capture the two types of classification errors. Variance components of the estimated transition parameters are calculated with three estimation procedures: observed information, jackknife, and bootstrap techniques. Simulation studies are used to compare variance estimates and reveal the effect of misclassification on transition parameter estimation. The three approaches generally provide similar variance estimates for large samples and moderate misclassification. In these situations, the resampling methods are reasonable alternatives when programming partial derivatives is not appealing. With smaller chains or higher misclassification probabilities, the bootstrap method appears to be the best choice.     

Self-reported health status and latent health dynamics

Many dynamic structural models that include health as a state variable use categorical self-reported health status (SRHS) as their sole empirical measure of health. Commonly, transition probabilities among discrete health states are calculated directly from one-wave-ahead transitions observed in panel data. The predicted dynamics of SRHS generated from these calculations rapidly deviate from empirical transitions more than one wave ahead, as the assumption of a Markovian process is not satisfied. Consequently, models that treat SRHS as a true representation of health are not well calibrated to accurately match the long term dynamics of individual health. This paper specifies a model in which SRHS is a noisy measure of a continuous latent health state. I estimate the model by maximum likelihood on several panel datasets, exploiting long sequences of SRHS rather than simple wave-to-wave transitions. By accounting for transitory reporting error, the latent health model is able to match both short-run and long-run SRHS transitions up to twenty years in the future, as well as other features of empirical SRHS dynamics. Moreover, the measure of latent health is at least as good at predicting other outcomes (including mortality, labor supply, and medical expenses) as SRHS itself. To aid structural modelers in replacing their representation of health, I provide simple software tools for generating a discretized latent health process, and for filtering observed sequences of SRHS into (a distribution of) latent health.     

Alternative approaches for econometric analysis of panel count data using dynamic latent class models (with application to doctor visits data)

Cross-sectional latent class regression models, also known as switching regressions or hidden Markov models, cannot identify transitions between classes that may occur over time. This limitation can potentially be overcome when panel data are available. For such data, we develop a sequence of models that combine features of the static cross-sectional latent class (finite mixture) models with those of hidden Markov models. We model the probability of movement between categories in terms of a Markovian structure, which links the current state with a previous state, where state may refer to the category of an individual. This article presents a suite of mixture models of varying degree of complexity and flexibility for use in a panel count data setting, beginning with a baseline model which is a two-component mixture of Poisson distribution in which latent classes are fixed and permanent. Sequentially, we extend this framework (i) to allow the mixing proportions to be smoothly varying continuous functions of time-varying covariates, (ii) to add time dependence to the benchmark model by modeling the class-indicator variable as a first-order Markov chain and (iii) to extend item (i) by making it dynamic and introducing covariate dependence in the transition probabilities. We develop and implement estimation algorithms for these models and provide an empirical illustration using 1995-1999 panel data on the number of doctor visits derived from the German Socio-Economic Panel.     

Joint modelling of longitudinal and multi‐state processes: application to clinical progressions in prostate cancer

Joint modelling of longitudinal and survival data is increasingly used in clinical trials on cancer. In prostate cancer for example, these models permit to account for the link between longitudinal measures of prostate‐specific antigen (PSA) and time of clinical recurrence when studying the risk of relapse. In practice, multiple types of relapse may occur successively. Distinguishing these transitions between health states would allow to evaluate, for example, how PSA trajectory and classical covariates impact the risk of dying after a distant recurrence post‐radiotherapy, or to predict the risk of one specific type of clinical recurrence post‐radiotherapy, from the PSA history. In this context, we present a joint model for a longitudinal process and a multi‐state process, which is divided into two sub‐models: a linear mixed sub‐model for longitudinal data and a multi‐state sub‐model with proportional hazards for transition times, both linked by a function of shared random effects. Parameters of this joint multi‐state model are estimated within the maximum likelihood framework using an EM algorithm coupled with a quasi‐Newton algorithm in case of slow convergence. It is implemented under R, by combining and extending mstate and JM packages. The estimation program is validated by simulations and applied on pooled data from two cohorts of men with localized prostate cancer. Thanks to the classical covariates available at baseline and the repeated PSA measurements, we are able to assess the biomarker's trajectory, define the risks of transitions between health states and quantify the impact of the PSA dynamics on each transition intensity. Copyright © 2016 John Wiley & Sons, Ltd.     

A global logrank test for adaptive treatment strategies based on observational studies

In studying adaptive treatment strategies, a natural question that is of paramount interest is whether there is any significant difference among all possible treatment strategies. When the outcome variable of interest is time‐to‐event, we propose an inverse probability weighted logrank test for testing the equivalence of a fixed set of pre‐specified adaptive treatment strategies based on data from an observational study. The weights take into account both the possible selection bias in an observational study and the fact that the same subject may be consistent with more than one treatment strategy. The asymptotic distribution of the weighted logrank statistic under the null hypothesis is obtained. We show that, in an observational study where the treatment selection probabilities need to be estimated, the estimation of these probabilities does not have an effect on the asymptotic distribution of the weighted logrank statistic, as long as the estimation of the parameters in the models for these probabilities is ‐consistent. Finite sample performance of the test is assessed via a simulation study. We also show in the simulation that the test can be pretty robust to misspecification of the models for the probabilities of treatment selection. The method is applied to analyze data on antidepressant adherence time from an observational database maintained at the Department of Veterans Affairs’ Serious Mental Illness Treatment Research and Evaluation Center. Copyright © 2013 John Wiley & Sons, Ltd.     

Estimation of markov chain transition probabilities and rates from fully and partially observed data: uncertainty propagation, evidence synthesis, and model calibration.

Markov transition models are frequently used to model disease progression. The authors show how the solution to Kolmogorov's forward equations can be exploited to map between transition rates and probabilities from probability data in multistate models. They provide a uniform, Bayesian treatment of estimation and propagation of uncertainty of transition rates and probabilities when 1) observations are available on all transitions and exact time at risk in each state (fully observed data) and 2) observations are on initial state and final state after a fixed interval of time but not on the sequence of transitions (partially observed data). The authors show how underlying transition rates can be recovered from partially observed data using Markov chain Monte Carlo methods in WinBUGS, and they suggest diagnostics to investigate inconsistencies between evidence from different starting states. An illustrative example for a 3-state model is given, which shows how the methods extend to more complex Markov models using the software WBDiff to compute solutions. Finally, the authors illustrate how to statistically combine data from multiple sources, including partially observed data at several follow-up times and also how to calibrate a Markov model to be consistent with data from one specific study.     

Network Meta-Analysis with Competing Risk Outcomes

BackgroundCost-effectiveness analysis often requires information on the effectiveness of interventions on multiple outcomes, and commonly these take the form of competing risks. Nevertheless, methods for synthesis of randomized controlled trials with competing risk outcomes are limited.ObjectiveThe aim of this study was to develop and illustrate flexible evidence synthesis methods for trials reporting competing risk results, which allow for studies with different follow-up times, and that take account of the statistical dependencies between outcomes, regardless of the number of outcomes and treatments.MethodsWe propose a competing risk meta-analysis based on hazards, rather than probabilities, estimated in a Bayesian Markov chain Monte Carlo (MCMC) framework using WinBUGS software. Our approach builds on existing work on mixed treatment comparison (network) meta-analysis, which can be applied to any number of treatments, and any number of competing outcomes, and to data sets with varying follow-up times. We show how a fixed effect model can be estimated, and two random treatment effect models with alternative structures for between-trial variation. We suggest methods for choosing between these alternative models.ResultsWe illustrate the methods by applying them to a data set involving 17 trials comparing nine antipsychotic treatments for schizophrenia including placebo, on three competing outcomes: relapse, discontinuation because of intolerable side effects, and discontinuation for other reasons.ConclusionsBayesian MCMC provides a flexible framework for synthesis of competing risk outcomes with multiple treatments, particularly suitable for embedding within probabilistic cost-effectiveness analysis.     

Prediction of hemoglobin in blood donors using a latent class mixed‐effects transition model

Blood donors experience a temporary reduction in their hemoglobin (Hb) value after donation. At each visit, the Hb value is measured, and a too low Hb value leads to a deferral for donation. Because of the recovery process after each donation as well as state dependence and unobserved heterogeneity, longitudinal data of Hb values of blood donors provide unique statistical challenges. To estimate the shape and duration of the recovery process and to predict future Hb values, we employed three models for the Hb value: (i) a mixed‐effects models; (ii) a latent‐class mixed‐effects model; and (iii) a latent‐class mixed‐effects transition model. In each model, a flexible function was used to model the recovery process after donation. The latent classes identify groups of donors with fast or slow recovery times and donors whose recovery time increases with the number of donations. The transition effect accounts for possible state dependence in the observed data. All models were estimated in a Bayesian way, using data of new entrant donors from the Donor InSight study. Informative priors were used for parameters of the recovery process that were not identified using the observed data, based on results from the clinical literature. The results show that the latent‐class mixed‐effects transition model fits the data best, which illustrates the importance of modeling state dependence, unobserved heterogeneity, and the recovery process after donation. The estimated recovery time is much longer than the current minimum interval between donations, suggesting that an increase of this interval may be warranted. Copyright © 2015 John Wiley & Sons, Ltd.     

Flexible parametric modelling of the cause‐specific cumulative incidence function

Competing risks arise with time‐to‐event data when individuals are at risk of more than one type of event and the occurrence of one event precludes the occurrence of all other events. A useful measure with competing risks is the cause‐specific cumulative incidence function (CIF), which gives the probability of experiencing a particular event as a function of follow‐up time, accounting for the fact that some individuals may have a competing event. When modelling the cause‐specific CIF, the most common model is a semi‐parametric proportional subhazards model. In this paper, we propose the use of flexible parametric survival models to directly model the cause‐specific CIF where the effect of follow‐up time is modelled using restricted cubic splines. The models provide smooth estimates of the cause‐specific CIF with the important advantage that the approach is easily extended to model time‐dependent effects. The models can be fitted using standard survival analysis tools by a combination of data expansion and introducing time‐dependent weights. Various link functions are available that allow modelling on different scales and have proportional subhazards, proportional odds and relative absolute risks as particular cases. We conduct a simulation study to evaluate how well the spline functions approximate subhazard functions with complex shapes. The methods are illustrated using data from the European Blood and Marrow Transplantation Registry showing excellent agreement between parametric estimates of the cause‐specific CIF and those obtained from a semi‐parametric model. We also fit models relaxing the proportional subhazards assumption using alternative link functions and/or including time‐dependent effects. Copyright © 2016 John Wiley & Sons, Ltd.     

Reduced-rank proportional hazards regression and simulation-based prediction for multi-state models

In this paper we address two issues arising in multi-state models with covariates. The first issue deals with how to obtain parsimony in the modeling of the effect of covariates. The standard way of incorporating covariates in multi-state models is by considering the transitions as separate building blocks, and modeling the effect of covariates for each transition separately, usually through a proportional hazards model for the transition hazard. This typically leads to a large number of regression coefficients to be estimated, and there is a real danger of over-fitting, especially when transitions with few events are present. We extend the reduced-rank ideas, proposed earlier in the context of competing risks, to multi-state models, in order to deal with this issue. The second issue addressed in this paper was motivated by the wish to obtain standard errors of the regression coefficients of the reduced-rank model. We propose a model-based resampling technique, based on repeatedly sampling trajectories through the multi-state model. The same ideas are also used for the estimation of prediction probabilities in general multi-state models and associated standard errors.We use data from the European Group for Blood and Marrow Transplantation to illustrate our techniques.     

Semiparametric methods for multistate survival models in randomised trials

Transform methods have proved effective for networks describing a progression of events. In semi‐Markov networks, we calculated the transform of time to a terminating event from corresponding transforms of intermediate steps. Saddlepoint inversion then provided survival and hazard functions, which integrated, and fully utilised, the network data. However, the presence of censored data introduces significant difficulties for these methods. Many participants in controlled trials commonly remain event‐free at study completion, a consequence of the limited period of follow‐up specified in the trial design. Transforms are not estimable using nonparametric methods in states with survival truncated by end‐of‐study censoring. We propose the use of parametric models specifying residual survival to next event. As a simple approach to extrapolation with competing alternative states, we imposed a proportional incidence (constant relative hazard) assumption beyond the range of study data. No proportional hazards assumptions are necessary for inferences concerning time to endpoint; indeed, estimation of survival and hazard functions can proceed in a single study arm. We demonstrate feasibility and efficiency of transform inversion in a large randomised controlled trial of cholesterol‐lowering therapy, the Long‐Term Intervention with Pravastatin in Ischaemic Disease study. Transform inversion integrates information available in components of multistate models: estimates of transition probabilities and empirical survival distributions. As a by‐product, it provides some ability to forecast survival and hazard functions forward, beyond the time horizon of available follow‐up. Functionals of survival and hazard functions provide inference, which proves sharper than that of log‐rank and related methods for survival comparisons ignoring intermediate events. Copyright © 2013 John Wiley & Sons, Ltd.     

Estimation of life‐years gained and cost effectiveness based on cause‐specific mortality

Cost‐effectiveness analysis is usually based on life‐years gained estimated from all‐cause mortality. When an intervention affects only a few causes of death accounting for a small fraction of all deaths, this approach may lack precision. We develop a novel technique for cost‐effectiveness analysis when life‐years gained are estimated from cause‐specific mortality, allowing for competing causes of death. In the context of randomised trial data, we adjust for other‐cause mortality combined across randomised groups. This method yields a greater precision than analysis based on total mortality, and we show application to life‐years gained, quality‐adjusted life‐years gained, incremental costs, and cost effectiveness. In multi‐state health economic models, however, mortality from competing causes is commonly derived from national statistics and is assumed to be known and equal across intervention groups. In such models, our method based on cause‐specific mortality and standard methods using total mortality give essentially identical estimates and precision. The methods are applied to a randomised trial and a health economic model, both of screening for abdominal aortic aneurysm. A gain in precision for cost‐effectiveness estimates is clearly helpful for decision making, but it is important to ensure that ‘cause‐specific mortality’ is defined to include all causes of death potentially affected by the intervention. Copyright © 2010 John Wiley & Sons, Ltd.     

Fitting and interpreting continuous‐time latent Markov models for panel data

Multistate models characterize disease processes within an individual. Clinical studies often observe the disease status of individuals at discrete time points, making exact times of transitions between disease states unknown. Such panel data pose considerable modeling challenges. Assuming the disease process progresses accordingly, a standard continuous‐time Markov chain (CTMC) yields tractable likelihoods, but the assumption of exponential sojourn time distributions is typically unrealistic. More flexible semi‐Markov models permit generic sojourn distributions yet yield intractable likelihoods for panel data in the presence of reversible transitions. One attractive alternative is to assume that the disease process is characterized by an underlying latent CTMC, with multiple latent states mapping to each disease state. These models retain analytic tractability due to the CTMC framework but allow for flexible, duration‐dependent disease state sojourn distributions. We have developed a robust and efficient expectation–maximization algorithm in this context. Our complete data state space consists of the observed data and the underlying latent trajectory, yielding computationally efficient expectation and maximization steps. Our algorithm outperforms alternative methods measured in terms of time to convergence and robustness. We also examine the frequentist performance of latent CTMC point and interval estimates of disease process functionals based on simulated data. The performance of estimates depends on time, functional, and data‐generating scenario. Finally, we illustrate the interpretive power of latent CTMC models for describing disease processes on a dataset of lung transplant patients. We hope our work will encourage wider use of these models in the biomedical setting. Copyright © 2013 John Wiley & Sons, Ltd.