Bias and Sensitivity Analysis When Estimating Treatment Effects from the Cox Model with Omitted Covariates

Summary

Omission of relevant covariates can lead to bias when estimating treatment or exposure effects from survival data in both randomized controlled trials and observational studies. This paper presents a general approach to assessing bias when covariates are omitted from the Cox model. The proposed method is applicable to both randomized and non‐randomized studies. We distinguish between the effects of three possible sources of bias: omission of a balanced covariate, data censoring and unmeasured confounding. Asymptotic formulae for determining the bias are derived from the large sample properties of the maximum likelihood estimator. A simulation study is used to demonstrate the validity of the bias formulae and to characterize the influence of the different sources of bias. It is shown that the bias converges to fixed limits as the effect of the omitted covariate increases, irrespective of the degree of confounding. The bias formulae are used as the basis for developing a new method of sensitivity analysis to assess the impact of omitted covariates on estimates of treatment or exposure effects. In simulation studies, the proposed method gave unbiased treatment estimates and confidence intervals with good coverage when the true sensitivity parameters were known. We describe application of the method to a randomized controlled trial and a non‐randomized study.     

The effect of treatment delay on time-to-recovery in the presence of unobserved heterogeneity

We study the effect of delaying treatment in the presence of (unobserved) heterogeneity. In a homogeneous population and assuming a proportional treatment effect, a treatment delay period will result in notably lower cumulative recovery percentages. We show in theoretical scenarios using frailty models that if the population is heterogeneous, the effect of a delay period is much smaller. This can be explained by the selection process that is induced by the frailty. Patient groups that start treatment later have already undergone more selection. The marginal hazard ratio for the treatment will act differently in such a more homogeneous patient group. We further discuss modeling approaches for estimating the effect of treatment delay in the presence of heterogeneity, and compare their performance in a simulation study. The conventional Cox model that fails to account for heterogeneity overestimates the effect of treatment delay. Including interaction terms between treatment and starting time of treatment or between treatment and follow up time gave no improvement. Estimating a frailty term can improve the estimation, but is sensitive to misspecification of the frailty distribution. Therefore, multiple frailty distributions should be used and the results should be compared using the Akaike Information Criterion. Non-parametric estimation of the cumulative recovery percentages can be considered if the dataset contains sufficient long term follow up for each of the delay strategies. The methods are demonstrated on a motivating application evaluating the effect of delaying the start of treatment with assisted reproductive techniques on time-to-pregnancy in couples with unexplained subfertility.     

Marginal hazard ratio estimates in joint frailty models for heart failure trials

This work is motivated by clinical trials in chronic heart failure disease, where treatment has effects both on morbidity (assessed as recurrent non‐fatal hospitalisations) and on mortality (assessed as cardiovascular death, CV death). Recently, a joint frailty proportional hazards model has been proposed for these kind of efficacy outcomes to account for a potential association between the risk rates for hospital admissions and CV death. However, more often clinical trial results are presented by treatment effect estimates that have been derived from marginal proportional hazards models, that is, a Cox model for mortality and an Andersen–Gill model for recurrent hospitalisations. We show how these marginal hazard ratios and their estimates depend on the association between the risk processes, when these are actually linked by shared or dependent frailty terms. First we derive the marginal hazard ratios as a function of time. Then, applying least false parameter theory, we show that the marginal hazard ratio estimate for the hospitalisation rate depends on study duration and on parameters of the underlying joint frailty model. In particular, we identify parameters, for example the treatment effect on mortality, that determine if the marginal hazard ratio estimate for hospitalisations is smaller, equal or larger than the conditional one. How this affects rejection probabilities is further investigated in simulation studies. Our findings can be used to interpret marginal hazard ratio estimates in heart failure trials and are illustrated by the results of the CHARM‐Preserved trial (where CHARM is the ‘Candesartan in Heart failure Assessment of Reduction in Mortality and morbidity’ programme).     

A causal proportional hazards estimator for the effect of treatment actually received in a randomized trial with all-or-nothing compliance

Survival data from randomized trials are most often analyzed in a proportional hazards (PH) framework that follows the intention-to-treat (ITT) principle. When not all the patients on the experimental arm actually receive the assigned treatment, the ITT-estimator mixes its effect on treatment compliers with its absence of effect on noncompliers. The structural accelerated failure time (SAFT) models of Robins and Tsiatis are designed to consistently estimate causal effects on the treated, without direct assumptions about the compliance selection mechanism. The traditional PH-model, however, has not yet led to such causal interpretation. In this article, we examine a PH-model of treatment effect on the treated subgroup. While potential treatment compliance is unobserved in the control arm, we derive an estimating equation for the Compliers PROPortional Hazards Effect of Treatment (C-PROPHET). The jackknife is used for bias correction and variance estimation. The method is applied to data from a recently finished clinical trial in cancer patients with liver metastases.     

Competing risks regression for stratified data

For competing risks data, the Fine-Gray proportional hazards model for subdistribution has gained popularity for its convenience in directly assessing the effect of covariates on the cumulative incidence function. However, in many important applications, proportional hazards may not be satisfied, including multicenter clinical trials, where the baseline subdistribution hazards may not be common due to varying patient populations. In this article, we consider a stratified competing risks regression, to allow the baseline hazard to vary across levels of the stratification covariate. According to the relative size of the number of strata and strata sizes, two stratification regimes are considered. Using partial likelihood and weighting techniques, we obtain consistent estimators of regression parameters. The corresponding asymptotic properties and resulting inferences are provided for the two regimes separately. Data from a breast cancer clinical trial and from a bone marrow transplantation registry illustrate the potential utility of the stratified Fine-Gray model.     

Regression models for infant mortality data in Norwegian siblings, using a compound Poisson frailty distribution with random scale

The power variance function distributions, which include the gamma and compound Poisson (CP) distributions among others, are commonly used in frailty models for family data. In a previous paper, we presented a frailty model constructed by randomizing the scale parameter in a CP distribution. When combined with a parametric baseline hazard, this yields a model with heterogeneity on both the individual and the family level and a subgroup with zero frailty, corresponding to people not experiencing the event. In this paper, we discuss covariates in the model. Depending on where the covariates are inserted in the model, one may have proportional hazards at the individual level, the family level, and a larger group level (for covariates shared by many families, e.g. ethnic groups) or get accelerated failure times. Each of these alternatives gives a specific interpretation of the covariate effects. An application to data infant mortality in siblings from the Medical Birth Registry of Norway is included. We compare the results for some of the different covariate modeling options.     

Bias correction in the hierarchical likelihood approach to the analysis of multivariate survival data

Frailty models are useful for measuring unobserved heterogeneity in risk of failures across clusters, providing cluster-specific risk prediction. In a frailty model, the latent frailties shared by members within a cluster are assumed to act multiplicatively on the hazard function. In order to obtain parameter and frailty variate estimates, we consider the hierarchical likelihood (H-likelihood) approach (Ha, Lee and Song, 2001. Hierarchical-likelihood approach for frailty models. Biometrika 88, 233-243) in which the latent frailties are treated as "parameters" and estimated jointly with other parameters of interest. We find that the H-likelihood estimators perform well when the censoring rate is low, however, they are substantially biased when the censoring rate is moderate to high. In this paper, we propose a simple and easy-to-implement bias correction method for the H-likelihood estimators under a shared frailty model. We also extend the method to a multivariate frailty model, which incorporates complex dependence structure within clusters. We conduct an extensive simulation study and show that the proposed approach performs very well for censoring rates as high as 80%. We also illustrate the method with a breast cancer data set. Since the H-likelihood is the same as the penalized likelihood function, the proposed bias correction method is also applicable to the penalized likelihood estimators.     

A General Class of Semiparametric Transformation Frailty Models for Nonproportional Hazards Survival Data

Summary We propose a semiparametrically efficient estimation of a broad class of transformation regression models for nonproportional hazards data. Classical transformation models are to be viewed from a frailty model paradigm, and the proposed method provides a unified approach that is valid for both continuous and discrete frailty models. The proposed models are shown to be flexible enough to model long‐term follow‐up survival data when the treatment effect diminishes over time, a case for which the PH or proportional odds assumption is violated, or a situation in which a substantial proportion of patients remains cured after treatment. Estimation of the link parameter in frailty distribution, considered to be unknown and possibly dependent on a time‐independent covariates, is automatically included in the proposed methods. The observed information matrix is computed to evaluate the variances of all the parameter estimates. Our likelihood‐based approach provides a natural way to construct simple statistics for testing the PH and proportional odds assumptions for usual survival data or testing the short‐ and long‐term effects for survival data with a cure fraction. Simulation studies demonstrate that the proposed inference procedures perform well in realistic settings. Applications to two medical studies are provided.     

Incorporation of inter-individual heterogeneity into the multistage carcinogenesis model: approach to the analysis of cancer incidence data

We investigate a multistage carcinogenesis frailty model to incorporate inter-individual heterogeneity into carcinogenic response. Attention is focused on inference concerning the effects of different sources of population heterogeneity on cancer rates. The authors consider unobserved variability arising from either carcinogen exposure or background characteristics. Gamma and Inverse-Gaussian distributions are selected for frailty models, and the baseline hazard function is the generalized Armitage-Doll model (i.e. non-frailty model) in which exposure effects shift the age scale instead of acting multiplicatively on cancer rates. For illustration, we apply the method to solid cancer data from a cohort of atomic bomb survivors to examine some features of proposed models. The results show that the Gamma frailty model for the heterogeneity of baseline rates provides the best goodness-of-fit of the model and a non-zero frailty variance. Parameter estimates are, for the most part, comparable between the Gamma and Inverse-Gaussian frailty models. In a heterogeneous population the exposure effects on young adulthood cancer rates might be underestimated for the non-frailty model. Meaningful information regarding each source of heterogeneity has been provided by the proposed method. Therefore, the multistage carcinogenesis frailty model approach is useful for analyses of epidemiological cancer data to assess population heterogeneity and heterogeneity-influenced exposure effects.     

Bayesian Semiparametric Frailty Selection in Multivariate Event Time Data

Biomedical studies often collect multivariate event time data from multiple clusters (either subjects or groups) within each of which event times for individuals are correlated and the correlation may vary in different classes. In such survival analyses, heterogeneity among clusters for shared and specific classes can be accommodated by incorporating parametric frailty terms into the model. In this article, we propose a Bayesian approach to relax the parametric distribution assumption for shared and specific‐class frailties by using a Dirichlet process prior while also allowing for the uncertainty of heterogeneity for different classes. Multiple cluster‐specific frailty selections rely on variable selection‐type mixture priors by applying mixtures of point masses at zero and inverse gamma distributions to the variance of log frailties. This selection allows frailties with zero variance to effectively drop out of the model. A reparameterization of log‐frailty terms is performed to reduce the potential bias of fixed effects due to variation of the random distribution and dependence among the parameters resulting in easy interpretation and faster Markov chain Monte Carlo convergence. Simulated data examples and an application to a lung cancer clinical trial are used for illustration.     

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.     

A sensitivity analysis approach for the causal hazard ratio in randomized and observational studies

The hazard ratio (HR) is often reported as the main causal effect when studying survival data. Despite its popularity, the HR suffers from an unclear causal interpretation. As already pointed out in the literature, there is a built-in selection bias in the HR, because similarly to the truncation by death problem, the HR conditions on post-treatment survival. A recently proposed alternative, inspired by the Survivor Average Causal Effect, is the causal HR, defined as the ratio between hazards across treatment groups among the study participants that would have survived regardless of their treatment assignment. We discuss the challenge in identifying the causal HR and present a sensitivity analysis identification approach in randomized controlled trials utilizing a working frailty model. We further extend our framework to adjust for potential confounders using inverse probability of treatment weighting. We present a Cox-based and a flexible non-parametric kernel-based estimation under right censoring. We study the finite-sample properties of the proposed estimation methods through simulations. We illustrate the utility of our framework using two real-data examples.     

Periodontal Treatment for Preventing Adverse Pregnancy Outcomes: A Meta- and Trial Sequential Analysis

Objectives

Periodontal treatment might reduce adverse pregnancy outcomes. The efficacy of periodontal treatment to prevent preterm birth, low birth weight, and perinatal mortality was evaluated using meta-analysis and trial sequential analysis.

Methods

An existing systematic review was updated and meta-analyses performed. Risk of bias, heterogeneity, and publication bias were evaluated, and meta-regression performed. Subgroup analysis was used to compare different studies with low and high risk of bias and different populations, i.e., risk groups. Trial sequential analysis was used to assess risk of random errors.

Results

Thirteen randomized clinical trials evaluating 6283 pregnant women were meta-analyzed. Four and nine trials had low and high risk of bias, respectively. Overall, periodontal treatment had no significant effect on preterm birth (odds ratio [95% confidence interval] 0.79 [0.57-1.10]) or low birth weight (0.69 [0.43-1.13]). Trial sequential analysis demonstrated that futility was not reached for any of the outcomes. For populations with moderate occurrence (<20%) of preterm birth or low birth weight, periodontal treatment was not efficacious for any of the outcomes, and trial sequential analyses indicated that further trials might be futile. For populations with high occurrence (≥20%) of preterm birth and low birth weight, periodontal treatment seemed to reduce the risk of preterm birth (0.42 [0.24-0.73]) and low birth weight (0.32 [0.15-0.67]), but trial sequential analyses showed that firm evidence was not reached. Periodontal treatment did not significantly affect perinatal mortality, and firm evidence was not reached. Risk of bias, but not publication bias or patients’ age modified the effect estimates.

Conclusions

Providing periodontal treatment to pregnant women could potentially reduce the risks of perinatal outcomes, especially in mothers with high risks. Conclusive evidence could not be reached due to risks of bias, risks of random errors, and unclear effects of confounding. Further randomized clinical trials are required.     

Selecting factors predictive of heterogeneity in multivariate event time data

In multivariate survival analysis, investigators are often interested in testing for heterogeneity among clusters, both overall and within specific classes. We represent different hypotheses about the heterogeneity structure using a sequence of gamma frailty models, ranging from a null model with no random effects to a full model having random effects for each class. Following a Bayesian approach, we define prior distributions for the frailty variances consisting of mixtures of point masses at zero and inverse-gamma densities. Since frailties with zero variance effectively drop out of the model, this prior allocates probability to each model in the sequence, including the overall null hypothesis of homogeneity. Using a counting process formulation, the conditional posterior distributions of the frailties and proportional hazards regression coefficients have simple forms. Posterior computation proceeds via a data augmentation Gibbs sampling algorithm, a single run of which can be used to obtain model-averaged estimates of the population parameters and posterior model probabilities for testing hypotheses about the heterogeneity structure. The methods are illustrated using data from a lung cancer trial.     

Estimating haplotype relative risks on human survival in population-based association studies

Association-based linkage disequilibrium (LD) mapping is an increasingly important tool for localizing genes that show potential influence on human aging and longevity. As haplotypes contain more LD information than single markers, a haplotype-based LD approach can have increased power in detecting associations as well as increased robustness in statistical testing. In this paper, we develop a new statistical model to estimate haplotype relative risks (HRRs) on human survival using unphased multilocus genotype data from unrelated individuals in cross-sectional studies. Based on the proportional hazard assumption, the model can estimate haplotype risk and frequency parameters, incorporate observed covariates, assess interactions between haplotypes and the covariates, and investigate the modes of gene function. By introducing population survival information available from population statistics, we are able to develop a procedure that carries out the parameter estimation using a nonparametric baseline hazard function and estimates sex-specific HRRs to infer gene-sex interaction. We also evaluate the haplotype effects on human survival while taking into account individual heterogeneity in the unobserved genetic and nongenetic factors or frailty by introducing the gamma-distributed frailty into the survival function. After model validation by computer simulation, we apply our method to an empirical data set to measure haplotype effects on human survival and to estimate haplotype frequencies at birth and over the observed ages. Results from both simulation and model application indicate that our survival analysis model is an efficient method for inferring haplotype effects on human survival in population-based association studies.     

The use of frailty hazard models for unrecognized heterogeneity that interacts with treatment: considerations of efficiency and power

Increasingly, genetic studies of tumors of the same histologic diagnosis are elucidating subtypes that are distinct with respect to clinical endpoints such as response to treatment and survival. This raises concerns about the efficiency of using the simple log-rank test for analysis of treatment effect on survival in studies of possibly heterogeneous tumors. Furthermore, such studies, designed under the assumption of homogeneity, may be severely underpowered. We derive analytic approximations for the asymptotic relative efficiency of the simple log-rank test relative to the optimally weighted log-rank test and for the power of the simple log-rank test when applied to subjects with unobserved heterogeneity, as reflected in a continuous frailty, that may interact with treatment. Numerical studies demonstrate that the simple log-rank test may be quite inefficient if the frailty interacts with treatment. Further, there may be a substantial loss of power in the presence of the frailty with or without an interaction with treatment.     

Merits of Modelling Multivariate Survival Data Using Random Effects Proportional Odds Model

Recently, there has been a great deal of interest in the analysis of multivariate survival data. In most epidemiological studies, survival times of the same cluster are related because of some unobserved risk factors such as the environmental or genetic factors. Therefore, modelling of dependence between events of correlated individuals is required to ensure a correct inference on the effects of treatments or covariates on the survival times. In the past decades, extension of proportional hazards model has been widely considered for modelling multivariate survival data by incorporating a random effect which acts multiplicatively on the hazard function. In this article, we consider the proportional odds model, which is an alternative to the proportional hazards model at which the hazard ratio between individuals converges to unity eventually. This is a reasonable property particularly when the treatment effect fades out gradually and the homogeneity of the population increases over time. The objective of this paper is to assess the influence of the random effect on the within‐subject correlation and the population heterogeneity. We are particularly interested in the properties of the proportional odds model with univariate random effect and correlated random effect. The correlations between survival times are derived explicitly for both choices of mixing distributions and are shown to be independent of the covariates. The time path of the odds function among the survivors are also examined to study the effect of the choice of mixing distribution. Modelling multivariate survival data using a univariate mixing distribution may be inadequate as the random effect not only characterises the dependence of the survival times, but also the conditional heterogeneity among the survivors. A robust estimate for the correlation of the logarithm of the survival times within a cluster is obtained disregarding the choice of the mixing distributions. The sensitivity of the estimate of the regression parameter under a misspecification of the mixing distribution is studied through simulation. (© 2004 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

Effect of covariate omission in Weibull accelerated failure time model: A caution

The accelerated failure time model is presented as an alternative to the proportional hazard model in the analysis of survival data. We investigate the effect of covariates omission in the case of applying a Weibull accelerated failure time model. In an uncensored setting, the asymptotic bias of the treatment effect is theoretically zero when important covariates are omitted; however, the asymptotic variance estimator of the treatment effect could be biased and then the size of the Wald test for the treatment effect is likely to exceed the nominal level. In some cases, the test size could be more than twice the nominal level. In a simulation study, in both censored and uncensored settings, Type I error for the test of the treatment effect was likely inflated when the prognostic covariates are omitted. This work remarks the careless use of the accelerated failure time model. We recommend the use of the robust sandwich variance estimator in order to avoid the inflation of the Type I error in the accelerated failure time model, although the robust variance is not commonly used in the survival data analyses.     

A Bayesian joint model of recurrent events and a terminal event

Recurrent events could be stopped by a terminal event, which commonly occurs in biomedical and clinical studies. In this situation, dependent censoring is encountered because of potential dependence between these two event processes, leading to invalid inference if analyzing recurrent events alone. The joint frailty model is one of the widely used approaches to jointly model these two processes by sharing the same frailty term. One important assumption is that recurrent and terminal event processes are conditionally independent given the subject‐level frailty; however, this could be violated when the dependency may also depend on time‐varying covariates across recurrences. Furthermore, marginal correlation between two event processes based on traditional frailty modeling has no closed form solution for estimation with vague interpretation. In order to fill these gaps, we propose a novel joint frailty‐copula approach to model recurrent events and a terminal event with relaxed assumptions. Metropolis–Hastings within the Gibbs Sampler algorithm is used for parameter estimation. Extensive simulation studies are conducted to evaluate the efficiency, robustness, and predictive performance of our proposal. The simulation results show that compared with the joint frailty model, the bias and mean squared error of the proposal is smaller when the conditional independence assumption is violated. Finally, we apply our method into a real example extracted from the MarketScan database to study the association between recurrent strokes and mortality.     

Systematic review of long term anticoagulation or antiplatelet treatment in patients with non-rheumatic atrial fibrillation

ObjectiveTo examine the benefits and risks of long term anticoagulation (warfarin) compared with antiplatelet treatment (aspirin/indoprofen) in patients with non-rheumatic atrial fibrillation.MethodsMeta-analysis of randomised controlled trials from Cochrane library, Medline, Embase, Cinhal, and Sigle from 1966 to December 1999. Odds ratios (95% confidence intervals) calculated to estimate treatment effects.ResultsNo trials were found from before 1989. There were five randomised controlled trials published between 1989-99. There were no significant differences in mortality between the two treatment options (fixed effects model: odd ratio 0.74 (95% confidence interval 0.39 to 1.40) for stroke deaths; 0.86 (0.63 to 1.17) for vascular deaths). There was a borderline significant difference in non-fatal stroke in favour of anticoagulation (0.68 (0.46 to 0.99)); and 0.75 (0.50 to 1.13) after exclusion of one trial with weak methodological design. A random effects model showed no significant difference in combined fatal and non-fatal events (odds ratio 0.79 (0.61 to 1.02)). There were more major bleeding events among patients on anticoagulation than on antiplatelet treatment (odds ratio 1.45 (0.93 to 2.27)). One trial was stopped prematurely after a significant difference in favour of anticoagulation was observed. The only trial to show a significant difference in effect (favouring anticoagulation) was methodologically weaker in design than the others.ConclusionsThe heterogeneity between the trials and the limited data result in considerable uncertainty about the value of long term anticoagulation compared with antiplatelet treatment. The risks of bleeding and the higher cost of anticoagulation make it an even less convincing treatment option.