Bayesian nonparametric analysis of restricted mean survival time

The restricted mean survival time (RMST) evaluates the expectation of survival time truncated by a prespecified time point, because the mean survival time in the presence of censoring is typically not estimable. The frequentist inference procedure for RMST has been widely advocated for comparison of two survival curves, while research from the Bayesian perspective is rather limited. For the RMST of both right- and interval-censored data, we propose Bayesian nonparametric estimation and inference procedures. By assigning a mixture of Dirichlet processes (MDP) prior to the distribution function, we can estimate the posterior distribution of RMST. We also explore another Bayesian nonparametric approach using the Dirichlet process mixture model and make comparisons with the frequentist nonparametric method. Simulation studies demonstrate that the Bayesian nonparametric RMST under diffuse MDP priors leads to robust estimation and under informative priors it can incorporate prior knowledge into the nonparametric estimator. Analysis of real trial examples demonstrates the flexibility and interpretability of the Bayesian nonparametric RMST for both right- and interval-censored data.     

Sequential methods for comparing years of life saved in the two-sample censored data problem

This research develops nonparametric strategies for sequentially monitoring clinical trial data where detecting years of life saved is of interest. The recommended test statistic looks at integrated differences in survival estimates during the time frame of interest. In many practical situations, the test statistic presented has an independent increments covariance structure. Hence, with little additional work, we may apply these testing procedures using available methodology. In the case where an independent increments covariance structure is present, we suggest how clinical trial data might be monitored using these statistics in an information-based design. The resulting study design maintains the desired stochastic operating characteristics regardless of the shapes of the survival curves being compared. This offers an advantage over the popular log-rank-based design strategy since more restrictive assumptions relating to the behavior of the hazards are required to guarantee the planned power of the test. Recommendations for how to sequentially monitor clinical trial progress in the nonindependent increments case are also provided along with an example.     

A Bayesian semiparametric joint hierarchical model for longitudinal and survival data

This article proposes a new semiparametric Bayesian hierarchical model for the joint modeling of longitudinal and survival data. We relax the distributional assumptions for the longitudinal model using Dirichlet process priors on the parameters defining the longitudinal model. The resulting posterior distribution of the longitudinal parameters is free of parametric constraints, resulting in more robust estimates. This type of approach is becoming increasingly essential in many applications, such as HIV and cancer vaccine trials, where patients' responses are highly diverse and may not be easily modeled with known distributions. An example will be presented from a clinical trial of a cancer vaccine where the survival outcome is time to recurrence of a tumor. Immunologic measures believed to be predictive of tumor recurrence were taken repeatedly during follow-up. We will present an analysis of this data using our new semiparametric Bayesian hierarchical joint modeling methodology to determine the association of these longitudinal immunologic measures with time to tumor recurrence.     

Nonparametric bayes testing of changes in a response distribution with an ordinal predictor

Summary .   In certain biomedical studies, one may anticipate changes in the shape of a response distribution across the levels of an ordinal predictor. For instance, in toxicology studies, skewness and modality might change as dose increases. To address this issue, we propose a Bayesian nonparametric method for testing for distribution changes across an ordinal predictor. Using a dynamic mixture of Dirichlet processes, we allow the response distribution to change flexibly at each level of the predictor. In addition, by assigning mixture priors to the hyperparameters, we can obtain posterior probabilities of no effect of the predictor and identify the lowest dose level for which there is an appreciable change in distribution. The method also provides a natural framework for performing tests across multiple outcomes. We apply our method to data from a genotoxicity experiment.     

Sample size and power for the weighted log-rank test and Kaplan-Meier based tests with allowance for nonproportional hazards

Asymptotic distributions under alternative hypotheses and their corresponding sample size and power equations are derived for nonparametric test statistics commonly used to compare two survival curves. Test statistics include the weighted log-rank test and the Wald test for difference in (or ratio of) Kaplan-Meier survival probability, percentile survival, and restricted mean survival time. Accrual, survival, and loss to follow-up are allowed to follow any arbitrary continuous distribution. We show that Schoenfeld's equation—often used by practitioners to calculate the required number of events for the unweighted log-rank test—can be inaccurate even when the proportional hazards (PH) assumption holds. In fact, it can mislead one to believe that 1:1 is the optimal randomization ratio (RR), when actually power can be gained by assigning more patients to the active arm. Meaningful improvements to Schoenfeld's equation are made. The present theory should be useful in designing clinical trials, particularly in immuno-oncology where nonproportional hazards are frequently encountered. We illustrate the application of our theory with an example exploring optimal RR under PH and a second example examining the impact of delayed treatment effect. A companion R package npsurvSS is available for download on CRAN.     

Bayesian Hierarchical Functional Data Analysis Via Contaminated Informative Priors

Summary .  A variety of flexible approaches have been proposed for functional data analysis, allowing both the mean curve and the distribution about the mean to be unknown. Such methods are most useful when there is limited prior information. Motivated by applications to modeling of temperature curves in the menstrual cycle, this article proposes a flexible approach for incorporating prior information in semiparametric Bayesian analyses of hierarchical functional data. The proposed approach is based on specifying the distribution of functions as a mixture of a parametric hierarchical model and a nonparametric contamination. The parametric component is chosen based on prior knowledge, while the contamination is characterized as a functional Dirichlet process. In the motivating application, the contamination component allows unanticipated curve shapes in unhealthy menstrual cycles. Methods are developed for posterior computation, and the approach is applied to data from a European fecundability study.     

A non-parametric mixture model for genome-enabled prediction of genetic value for a quantitative trait

A Bayesian nonparametric form of regression based on Dirichlet process priors is adapted to the analysis of quantitative traits possibly affected by cryptic forms of gene action, and to the context of SNP-assisted genomic selection, where the main objective is to predict a genomic signal on phenotype. The procedure clusters unknown genotypes into groups with distinct genetic values, but in a setting in which the number of clusters is unknown a priori, so that standard methods for finite mixture analysis do not work. The central assumption is that genetic effects follow an unknown distribution with some “baseline” family, which is a normal process in the cases considered here. A Bayesian analysis based on the Gibbs sampler produces estimates of the number of clusters, posterior means of genetic effects, a measure of credibility in the baseline distribution, as well as estimates of parameters of the latter. The procedure is illustrated with a simulation representing two populations. In the first one, there are 3 unknown QTL, with additive, dominance and epistatic effects; in the second, there are 10 QTL with additive, dominance and additive × additive epistatic effects. In the two populations, baseline parameters are inferred correctly. The Dirichlet process model infers the number of unique genetic values correctly in the first population, but it produces an understatement in the second one; here, the true number of clusters is over 900, and the model gives a posterior mean estimate of about 140, probably because more replication of genotypes is needed for correct inference. The impact on inferences of the prior distribution of a key parameter (M), and of the extent of replication, was examined via an analysis of mean body weight in 192 paternal half-sib families of broiler chickens, where each sire was genotyped for nearly 7,000 SNPs. In this small sample, it was found that inference about the number of clusters was affected by the prior distribution of M. For a set of combinations of parameters of a given prior distribution, the effects of the prior dissipated when the number of replicate samples per genotype was increased. Thus, the Dirichlet process model seems to be useful for gauging the number of QTLs affecting the trait: if the number of clusters inferred is small, probably just a few QTLs code for the trait. If the number of clusters inferred is large, this may imply that standard parametric models based on the baseline distribution may suffice. However, priors may be influential, especially if sample size is not large and if only a few genotypic configurations have replicate phenotypes in the sample.     

Bayesian estimation in animal breeding using the Dirichlet process prior for correlated random effects

In the case of the mixed linear model the random effects are usually assumed to be normally distributed in both the Bayesian and classical frameworks. In this paper, the Dirichlet process prior was used to provide nonparametric Bayesian estimates for correlated random effects. This goal was achieved by providing a Gibbs sampler algorithm that allows these correlated random effects to have a nonparametric prior distribution. A sampling based method is illustrated. This method which is employed by transforming the genetic covariance matrix to an identity matrix so that the random effects are uncorrelated, is an extension of the theory and the results of previous researchers. Also by using Gibbs sampling and data augmentation a simulation procedure was derived for estimating the precision parameter M associated with the Dirichlet process prior. All needed conditional posterior distributions are given. To illustrate the application, data from the Elsenburg Dormer sheep stud were analysed. A total of 3325 weaning weight records from the progeny of 101 sires were used.     

Bayesian monitoring of event rates with censored data

A Bayesian approach to monitoring event rates with censored data is proposed. A Dirichlet prior for discrete time event probabilities is blended with discrete survival times to provide a posterior distribution that is a mixture of Dirichlets. Approximation of the posterior distribution via data augmentation is discussed. Practical issues involved in implementing the procedure are discussed and illustrated with a simulation of the single arm Cord Blood Transplantation Study where 6-month survival is monitored.     

Nonparametric estimation of the joint distribution of a survival time subject to interval censoring and a continuous mark variable

This article considers three nonparametric estimators of the joint distribution function for a survival time and a continuous mark variable when the survival time is interval censored and the mark variable may be missing for interval-censored observations. Finite and large sample properties are described for the nonparametric maximum likelihood estimator (NPMLE) as well as estimators based on midpoint imputation (MIDMLE) and coarsening the mark variable (CMLE). The estimators are compared using data from a simulation study and a recent phase III HIV vaccine efficacy trial where the survival time is the time from enrollment to infection and the mark variable is the genetic distance from the infecting HIV sequence to the HIV sequence in the vaccine. Theoretical and empirical evidence are presented indicating the NPMLE and MIDMLE are inconsistent. Conversely, the CMLE is shown to be consistent in general and thus is preferred.     

Bayesian Nonparametric Nonproportional Hazards Survival Modeling

Summary .  We develop a dependent Dirichlet process model for survival analysis data. A major feature of the proposed approach is that there is no necessity for resulting survival curve estimates to satisfy the ubiquitous proportional hazards assumption. An illustration based on a cancer clinical trial is given, where survival probabilities for times early in the study are estimated to be lower for those on a high-dose treatment regimen than for those on the low dose treatment, while the reverse is true for later times, possibly due to the toxic effect of the high dose for those who are not as healthy at the beginning of the study.     

Bayesian nonparametric centered random effects models with variable selection

In a linear mixed effects model, it is common practice to assume that the random effects follow a parametric distribution such as a normal distribution with mean zero. However, in the case of variable selection, substantial violation of the normality assumption can potentially impact the subset selection and result in poor interpretation and even incorrect results. In nonparametric random effects models, the random effects generally have a nonzero mean, which causes an identifiability problem for the fixed effects that are paired with the random effects. In this article, we focus on a Bayesian method for variable selection. We characterize the subject‐specific random effects nonparametrically with a Dirichlet process and resolve the bias simultaneously. In particular, we propose flexible modeling of the conditional distribution of the random effects with changes across the predictor space. The approach is implemented using a stochastic search Gibbs sampler to identify subsets of fixed effects and random effects to be included in the model. Simulations are provided to evaluate and compare the performance of our approach to the existing ones. We then apply the new approach to a real data example, cross‐country and interlaboratory rodent uterotrophic bioassay.     

Modeling disease incidence data with spatial and spatio temporal dirichlet process mixtures

Disease incidence or mortality data are typically available as rates or counts for specified regions, collected over time. We propose Bayesian nonparametric spatial modeling approaches to analyze such data. We develop a hierarchical specification using spatial random effects modeled with a Dirichlet process prior. The Dirichlet process is centered around a multivariate normal distribution. This latter distribution arises from a log-Gaussian process model that provides a latent incidence rate surface, followed by block averaging to the areal units determined by the regions in the study. With regard to the resulting posterior predictive inference, the modeling approach is shown to be equivalent to an approach based on block averaging of a spatial Dirichlet process to obtain a prior probability model for the finite dimensional distribution of the spatial random effects. We introduce a dynamic formulation for the spatial random effects to extend the model to spatio-temporal settings. Posterior inference is implemented through Gibbs sampling. We illustrate the methodology with simulated data as well as with a data set on lung cancer incidences for all 88 counties in the state of Ohio over an observation period of 21 years.     

Time-dependent ROC curves for censored survival data and a diagnostic marker

ROC curves are a popular method for displaying sensitivity and specificity of a continuous diagnostic marker, X, for a binary disease variable, D. However, many disease outcomes are time dependent, D(t), and ROC curves that vary as a function of time may be more appropriate. A common example of a time-dependent variable is vital status, where D(t) = 1 if a patient has died prior to time t and zero otherwise. We propose summarizing the discrimination potential of a marker X, measured at baseline (t = 0), by calculating ROC curves for cumulative disease or death incidence by time t, which we denote as ROC(t). A typical complexity with survival data is that observations may be censored. Two ROC curve estimators are proposed that can accommodate censored data. A simple estimator is based on using the Kaplan-Meier estimator for each possible subset X > c. However, this estimator does not guarantee the necessary condition that sensitivity and specificity are monotone in X. An alternative estimator that does guarantee monotonicity is based on a nearest neighbor estimator for the bivariate distribution function of (X, T), where T represents survival time (Akritas, M. J., 1994, Annals of Statistics 22, 1299-1327). We present an example where ROC(t) is used to compare a standard and a modified flow cytometry measurement for predicting survival after detection of breast cancer and an example where the ROC(t) curve displays the impact of modifying eligibility criteria for sample size and power in HIV prevention trials.     

Nonparametric Analysis of Ordered Categorical Data in Designs with Longitudinal Observations and Small Sample Sizes

For designs with longitudinal observations of ordered categorical data, a nonparametric model is considered where treatment effects and interactions are defined by means of the marginal distributions. These treatment effects are estimated consistently by ranking methods. The hypotheses in this nonparametric setup are formulated by means of the distribution functions. The asymptotic distribution of the estimators for the nonparametric effects are given under the hypotheses. For small samples, a rather accurate approximation is suggested. A clinical trial with ordered categorical data is used to motivate the ideas and to explain the procedures which are extensions of the Wilcoxon‐Mann‐Whitney test to factorial designs with longitudinal observations. The application of the procedures requires only some trivial regularity assumptions.     

The Nonparametric Behrens‐Fisher Problem: Asymptotic Theory and a Small‐Sample Approximation

A generalization of the Behrens‐Fisher problem for two samples is examined in a nonparametric model. It is not assumed that the underlying distribution functions are continuous so that data with arbitrary ties can be handled. A rank test is considered where the asymptotic variance is estimated consistently by using the ranks over all observations as well as the ranks within each sample. The consistency of the estimator is derived in the appendix. For small samples (n₁, n₂ ≥ 10), a simple approximation by a central t‐distribution is suggested where the degrees of freedom are taken from the Satterthwaite‐Smith‐Welch approximation in the parametric Behrens‐Fisher problem. It is demonstrated by means of a simulation study that the Wilcoxon‐Mann‐Whitney‐test may be conservative or liberal depending on the ratio of the sample sizes and the variances of the underlying distribution functions. For the suggested approximation, however, it turns out that the nominal level is maintained rather accurately. The suggested nonparametric procedure is applied to a data set from a clinical trial. Moreover, a confidence interval for the nonparametric treatment effect is given.     

Bayesian case influence diagnostics for survival models

We propose Bayesian case influence diagnostics for complex survival models. We develop case deletion influence diagnostics for both the joint and marginal posterior distributions based on the Kullback-Leibler divergence (K-L divergence). We present a simplified expression for computing the K-L divergence between the posterior with the full data and the posterior based on single case deletion, as well as investigate its relationships to the conditional predictive ordinate. All the computations for the proposed diagnostic measures can be easily done using Markov chain Monte Carlo samples from the full data posterior distribution. We consider the Cox model with a gamma process prior on the cumulative baseline hazard. We also present a theoretical relationship between our case-deletion diagnostics and diagnostics based on Cox's partial likelihood. A simulated data example and two real data examples are given to demonstrate the methodology.     

Bayesian functional integral method for inferring continuous data from discrete measurements

Inference of the insulin secretion rate (ISR) from C-peptide measurements as a quantification of pancreatic β-cell function is clinically important in diseases related to reduced insulin sensitivity and insulin action. ISR derived from C-peptide concentration is an example of nonparametric Bayesian model selection where a proposed ISR time-course is considered to be a "model". An inferred value of inaccessible continuous variables from discrete observable data is often problematic in biology and medicine, because it is a priori unclear how robust the inference is to the deletion of data points, and a closely related question, how much smoothness or continuity the data actually support. Predictions weighted by the posterior distribution can be cast as functional integrals as used in statistical field theory. Functional integrals are generally difficult to evaluate, especially for nonanalytic constraints such as positivity of the estimated parameters. We propose a computationally tractable method that uses the exact solution of an associated likelihood function as a prior probability distribution for a Markov-chain Monte Carlo evaluation of the posterior for the full model. As a concrete application of our method, we calculate the ISR from actual clinical C-peptide measurements in human subjects with varying degrees of insulin sensitivity. Our method demonstrates the feasibility of functional integral Bayesian model selection as a practical method for such data-driven inference, allowing the data to determine the smoothing timescale and the width of the prior probability distribution on the space of models. In particular, our model comparison method determines the discrete time-step for interpolation of the unobservable continuous variable that is supported by the data. Attempts to go to finer discrete time-steps lead to less likely models.     

Time-dependent mediators in survival analysis: Modeling direct and indirect effects with the additive hazards model

We discuss causal mediation analyses for survival data and propose a new approach based on the additive hazards model. The emphasis is on a dynamic point of view, that is, understanding how the direct and indirect effects develop over time. Hence, importantly, we allow for a time varying mediator. To define direct and indirect effects in such a longitudinal survival setting we take an interventional approach (Didelez, 2018) where treatment is separated into one aspect affecting the mediator and a different aspect affecting survival. In general, this leads to a version of the nonparametric g-formula (Robins, 1986). In the present paper, we demonstrate that combining the g-formula with the additive hazards model and a sequential linear model for the mediator process results in simple and interpretable expressions for direct and indirect effects in terms of relative survival as well as cumulative hazards. Our results generalize and formalize the method of dynamic path analysis (Fosen, Ferkingstad, Borgan, & Aalen, 2006; Strohmaier et al., 2015). An application to data from a clinical trial on blood pressure medication is given.     

A modification of Peto's nonparametric estimation of survival curves for interval-censored data

Peto (1973, Applied Statistics, 22, 86-91) gave a nonparametric generalized maximum-likelihood estimate of the survival function for interval-censored data. His method has a tendency to concentrate probability masses at the endpoints of the intervals, even for the ordinary grouped data, instead of spreading them through the intervals, as one might expect them to be in the underlying distribution. We describe a modification that overcomes this. The new estimate reduces to the standard binomial estimate when applied to grouped data. It also reduces to the Kaplan-Meier estimate when applied to survival data that consist of only exact or right-censored observations. Both estimates are maximum-likelihood estimates but are based on different interpretations of the endpoints of the intervals.