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.     

Indirect estimation of a discrete‐state discrete‐time model using secondary data analysis of regression data

Multi‐state models of chronic disease are becoming increasingly important in medical research to describe the progression of complicated diseases. However, studies seldom observe health outcomes over long time periods. Therefore, current clinical research focuses on the secondary data analysis of the published literature to estimate a single transition probability within the entire model. Unfortunately, there are many difficulties when using secondary data, especially since the states and transitions of published studies may not be consistent with the proposed multi‐state model. Early approaches to reconciling published studies with the theoretical framework of a multi‐state model have been limited to data available as cumulative counts of progression. This paper presents an approach that allows the use of published regression data in a multi‐state model when the published study may have ignored intermediary states in the multi‐state model. Colloquially, we call this approach the Lemonade Method since when study data give you lemons, make lemonade. The approach uses maximum likelihood estimation. An example is provided for the progression of heart disease in people with diabetes. Copyright © 2009 John Wiley & Sons, Ltd.     

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.     

Incorporation of nested frailties into semiparametric multi‐state models

Proportional hazards models are among the most popular regression models in survival analysis. Multi‐state models generalize them by jointly considering different types of events and their interrelations, whereas frailty models incorporate random effects to account for unobserved risk factors, possibly shared by clusters of subjects. The integration of multi‐state and frailty methodology is an interesting way to control for unobserved heterogeneity in the presence of complex event history structures and is particularly appealing for multicenter clinical trials. We propose the incorporation of correlated frailties in the transition‐specific hazard function, thanks to a nested hierarchy. We studied a semiparametric estimation approach based on maximum integrated partial likelihood. We show in a simulation study that the nested frailty multi‐state model improves the estimation of the effect of covariates, as well as the coverage probability of their confidence intervals. We present a case study concerning a prostate cancer multicenter clinical trial. The multi‐state nature of the model allows us to evidence the effect of treatment on death taking into account intermediate events. Copyright © 2015 JohnWiley & Sons, Ltd.     

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.     

Modelling two cause‐specific hazards of competing risks in one cumulative proportional odds model?

Competing risks extend standard survival analysis to considering time‐to‐first‐event and type‐of‐first‐event, where the event types are called competing risks. The competing risks process is completely described by all cause‐specific hazards, ie, the hazard marked by the event type. Separate Cox models for each cause‐specific hazard are the standard approach to regression modelling, but they come with the interpretational challenge that there are as many regression coefficients as there are competing risks. An alternative approach is to directly model the cumulative event probabilities, but again, there will be as many models as there are competing risks. The aim of this paper is to investigate the usefulness of a third alternative. Proportional odds modelling of all cause‐specific hazards summarizes the effect of one covariate on “opposing” competing outcomes in one regression coefficient. For instance, if the competing outcomes are hospital death and alive discharge from hospital, the modelling assumption is that a covariate affects both outcomes in opposing directions, but the effect size is of the same absolute magnitude. We will investigate the interpretational aspects of the approach analysing a data set on intensive care unit patients using parametric methods.     

Bayesian adaptive group lasso with semiparametric hidden Markov models

This paper presents a Bayesian adaptive group least absolute shrinkage and selection operator method to conduct simultaneous model selection and estimation under semiparametric hidden Markov models. We specify the conditional regression model and the transition probability model in the hidden Markov model into additive nonparametric functions of covariates. A basis expansion is adopted to approximate the nonparametric functions. We introduce multivariate conditional Laplace priors to impose adaptive penalties on regression coefficients and different groups of basis expansions under the Bayesian framework. An efficient Markov chain Monte Carlo algorithm is then proposed to identify the nonexistent, constant, linear, and nonlinear forms of covariate effects in both conditional and transition models. The empirical performance of the proposed methodology is evaluated via simulation studies. We apply the proposed model to analyze a real data set that was collected from the Alzheimer's Disease Neuroimaging Initiative study. The analysis identifies important risk factors on cognitive decline and the transition from cognitive normal to Alzheimer's disease.     

Competing‐risks duration models with correlated random effects: an application to dementia patients’ transition histories

Multi‐state transition models are widely applied tools to analyze individual event histories in the medical or social sciences. In this paper, we propose the use of (discrete‐time) competing‐risks duration models to analyze multi‐transition data. Unlike conventional Markov transition models, these models allow the estimated transition probabilities to depend on the time spent in the current state. Moreover, the models can be readily extended to allow for correlated transition probabilities. A further virtue of these models is that they can be estimated using conventional regression tools for discrete‐response data, such as the multinomial logit model. The latter is implemented in many statistical software packages and can be readily applied by empirical researchers. Moreover, model estimation is feasible, even when dealing with very large data sets, and simultaneously allowing for a flexible form of duration dependence and correlation between transition probabilities. We derive the likelihood function for a model with three competing target states and discuss a feasible and readily applicable estimation method. We also present the results from a simulation study, which indicate adequate performance of the proposed approach. In an empirical application, we analyze dementia patients’ transition probabilities from the domestic setting, taking into account several, partly duration‐dependent covariates. Copyright © 2014 John Wiley & Sons, Ltd.     

Spatially regularized estimation for the analysis of dynamic contrast‐enhanced magnetic resonance imaging data

Competing compartment models of different complexities have been used for the quantitative analysis of dynamic contrast‐enhanced magnetic resonance imaging data. We present a spatial elastic net approach that allows to estimate the number of compartments for each voxel such that the model complexity is not fixed a priori. A multi‐compartment approach is considered, which is translated into a restricted least square model selection problem. This is done by using a set of basis functions for a given set of candidate rate constants. The form of the basis functions is derived from a kinetic model and thus describes the contribution of a specific compartment. Using a spatial elastic net estimator, we chose a sparse set of basis functions per voxel, and hence, rate constants of compartments. The spatial penalty takes into account the voxel structure of an image and performs better than a penalty treating voxels independently. The proposed estimation method is evaluated for simulated images and applied to an in vivo dataset. Copyright © 2013 John Wiley & Sons, Ltd.     

Bayesian mixed hidden Markov models: a multi‐level approach to modeling categorical outcomes with differential misclassification

Questionnaire‐based health status outcomes are often prone to misclassification. When studying the effect of risk factors on such outcomes, ignoring any potential misclassification may lead to biased effect estimates. Analytical challenges posed by these misclassified outcomes are further complicated when simultaneously exploring factors for both the misclassification and health processes in a multi‐level setting. To address these challenges, we propose a fully Bayesian mixed hidden Markov model (BMHMM) for handling differential misclassification in categorical outcomes in a multi‐level setting. The BMHMM generalizes the traditional hidden Markov model (HMM) by introducing random effects into three sets of HMM parameters for joint estimation of the prevalence, transition, and misclassification probabilities. This formulation not only allows joint estimation of all three sets of parameters but also accounts for cluster‐level heterogeneity based on a multi‐level model structure. Using this novel approach, both the true health status prevalence and the transition probabilities between the health states during follow‐up are modeled as functions of covariates. The observed, possibly misclassified, health states are related to the true, but unobserved, health states and covariates. Results from simulation studies are presented to validate the estimation procedure, to show the computational efficiency due to the Bayesian approach and also to illustrate the gains from the proposed method compared to existing methods that ignore outcome misclassification and cluster‐level heterogeneity. We apply the proposed method to examine the risk factors for both asthma transition and misclassification in the Southern California Children's Health Study. Copyright © 2013 John Wiley & Sons, Ltd.     

Variance reduction in randomised trials by inverse probability weighting using the propensity score

In individually randomised controlled trials, adjustment for baseline characteristics is often undertaken to increase precision of the treatment effect estimate. This is usually performed using covariate adjustment in outcome regression models. An alternative method of adjustment is to use inverse probability‐of‐treatment weighting (IPTW), on the basis of estimated propensity scores. We calculate the large‐sample marginal variance of IPTW estimators of the mean difference for continuous outcomes, and risk difference, risk ratio or odds ratio for binary outcomes. We show that IPTW adjustment always increases the precision of the treatment effect estimate. For continuous outcomes, we demonstrate that the IPTW estimator has the same large‐sample marginal variance as the standard analysis of covariance estimator. However, ignoring the estimation of the propensity score in the calculation of the variance leads to the erroneous conclusion that the IPTW treatment effect estimator has the same variance as an unadjusted estimator; thus, it is important to use a variance estimator that correctly takes into account the estimation of the propensity score. The IPTW approach has particular advantages when estimating risk differences or risk ratios. In this case, non‐convergence of covariate‐adjusted outcome regression models frequently occurs. Such problems can be circumvented by using the IPTW adjustment approach. © 2013 The authors. Statistics in Medicine published by John Wiley & Sons, Ltd.     

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

This paper presents a parametric method of fitting semi‐Markov models with piecewise‐constant hazards in the presence of left, right and interval censoring. We investigate transition intensities in a three‐state illness–death model with no recovery. We relax the Markov assumption by adjusting the intensity for the transition from state 2 (illness) to state 3 (death) for the time spent in state 2 through a time‐varying covariate. This involves the exact time of the transition from state 1 (healthy) to state 2. When the data are subject to left or interval censoring, this time is unknown. In the estimation of the likelihood, we take into account interval censoring by integrating out all possible times for the transition from state 1 to state 2. For left censoring, we use an Expectation–Maximisation inspired algorithm. A simulation study reflects the performance of the method. The proposed combination of statistical procedures provides great flexibility. We illustrate the method in an application by using data on stroke onset for the older population from the UK Medical Research Council Cognitive Function and Ageing Study. Copyright © 2012 John Wiley & Sons, Ltd.     

Simultaneous estimation of parameters in the bivariate Emax model

In this paper, we explore inference in multi‐response, nonlinear models. By multi‐response, we mean models with m > 1 response variables and accordingly m relations. Each parameter/explanatory variable may appear in one or more of the relations. We study a system estimation approach for simultaneous computation and inference of the model and (co)variance parameters. For illustration, we fit a bivariate Emax model to diabetes dose‐response data. Further, the bivariate Emax model is used in a simulation study that compares the system estimation approach to equation‐by‐equation estimation. We conclude that overall, the system estimation approach performs better for the bivariate Emax model when there are dependencies among relations. The stronger the dependencies, the more we gain in precision by using system estimation rather than equation‐by‐equation estimation. Copyright © 2015 John Wiley & Sons, Ltd.     

Flexible modelling of discrete failure time including time-varying smooth effects

Discrete survival models have been extended in several ways. More flexible models are obtained by including time-varying coefficients and covariates which determine the hazard rate in an additive but not further specified form. In this paper, a general model is considered which comprises both types of covariate effects. An additional extension is the incorporation of smooth interaction between time and covariates. Thus, in the linear predictor smooth effects of covariates which may vary across time are allowed. It is shown how simple duration models produce artefacts which may be avoided by flexible models. For the general model which includes parametric terms, time-varying coefficients in parametric terms and time-varying smooth effects estimation procedures are derived which are based on the regularized expansion of smooth effects in basis functions. The approach is used to model the sojourn time in a psychiatric hospital. It is demonstrated how initial conditions which have non-linear influence are damped over time.     

Adjusting for time‐dependent sensitivity in an illness‐death model,with application to mother‐to‐child transmission of HIV

In mother‐to‐child transmission of HIV, identifying infected infants relies on a diagnostic test with imperfect sensitivity that is administered at scheduled visits. Under this scenario, a participant's true state may be unknown at the start and end times of the study, and the detection of transitions into illness may be delayed or missed altogether. This could lead to biased estimates of the risk of transmission and covariate associations. When a test has imperfect sensitivity, but perfect specificity, the additional uncertainty can be captured as a random variable measuring delay in detection. The cumulative distribution then defines a time‐dependent sensitivity function that increases over time. We present a maximum likelihood based illness‐death model that accounts for imperfect sensitivity by including the delay as an exponential distribution. We specify transition rates as penalized B‐splines to allow for nonhomogeneity of risk and discuss the model under Markov and semi‐Markov assumptions. We apply this method to our motivating data set, a study of 1499 mother and infant pairs at three sites in Africa. Copyright © 2014 John Wiley & Sons, Ltd.     

Fractional polynomial adjustment for time‐varying covariates in a self‐controlled case series analysis

The self‐controlled case series method is a statistical approach to investigating associations between acute outcomes and transient exposures. The method uses cases only and compares time at risk after the transient exposure with time at risk outside the exposure period within an individual, using conditional Poisson regression. The risk of outcome and exposure often varies over time, for example, with age, and it is important to allow for such time dependence within the analysis. The standard approach for modelling time‐varying covariates is to split observation periods into blocks according to categories of the covariate and then to model the relationship using indicators for each category. However, this can be inefficient and can lead to problems with collinearity if the exposure occurs at approximately the same time in all individuals. As an alternative, we propose using fractional polynomials to model the relationship between the time‐varying covariate and incidence of the outcome. We present the results from an analysis exploring the association between rotavirus vaccination and intussusception risk as well as a simulation study. We conclude that fractional polynomials provide a useful approach to adjusting for time‐varying covariates but that it is important to explore the sensitivity of the results to the number of categories and the method of adjustment. Copyright © 2013 John Wiley & Sons, Ltd.     

A case‐cohort approach for multi‐state models in hospital epidemiology

Analysing the determinants and consequences of hospital‐acquired infections involves the evaluation of large cohorts. Infected patients in the cohort are often rare for specific pathogens, because most of the patients admitted to the hospital are discharged or die without such an infection. Death and discharge are competing events to acquiring an infection, because these individuals are no longer at risk of getting a hospital‐acquired infection. Therefore, the data is best analysed with an extended survival model – the extended illness‐death model. A common problem in cohort studies is the costly collection of covariate values. In order to provide efficient use of data from infected as well as uninfected patients, we propose a tailored case‐cohort approach for the extended illness‐death model. The basic idea of the case‐cohort design is to only use a random sample of the full cohort, referred to as subcohort, and all cases, namely the infected patients. Thus, covariate values are only obtained for a small part of the full cohort. The method is based on existing and established methods and is used to perform regression analysis in adapted Cox proportional hazards models. We propose estimation of all cause‐specific cumulative hazards and transition probabilities in an extended illness‐death model based on case‐cohort sampling. As an example, we apply the methodology to infection with a specific pathogen using a large cohort from Spanish hospital data. The obtained results of the case‐cohort design are compared with the results in the full cohort to investigate the performance of the proposed method. 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.     

A population‐based approach to analyzing pulses in time series of hormone data

Studies of reproductive physiology involve rapid sampling protocols that result in time series of hormone concentrations. The signature pattern in these times series is pulses of hormone release. Various statistical models for quantifying the pulsatile release features exist. Currently these models are fitted separately to each individual and the resulting estimates averaged to arrive at post hoc population‐level estimates. When the signal‐to‐noise ratio is small or the time of observation is short (e.g., 6 h), this two‐stage estimation approach can fail. This work extends the single‐subject modelling framework to a population framework similar to what exists for complex pharamacokinetics data. The goal is to leverage information across subjects to more clearly identify pulse locations and improve estimation of other model parameters. This modelling extension has proven difficult because the pulse number and locations are unknown. Here, we show that simultaneously modelling a group of subjects is computationally feasible in a Bayesian framework using a birth–death Markov chain Monte Carlo estimation algorithm. Via simulation, we show that this population‐based approach reduces the false positive and negative pulse detection rates and results in less biased estimates of population‐level parameters of frequency, pulse size, and hormone elimination. We then apply the approach to a reproductive study in healthy women where approximately one‐third of the 21 subjects in the study did not have appropriate fits using the single‐subject fitting approach. Using the population model produced more precise, biologically plausible estimates of all model parameters. Copyright © 2017 John Wiley & Sons, Ltd.