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.     

Analysis of a molecular genetic neuro-oncology study with partially biased selection

Oligodendrogliomas are a common variant of malignant brain tumors, and are unique for their relative sensitivity to chemotherapy and better prognosis. For these reasons, the identification of an objective oligodendroglial marker has been a long sought-after goal in the field of neuro-oncology. To this end, 75 patients who received chemotherapy at the London Regional Cancer Centre between 1984 and 1999 were studied (Ino et al., Clinical Cancer Research, 7, 839-845, 2001). Of these 75 patients, 50 were initially treated with chemotherapy (the current practice) and comprise a population-based sample. The remaining 25 patients were initially treated with radiation and were included in the study only because their tumor recurred, at which time they received chemotherapy. Because this group of 25 patients included neither those radiation patients whose tumors never recurred nor those radiation patients whose tumors recurred but were not treated with chemotherapy, issues of selection bias were of concern. For this reason, the initial analysis of these data included only the 50 population-based patients. This was unsatisfying given the rarity of this disease and of genetic information on this disease and led us to question whether we could undertake an analysis that includes all of the patients. Here we examine approaches for utilizing the entire study population, as well as the assumptions required for doing so. We illustrate that there are both costs and benefits to using the 25 selected patients.     

Analysis of a partially observed binary covariate process and a censored failure time in the presence of truncation and competing risks

We develop methods for assessing the association between a binary time-dependent covariate process and a failure time endpoint when the former is observed only at a single time point and the latter is right censored, and when the observations are subject to truncation and competing causes of failure. Using a proportional hazards model for the effect of the covariate process on the failure time of interest, we develop an approach utilizing EM algorithm and profile likelihood for estimating the relative risk parameter and cause-specific hazards for failure. The methods are extended to account for other covariates that can influence the time-dependent covariate process and cause-specific risks of failure. We illustrate the methods with data from a recent study on the association between loss of hepatitis B e antigen and the development of hepatocellular carcinoma in a population of chronic carriers of hepatitis B.     

Parametric survival models for interval-censored data with time-dependent covariates

We present a parametric family of regression models for interval-censored event-time (survival) data that accomodates both fixed (e.g. baseline) and time-dependent covariates. The model employs a three-parameter family of survival distributions that includes the Weibull, negative binomial, and log-logistic distributions as special cases, and can be applied to data with left, right, interval, or non-censored event times. Standard methods, such as Newton-Raphson, can be employed to estimate the model and the resulting estimates have an asymptotically normal distribution about the true values with a covariance matrix that is consistently estimated by the information function. The deviance function is described to assess model fit and a robust sandwich estimate of the covariance may also be employed to provide asymptotically robust inferences when the model assumptions do not apply. Spline functions may also be employed to allow for non-linear covariates. The model is applied to data from a long-term study of type 1 diabetes to describe the effects of longitudinal measures of glycemia (HbA1c) over time (the time-dependent covariate) on the risk of progression of diabetic retinopathy (eye disease), an interval-censored event-time outcome.     

Regression Analysis with a Misclassified Covariate from a Current Status Observation Scheme

Summary Naive use of misclassified covariates leads to inconsistent estimators of covariate effects in regression models. A variety of methods have been proposed to address this problem including likelihood, pseudo‐likelihood, estimating equation methods, and Bayesian methods, with all of these methods typically requiring either internal or external validation samples or replication studies. We consider a problem arising from a series of orthopedic studies in which interest lies in examining the effect of a short‐term serological response and other covariates on the risk of developing a longer term thrombotic condition called deep vein thrombosis. The serological response is an indicator of whether the patient developed antibodies following exposure to an antithrombotic drug, but the seroconversion status of patients is only available at the time of a blood sample taken upon the discharge from hospital. The seroconversion time is therefore subject to a current status observation scheme, or Case I interval censoring, and subjects tested before seroconversion are misclassified as nonseroconverters. We develop a likelihood‐based approach for fitting regression models that accounts for misclassification of the seroconversion status due to early testing using parametric and nonparametric estimates of the seroconversion time distribution. The method is shown to reduce the bias resulting from naive analyses in simulation studies and an application to the data from the orthopedic studies provides further illustration.     

Borrowing strength with nonexchangeable priors over subpopulations

We introduce a nonparametric Bayesian model for a phase II clinical trial with patients presenting different subtypes of the disease under study. The objective is to estimate the success probability of an experimental therapy for each subtype. We consider the case when small sample sizes require extensive borrowing of information across subtypes, but the subtypes are not a priori exchangeable. The lack of a priori exchangeability hinders the straightforward use of traditional hierarchical models to implement borrowing of strength across disease subtypes. We introduce instead a random partition model for the set of disease subtypes. This is a variation of the product partition model that allows us to model a nonexchangeable prior structure. Like a hierarchical model, the proposed clustering approach considers all observations, across all disease subtypes, to estimate individual success probabilities. But in contrast to standard hierarchical models, the model considers disease subtypes a priori nonexchangeable. This implies that when assessing the success probability for a particular type our model borrows more information from the outcome of the patients sharing the same prognosis than from the others. Our data arise from a phase II clinical trial of patients with sarcoma, a rare type of cancer affecting connective or supportive tissues and soft tissue (e.g., cartilage and fat). Each patient presents one subtype of the disease and subtypes are grouped by good, intermediate, and poor prognosis. The prior model should respect the varying prognosis across disease subtypes. The practical motivation for the proposed approach is that the number of accrued patients within each disease subtype is small. Thus it is not possible to carry out a clinical study of possible new therapies for rare conditions, because it would be impossible to plan for sufficiently large sample size to achieve the desired power. We carry out a simulation study to compare the proposed model with a model that assumes similar success probabilities for all subtypes with the same prognosis, i.e., a fixed partition of subtypes by prognosis. When the assumption is satisfied the two models perform comparably. But the proposed model outperforms the competing model when the assumption is incorrect.     

Bayesian estimators for conditional hazard functions

This article introduces a new Bayesian approach to the analysis of right-censored survival data. The hazard rate of interest is modeled as a product of conditionally independent stochastic processes corresponding to (1) a baseline hazard function and (2) a regression function representing the temporal influence of the covariates. These processes jump at times that form a time-homogeneous Poisson process and have a pairwise dependency structure for adjacent values. The two processes are assumed to be conditionally independent given their jump times. Features of the posterior distribution, such as the mean covariate effects and survival probabilities (conditional on the covariate), are evaluated using the Metropolis-Hastings-Green algorithm. We illustrate our methodology by an application to nasopharynx cancer survival data.     

Estimated Pseudopartial‐Likelihood Method for Correlated Failure Time Data with Auxiliary Covariates

Summary As biological studies become more expensive to conduct, statistical methods that take advantage of existing auxiliary information about an expensive exposure variable are desirable in practice. Such methods should improve the study efficiency and increase the statistical power for a given number of assays. In this article, we consider an inference procedure for multivariate failure time with auxiliary covariate information. We propose an estimated pseudopartial likelihood estimator under the marginal hazard model framework and develop the asymptotic properties for the proposed estimator. We conduct simulation studies to evaluate the performance of the proposed method in practical situations and demonstrate the proposed method with a data set from the studies of left ventricular dysfunction ( SOLVD Investigators, 1991 , New England Journal of Medicine 325 , 293–302).     

Identifiability of causal effects for binary variables with baseline data missing due to death

We discuss identifiability and estimation of causal effects of a treatment in subgroups defined by a covariate that is sometimes missing due to death, which is different from a problem with outcomes censored by death. Frangakis et al. (2007, Biometrics 63, 641-662) proposed an approach for estimating the causal effects under a strong monotonicity (SM) assumption. In this article, we focus on identifiability of the joint distribution of the covariate, treatment and potential outcomes, show sufficient conditions for identifiability, and relax the SM assumption to monotonicity (M) and no-interaction (NI) assumptions. We derive expectation-maximization algorithms for finding the maximum likelihood estimates of parameters of the joint distribution under different assumptions. Further we remove the M and NI assumptions, and prove that signs of the causal effects of a treatment in the subgroups are identifiable, which means that their bounds do not cover zero. We perform simulations and a sensitivity analysis to evaluate our approaches. Finally, we apply the approaches to the National Study on the Costs and Outcomes of Trauma Centers data, which are also analyzed by Frangakis et al. (2007) and Xie and Murphy (2007, Biometrics 63, 655-658).     

ROC-guided survival trees and ensembles

Tree-based methods are popular nonparametric tools in studying time-to-event outcomes. In this article, we introduce a novel framework for survival trees and ensembles, where the trees partition the dynamic survivor population and can handle time-dependent covariates. Using the idea of randomized tests, we develop generalized time-dependent receiver operating characteristic (ROC) curves for evaluating the performance of survival trees. The tree-building algorithm is guided by decision-theoretic criteria based on ROC, targeting specifically for prediction accuracy. To address the instability issue of a single tree, we propose a novel ensemble procedure based on averaging martingale estimating equations, which is different from existing methods that average the predicted survival or cumulative hazard functions from individual trees. Extensive simulation studies are conducted to examine the performance of the proposed methods. We apply the methods to a study on AIDS for illustration.     

Cox regression model with randomly censored covariates

This paper deals with a Cox proportional hazards regression model, where some covariates of interest are randomly right‐censored. While methods for censored outcomes have become ubiquitous in the literature, methods for censored covariates have thus far received little attention and, for the most part, dealt with the issue of limit‐of‐detection. For randomly censored covariates, an often‐used method is the inefficient complete‐case analysis (CCA) which consists in deleting censored observations in the data analysis. When censoring is not completely independent, the CCA leads to biased and spurious results. Methods for missing covariate data, including type I and type II covariate censoring as well as limit‐of‐detection do not readily apply due to the fundamentally different nature of randomly censored covariates. We develop a novel method for censored covariates using a conditional mean imputation based on either Kaplan–Meier estimates or a Cox proportional hazards model to estimate the effects of these covariates on a time‐to‐event outcome. We evaluate the performance of the proposed method through simulation studies and show that it provides good bias reduction and statistical efficiency. Finally, we illustrate the method using data from the Framingham Heart Study to assess the relationship between offspring and parental age of onset of cardiovascular events.     

On the relevance of prognostic information for clinical trials: A theoretical quantification

The question of how individual patient data from cohort studies or historical clinical trials can be leveraged for designing more powerful, or smaller yet equally powerful, clinical trials becomes increasingly important in the era of digitalization. Today, the traditional statistical analyses approaches may seem questionable to practitioners in light of ubiquitous historical prognostic information. Several methodological developments aim at incorporating historical information in the design and analysis of future clinical trials, most importantly Bayesian information borrowing, propensity score methods, stratification, and covariate adjustment. Adjusting the analysis with respect to a prognostic score, which was obtained from some model applied to historical data, received renewed interest from a machine learning perspective, and we study the potential of this approach for randomized clinical trials. In an idealized situation of a normal outcome in a two-arm trial with 1:1 allocation, we derive a simple sample size reduction formula as a function of two criteria characterizing the prognostic score: (1) the coefficient of determination R² on historical data and (2) the correlation ρ between the estimated and the true unknown prognostic scores. While maintaining the same power, the original total sample size n planned for the unadjusted analysis reduces to $(1-R^2{\rho }^2)\times n$ in an adjusted analysis. Robustness in less ideal situations was assessed empirically. We conclude that there is potential for substantially more powerful or smaller trials, but only when prognostic scores can be accurately estimated.     

Dose–response relationship from longitudinal data with response‐dependent dose modification using likelihood methods

In some clinical trials or clinical practice, the therapeutic agent is administered repeatedly, and doses are adjusted in each patient based on repeatedly measured continuous responses, to maintain the response levels in a target range. Because a lower dose tends to be selected for patients with a better outcome, simple summarizations may wrongly show a better outcome for the lower dose, producing an incorrect dose–response relationship. In this study, we consider the dose–response relationship under these situations. We show that maximum‐likelihood estimates are consistent without modeling the dose‐modification mechanisms when the selection of the dose as a time‐dependent covariate is based only on observed, but not on unobserved, responses, and measurements are generated based on administered doses. We confirmed this property by performing simulation studies under several dose‐modification mechanisms. We examined an autoregressive linear mixed effects model. The model represents profiles approaching each patient's asymptote when identical doses are repeatedly administered. The model takes into account the previous dose history and provides a dose–response relationship of the asymptote as a summary measure. We also examined a linear mixed effects model assuming all responses are measured at steady state. In the simulation studies, the estimates of both the models were unbiased under the dose modification based on observed responses, but biased under the dose modification based on unobserved responses. In conclusion, the maximum‐likelihood estimates of the dose–response relationship are consistent under the dose modification based only on observed responses.     

A matching procedure for sequential experiments that iteratively learns which covariates improve power

We propose a dynamic allocation procedure that increases power and efficiency when measuring an average treatment effect in sequential randomized trials exploiting some subjects' previous assessed responses. Subjects arrive sequentially and are either randomized or paired to a previously randomized subject and administered the alternate treatment. The pairing is made via a dynamic matching criterion that iteratively learns which specific covariates are important to the response. We develop estimators for the average treatment effect as well as an exact test. We illustrate our method's increase in efficiency and power over other allocation procedures in both simulated scenarios and a clinical trial dataset. An R package “SeqExpMatch ” for use by practitioners is available on CRAN .     

Estimating treatment effects with partially observed covariates using outcome regression with missing indicators

Missing data is a common issue in research using observational studies to investigate the effect of treatments on health outcomes. When missingness occurs only in the covariates, a simple approach is to use missing indicators to handle the partially observed covariates. The missing indicator approach has been criticized for giving biased results in outcome regression. However, recent papers have suggested that the missing indicator approach can provide unbiased results in propensity score analysis under certain assumptions. We consider assumptions under which the missing indicator approach can provide valid inferences, namely, (1) no unmeasured confounding within missingness patterns; either (2a) covariate values of patients with missing data were conditionally independent of treatment or (2b) these values were conditionally independent of outcome; and (3) the outcome model is correctly specified: specifically, the true outcome model does not include interactions between missing indicators and fully observed covariates. We prove that, under the assumptions above, the missing indicator approach with outcome regression can provide unbiased estimates of the average treatment effect. We use a simulation study to investigate the extent of bias in estimates of the treatment effect when the assumptions are violated and we illustrate our findings using data from electronic health records. In conclusion, the missing indicator approach can provide valid inferences for outcome regression, but the plausibility of its assumptions must first be considered carefully.     

Inference of haplotype effects in case-control studies using unphased genotype and environmental data

A retrospective likelihood-based approach was proposed to test and estimate the effect of haplotype on disease risk using unphased genotype data with adjustment for environmental covariates. The proposed method was also extended to handle the data in which the haplotype and environmental covariates are not independent. Likelihood ratio tests were constructed to test the effects of haplotype and gene-environment interaction. The model parameters such as haplotype effect size was estimated using an Expectation Conditional-Maximization (ECM) algorithm developed by Meng and Rubin (1993). Model-based variance estimates were derived using the observed information matrix. Simulation studies were conducted for three different genetic effect models, including dominant effect, recessive effect, and additive effect. The results showed that the proposed method generated unbiased parameter estimates, proper type I error, and true beta coverage probabilities. The model performed well with small or large sample sizes, as well as short or long haplotypes.     

Valid Adjustment of Randomized Comparisons for Binary Covariates

In randomized trials, the treatment influences not only endpoints but also other variables measured after randomization which, when used as covariates to adjust for the observed imbalance, become pseudo‐covariates. There is a logical circularity in adjusting for a pseudo‐covariate because the variability in the endpoint that is attributed not to the treatment but rather to the pseudo‐covariate may actually represent an effect of the treatment modulated by the pseudo‐covariate. This potential bias is well known, but we offer new insight into how it can lead to reversals in the direction of the apparent treatment effect by way of stage migration. We then discuss a related problem that is not generally appreciated, specifically how the absence of allocation concealment can lead to this reversal of the direction of the apparent treatment effect even when adjustment is for a true covariate measured prior to randomization. (© 2004 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

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.     

Testing for bias in weighted estimating equations

It is very common in regression analysis to encounter incompletely observed covariate information. A recent approach to analyse such data is weighted estimating equations (Robins, J. M., Rotnitzky, A. and Zhao, L. P. (1994), JASA, 89, 846-866, and Zhao, L. P., Lipsitz, S. R. and Lew, D. (1996), Biometrics, 52, 1165-1182). With weighted estimating equations, the contribution to the estimating equation from a complete observation is weighted by the inverse of the probability of being observed. We propose a test statistic to assess if the weighted estimating equations produce biased estimates. Our test statistic is similar to the test statistic proposed by DuMouchel and Duncan (1983) for weighted least squares estimates for sample survey data. The method is illustrated using data from a randomized clinical trial on chemotherapy for multiple myeloma.     

Estimating the area under the ROC curve when transporting a prediction model to a target population

We propose methods for estimating the area under the receiver operating characteristic (ROC) curve (AUC) of a prediction model in a target population that differs from the source population that provided the data used for original model development. If covariates that are associated with model performance, as measured by the AUC, have a different distribution in the source and target populations, then AUC estimators that only use data from the source population will not reflect model performance in the target population. Here, we provide identification results for the AUC in the target population when outcome and covariate data are available from the sample of the source population, but only covariate data are available from the sample of the target population. In this setting, we propose three estimators for the AUC in the target population and show that they are consistent and asymptotically normal. We evaluate the finite-sample performance of the estimators using simulations and use them to estimate the AUC in a nationally representative target population from the National Health and Nutrition Examination Survey for a lung cancer risk prediction model developed using source population data from the National Lung Screening Trial.