Model selection for Cox models with time-varying coefficients

Summary Cox models with time-varying coefficients offer great flexibility in capturing the temporal dynamics of covariate effects on right-censored failure times. Because not all covariate coefficients are time varying, model selection for such models presents an additional challenge, which is to distinguish covariates with time-varying coefficient from those with time-independent coefficient. We propose an adaptive group lasso method that not only selects important variables but also selects between time-independent and time-varying specifications of their presence in the model. Each covariate effect is partitioned into a time-independent part and a time-varying part, the latter of which is characterized by a group of coefficients of basis splines without intercept. Model selection and estimation are carried out through a fast, iterative group shooting algorithm. Our approach is shown to have good properties in a simulation study that mimics realistic situations with up to 20 variables. A real example illustrates the utility of the method.     

Semiparametric approaches for joint modeling of longitudinal and survival data with time-varying coefficients

Summary .   We study joint modeling of survival and longitudinal data. There are two regression models of interest. The primary model is for survival outcomes, which are assumed to follow a time-varying coefficient proportional hazards model. The second model is for longitudinal data, which are assumed to follow a random effects model. Based on the trajectory of a subject's longitudinal data, some covariates in the survival model are functions of the unobserved random effects. Estimated random effects are generally different from the unobserved random effects and hence this leads to covariate measurement error. To deal with covariate measurement error, we propose a local corrected score estimator and a local conditional score estimator. Both approaches are semiparametric methods in the sense that there is no distributional assumption needed for the underlying true covariates. The estimators are shown to be consistent and asymptotically normal. However, simulation studies indicate that the conditional score estimator outperforms the corrected score estimator for finite samples, especially in the case of relatively large measurement error. The approaches are demonstrated by an application to data from an HIV clinical trial.     

Estimation for the optimal combination of markers without modeling the censoring distribution

Summary .  In the time-dependent receiver operating characteristic curve analysis with several baseline markers, research interest focuses on seeking appropriate composite markers to enhance the accuracy in predicting the vital status of individuals over time. Based on censored survival data, we proposed a more flexible estimation procedure for the optimal combination of markers under the validity of a time-varying coefficient generalized linear model for the event time without restrictive assumptions on the censoring pattern. The consistency of the proposed estimators is also established in this article. In contrast, the inverse probability weighting (IPW) approach might introduce a bias when the selection probabilities are misspecified in the estimating equations. The performance of both estimation procedures are examined and compared through a class of simulations. It is found from the simulation study that the proposed estimators are far superior to the IPW ones. Applying these methods to an angiography cohort, our estimation procedure is shown to be useful in predicting the time to all-cause and coronary artery disease related death.     

Analyzing complex capture-recapture data in the presence of individual and temporal covariates and model uncertainty

SUMMARY: We consider the issue of analyzing complex ecological data in the presence of covariate information and model uncertainty. Several issues can arise when analyzing such data, not least the need to take into account where there are missing covariate values. This is most acutely observed in the presence of time-varying covariates. We consider mark-recapture-recovery data, where the corresponding recapture probabilities are less than unity, so that individuals are not always observed at each capture event. This often leads to a large amount of missing time-varying individual covariate information, because the covariate cannot usually be recorded if an individual is not observed. In addition, we address the problem of model selection over these covariates with missing data. We consider a Bayesian approach, where we are able to deal with large amounts of missing data, by essentially treating the missing values as auxiliary variables. This approach also allows a quantitative comparison of different models via posterior model probabilities, obtained via the reversible jump Markov chain Monte Carlo algorithm. To demonstrate this approach we analyze data relating to Soay sheep, which pose several statistical challenges in fully describing the intricacies of the system.     

On the relation between the cause-specific hazard and the subdistribution rate for competing risks data: The Fine–Gray model revisited

The Fine–Gray proportional subdistribution hazards model has been puzzling many people since its introduction. The main reason for the uneasy feeling is that the approach considers individuals still at risk for an event of cause 1 after they fell victim to the competing risk of cause 2. The subdistribution hazard and the extended risk sets, where subjects who failed of the competing risk remain in the risk set, are generally perceived as unnatural . One could say it is somewhat of a riddle why the Fine–Gray approach yields valid inference. To take away these uneasy feelings, we explore the link between the Fine–Gray and cause-specific approaches in more detail. We introduce the reduction factor as representing the proportion of subjects in the Fine–Gray risk set that has not yet experienced a competing event. In the presence of covariates, the dependence of the reduction factor on a covariate gives information on how the effect of the covariate on the cause-specific hazard and the subdistribution hazard relate. We discuss estimation and modeling of the reduction factor, and show how they can be used in various ways to estimate cumulative incidences, given the covariates. Methods are illustrated on data of the European Society for Blood and Marrow Transplantation.     

Applying the Cox proportional hazards model when the change time of a binary time-varying covariate is interval censored

This paper develops methodology for estimation of the effect of a binary time-varying covariate on failure times when the change time of the covariate is interval censored. The motivating example is a study of cytomegalovirus (CMV) disease in patients with human immunodeficiency virus (HIV) disease. We are interested in determining whether CMV shedding predicts an increased hazard for developing active CMV disease. Since a clinical screening test is needed to detect CMV shedding, the time that shedding begins is only known to lie in an interval bounded by the patient's last negative and first positive tests. In a Cox proportional hazards model with a time-varying covariate for CMV shedding, the partial likelihood depends on the covariate status of every individual in the risk set at each failure time. Due to interval censoring, this is not always known. To solve this problem, we use a Monte Carlo EM algorithm with a Gibbs sampler embedded in the E-step. We generate multiple completed data sets by drawing imputed exact shedding times based on the joint likelihood of the shedding times and event times under the Cox model. The method is evaluated using a simulation study and is applied to the data set described above.     

Adjustment for Missingness Using Auxiliary Information in Semiparametric Regression

Summary .  In this article, we study the estimation of mean response and regression coefficient in semiparametric regression problems when response variable is subject to nonrandom missingness. When the missingness is independent of the response conditional on high-dimensional auxiliary information, the parametric approach may misspecify the relationship between covariates and response while the nonparametric approach is infeasible because of the curse of dimensionality. To overcome this, we study a model-based approach to condense the auxiliary information and estimate the parameters of interest nonparametrically on the condensed covariate space. Our estimators possess the double robustness property, i.e., they are consistent whenever the model for the response given auxiliary covariates or the model for the missingness given auxiliary covariate is correct. We conduct a number of simulations to compare the numerical performance between our estimators and other existing estimators in the current missing data literature, including the propensity score approach and the inverse probability weighted estimating equation. A set of real data is used to illustrate our approach.     

Modeling the hazard of transition into the absorbing state in the illness-death model

The illness-death model is the simplest multistate model where the transition from the initial state 0 to the absorbing state 2 may involve an intermediate state 1 (e.g., disease relapse). The impact of the transition into state 1 on the subsequent transition hazard to state 2 enables insight to be gained into the disease evolution. The standard approach of analysis is modeling the transition hazards from 0 to 2 and from 1 to 2, including time to illness as a time-varying covariate and measuring time from origin even after transition into state 1. The hazard from 1 to 2 can be also modeled separately using only patients in state 1, measuring time from illness and including time to illness as a fixed covariate. A recently proposed approach is a model where time after the transition into state 1 is measured in both scales and time to illness is included as a time-varying covariate. Another possibility is a model where time after transition into state 1 is measured only from illness and time to illness is included as a fixed covariate. Through theoretical reasoning and simulation protocols, we discuss the use of these models and we develop a practical strategy aiming to (a) validate the properties of the illness-death process, (b) estimate the impact of time to illness on the hazard from state 1 to 2, and (c) quantify the impact that the transition into state 1 has on the hazard of the absorbing state. The strategy is also applied to a literature dataset on diabetes.     

A censored quantile regression approach for relative survival analysis: Relative survival quantile regression

We propose a censored quantile regression model for the analysis of relative survival data. We create a hybrid data set consisting of the study observations and counterpart randomly sampled pseudopopulation observations imputed from population life tables that adjust for expected mortality. We then fit a censored quantile regression model to the hybrid data incorporating demographic variables (e.g., age, biologic sex, calendar time) corresponding to the population life tables of demographically-similar individuals, a population versus study covariate, and its interactions with the variables of interest. These latter variables can be interpreted as relative survival parameters that depict the differences in failure quantiles between the study participants and their population counterparts.     

Modelling paired survival data with covariates

The objective of this paper is to consider the parametric analysis of paired censored survival data when additional covariate information is available, as in the Diabetic Retinopathy Study, which assessed the effectiveness of laser photocoagulation in delaying loss of visual acuity. Our first approach is to extend the fully parametric model of Clayton (1978, Biometrika 65, 141-151) to incorporate covariate information. Our second approach is to obtain parameter estimates from an independence working model together with robust variance estimates. The approaches are compared in terms of efficiency and computational considerations. A fundamental consideration in choosing a strategy for the analysis of paired survival data is whether the correlation within a pair is a nuisance parameter or a parameter of intrinsic scientific interest. The approaches are illustrated with the Diabetic Retinopathy Study.     

Comparison of procedures to assess non-linear and time-varying effects in multivariable models for survival data

The focus of many medical applications is to model the impact of several factors on time to an event. A standard approach for such analyses is the Cox proportional hazards model. It assumes that the factors act linearly on the log hazard function (linearity assumption) and that their effects are constant over time (proportional hazards (PH) assumption). Variable selection is often required to specify a more parsimonious model aiming to include only variables with an influence on the outcome. As follow-up increases the effect of a variable often gets weaker, which means that it varies in time. However, spurious time-varying effects may also be introduced by mismodelling other parts of the multivariable model, such as omission of an important covariate or an incorrect functional form of a continuous covariate. These issues interact. To check whether the effect of a variable varies in time several tests for non-PH have been proposed. However, they are not sufficient to derive a model, as appropriate modelling of the shape of time-varying effects is required. In three examples we will compare five recently published strategies to assess whether and how the effects of covariates from a multivariable model vary in time. For practical use we will give some recommendations.     

Bayesian inference on age-specific survival for censored and truncated data

1. Traditional estimation of age-specific survival and mortality rates in vertebrates is limited to individuals with known age. Although this subject has been studied extensively using effective capture-recapture and capture-recovery models, inference remains challenging because of large numbers of incomplete records (i.e. unknown age of many individuals) and because of the inadequate duration of the studies. 2. Here, we present a hierarchical model for capture-recapture/recovery (CRR) data sets with large proportions of unknown times of birth and death. The model uses a Bayesian framework to draw inference on population-level age-specific demographic rates using parametric survival functions and applies this information to reconstruct times of birth and death for individuals with unknown age. 3. We simulated a set of CRR data sets with varying study span and proportions of individuals with known age, and varying recapture and recovery probabilities. We used these data sets to compare our method to a traditional CRR model, which requires knowledge of individual ages. Subsequently, we applied our method to a subset of a long-term CRR data set on Soay sheep. 4. Our results show that this method performs better than the common CRR model when sample sizes are low. Still, our model is sensitive to the choice of priors with low recapture probability and short studies. In such cases, priors that overestimate survival perform better than those that underestimate it. Also, the model was able to estimate accurately ages at death for Soay sheep, with an average error of 0.94 years and to identify differences in mortality rate between sexes. 5. Although many of the problems in the estimation of age-specific survival can be reduced through more efficient sampling schemes, most ecological data sets are still sparse and with a large proportion of missing records. Thus, improved sampling needs still to be combined with statistical models capable of overcoming the unavoidable limitations of any fieldwork. We show that our approach provides reliable estimates of parameters and unknown times of birth and death even with the most incomplete data sets while being flexible enough to accommodate multiple recapture probabilities and covariates.     

Graphical exploration of covariate effects on survival data through nonparametric quantile curves

Kaplan-Meier curves provide an effective means of presenting the distributional pattern in a sample of survival data. However, in order to assess the effect of a covariate, a standard scatterplot is often difficult to interpret because of the presence of censored observations. Several authors have proposed a running median as an effective way of indicating the effect of a covariate. This article proposes a form of kernel estimation, employing double smoothing, that can be applied in a simple and efficient manner to construct an estimator of a percentile of the survival distribution as a function of one or two covariates. Permutations and bootstrap samples can be used to construct reference bands that help identify whether particular features of the estimates indicate real features of the underlying curve or whether this may be due simply to random variation. The techniques are illustrated on data from a study of kidney transplant patients.     

Estimation in the Cox survival regression model with covariate measurement error and a changepoint

The Cox regression model is a popular model for analyzing the relationship between a covariate vector and a survival endpoint. The standard Cox model assumes a constant covariate effect across the entire covariate domain. However, in many epidemiological and other applications, the covariate of main interest is subject to a threshold effect: a change in the slope at a certain point within the covariate domain. Often, the covariate of interest is subject to some degree of measurement error. In this paper, we study measurement error correction in the case where the threshold is known. Several bias correction methods are examined: two versions of regression calibration (RC1 and RC2, the latter of which is new), two methods based on the induced relative risk under a rare event assumption (RR1 and RR2, the latter of which is new), a maximum pseudo-partial likelihood estimator (MPPLE), and simulation-extrapolation (SIMEX). We develop the theory, present simulations comparing the methods, and illustrate their use on data concerning the relationship between chronic air pollution exposure to particulate matter PM₁₀ and fatal myocardial infarction (Nurses Health Study (NHS)), and on data concerning the effect of a subject's long-term underlying systolic blood pressure level on the risk of cardiovascular disease death (Framingham Heart Study (FHS)). The simulations indicate that the best methods are RR2 and MPPLE.     

Inverse probability weighting methods for Cox regression with right-truncated data

Right-truncated data arise when observations are ascertained retrospectively, and only subjects who experience the event of interest by the time of sampling are selected. Such a selection scheme, without adjustment, leads to biased estimation of covariate effects in the Cox proportional hazards model. The existing methods for fitting the Cox model to right-truncated data, which are based on the maximization of the likelihood or solving estimating equations with respect to both the baseline hazard function and the covariate effects, are numerically challenging. We consider two alternative simple methods based on inverse probability weighting (IPW) estimating equations, which allow consistent estimation of covariate effects under a positivity assumption and avoid estimation of baseline hazards. We discuss problems of identifiability and consistency that arise when positivity does not hold and show that although the partial tests for null effects based on these IPW methods can be used in some settings even in the absence of positivity, they are not valid in general. We propose adjusted estimating equations that incorporate the probability of observation when it is known from external sources, which results in consistent estimation. We compare the methods in simulations and apply them to the analyses of human immunodeficiency virus latency.     

Repeated probit regression when covariates are measured with error

This paper develops a model for repeated binary regression when a covariate is measured with error. The model allows for estimating the effect of the true value of the covariate on a repeated binary response. The choice of a probit link for the effect of the error-free covariate, coupled with normal measurement error for the error-free covariate, results in a probit model after integrating over the measurement error distribution. We propose a two-stage estimation procedure where, in the first stage, a linear mixed model is used to fit the repeated covariate. In the second stage, a model for the correlated binary responses conditional on the linear mixed model estimates is fit to the repeated binary data using generalized estimating equations. The approach is demonstrated using nutrient safety data from the Diet Intervention of School Age Children (DISC) study.     

Estimation and interpretation of incidence rate difference for recurrent events when the estimation model is misspecified

Recurrent events data are common in experimental and observational studies. It is often of interest to estimate the effect of an intervention on the incidence rate of the recurrent events. The incidence rate difference is a useful measure of intervention effect. A weighted least squares estimator of the incidence rate difference for recurrent events was recently proposed for an additive rate model in which both the baseline incidence rate and the covariate effects were constant over time. In this article, we relax this model assumption and examine the properties of the estimator under the additive and multiplicative rate models assumption in which the baseline incidence rate and covariate effects may vary over time. We show analytically and numerically that the estimator gives an appropriate summary measure of the time‐varying covariate effects. In particular, when the underlying covariate effects are additive and time‐varying, the estimator consistently estimates the weighted average of the covariate effects over time. When the underlying covariate effects are multiplicative and time‐varying, and if there is only one binary covariate indicating the intervention status, the estimator consistently estimates the weighted average of the underlying incidence rate difference between the intervention and control groups over time. We illustrate the method with data from a randomized vaccine trial.     

Modeling Systematic Change in Stopover Duration Does Not Improve Bias in Trends Estimated from Migration Counts

The use of counts of unmarked migrating animals to monitor long term population trends assumes independence of daily counts and a constant rate of detection. However, migratory stopovers often last days or weeks, violating the assumption of count independence. Further, a systematic change in stopover duration will result in a change in the probability of detecting individuals once, but also in the probability of detecting individuals on more than one sampling occasion. We tested how variation in stopover duration influenced accuracy and precision of population trends by simulating migration count data with known constant rate of population change and by allowing daily probability of survival (an index of stopover duration) to remain constant, or to vary randomly, cyclically, or increase linearly over time by various levels. Using simulated datasets with a systematic increase in stopover duration, we also tested whether any resulting bias in population trend could be reduced by modeling the underlying source of variation in detection, or by subsampling data to every three or five days to reduce the incidence of recounting. Mean bias in population trend did not differ significantly from zero when stopover duration remained constant or varied randomly over time, but bias and the detection of false trends increased significantly with a systematic increase in stopover duration. Importantly, an increase in stopover duration over time resulted in a compounding effect on counts due to the increased probability of detection and of recounting on subsequent sampling occasions. Under this scenario, bias in population trend could not be modeled using a covariate for stopover duration alone. Rather, to improve inference drawn about long term population change using counts of unmarked migrants, analyses must include a covariate for stopover duration, as well as incorporate sampling modifications (e.g., subsampling) to reduce the probability that individuals will be detected on more than one occasion.     

A flexible cure rate model for spatially correlated survival data based on generalized extreme value distribution and Gaussian process priors

Our present work proposes a new survival model in a Bayesian context to analyze right‐censored survival data for populations with a surviving fraction, assuming that the log failure time follows a generalized extreme value distribution. Many applications require a more flexible modeling of covariate information than a simple linear or parametric form for all covariate effects. It is also necessary to include the spatial variation in the model, since it is sometimes unexplained by the covariates considered in the analysis. Therefore, the nonlinear covariate effects and the spatial effects are incorporated into the systematic component of our model. Gaussian processes (GPs) provide a natural framework for modeling potentially nonlinear relationship and have recently become extremely powerful in nonlinear regression. Our proposed model adopts a semiparametric Bayesian approach by imposing a GP prior on the nonlinear structure of continuous covariate. With the consideration of data availability and computational complexity, the conditionally autoregressive distribution is placed on the region‐specific frailties to handle spatial correlation. The flexibility and gains of our proposed model are illustrated through analyses of simulated data examples as well as a dataset involving a colon cancer clinical trial from the state of Iowa.