Distance correlation is a new class of multivariate dependence coefficients applicable to random vectors of arbitrary and not necessarily equal dimension. Distance covariance and distance correlation are analogous to product-moment covariance and correlation, but generalize and extend these classical bivariate measures of dependence. Distance correlation characterizes independence: it is zero if and only if the random vectors are independent. The notion of covariance with respect to a stochastic process is introduced, and it is shown that population distance covariance coincides with the covariance with respect to Brownian motion; thus, both can be called Brownian distance covariance. In the bivariate case, Brownian covariance is the natural extension of product-moment covariance, as we obtain Pearson product-moment covariance by replacing the Brownian motion in the definition with identity. The corresponding statistic has an elegantly simple computing formula. Advantages of applying Brownian covariance and correlation vs the classical Pearson covariance and correlation are discussed and illustrated.

We discuss briefly the very interesting concept of Brownian distance covariance developed by Székely and Rizzo [Ann. Appl. Statist. (2009), to appear] and describe two possible extensions. The first extension is for high dimensional data that can be coerced into a Hilbert space, including certain high throughput screening and functional data settings. The second extension involves very simple modifications that may yield increased power in some settings. We commend Székely and Rizzo for their very interesting work and recognize that this general idea has potential to have a large impact on the way in which statisticians evaluate dependency in data.

High-throughput biological assays such as microarrays let us ask very detailed questions about how diseases operate, and promise to let us personalize therapy. Data processing, however, is often not described well enough to allow for exact reproduction of the results, leading to exercises in “forensic bioinformatics” where aspects of raw data and reported results are used to infer what methods must have been employed. Unfortunately, poor documentation can shift from an inconvenience to an active danger when it obscures not just methods but errors. In this report we examine several related papers purporting to use microarray-based signatures of drug sensitivity derived from cell lines to predict patient response. Patients in clinical trials are currently being allocated to treatment arms on the basis of these results. However, we show in five case studies that the results incorporate several simple errors that may be putting patients at risk. One theme that emerges is that the most common errors are simple (e.g., row or column offsets); conversely, it is our experience that the most simple errors are common. We then discuss steps we are taking to avoid such errors in our own investigations.

A trend in all scientific disciplines, based on advances in technology, is the increasing availability of high dimensional data in which are buried important information. A current urgent challenge to statisticians is to develop effective methods of finding the useful information from the vast amounts of messy and noisy data available, most of which are noninformative. This paper presents a general computer intensive approach, based on a method pioneered by Lo and Zheng for detecting which, of many potential explanatory variables, have an influence on a dependent variable Y. This approach is suited to detect influential variables, where causal effects depend on the confluence of values of several variables. It has the advantage of avoiding a difficult direct analysis, involving possibly thousands of variables, by dealing with many randomly selected small subsets from which smaller subsets are selected, guided by a measure of influence I. The main objective is to discover the influential variables, rather than to measure their effects. Once they are detected, the problem of dealing with a much smaller group of influential variables should be vulnerable to appropriate analysis. In a sense, we are confining our attention to locating a few needles in a haystack.

Fiscal year-end balances of the Individual Indian Money System (a part of the Indian Trust) were constructed from data related to money collected in the system and disbursed by the system from 1887 to 2007. The data set of fiscal year accounting information had a high proportion of missing values, and much of the available data did not satisfy basic accounting relationships. Instead of just calculating a single estimate and arguing to the Court that the assumptions needed for the computation were reasonable, a distribution of calculated balances was developed using multiple imputation and time series models. These provided information to assess the uncertainty of the estimate due to missing and questionable data.

The T-wave of an electrocardiogram (ECG) represents the ventricular repolarization that is critical in restoration of the heart muscle to a pre-contractile state prior to the next beat. Alterations in the T-wave reflect various cardiac conditions; and links between abnormal (prolonged) ventricular repolarization and malignant arrhythmias have been documented. Cardiac safety testing prior to approval of any new drug currently relies on two points of the ECG waveform: onset of the Q-wave and termination of the T-wave; and only a few beats are measured. Using functional data analysis, a statistical approach extracts a common shape for each subject (reference curve) from a sequence of beats, and then models the deviation of each curve in the sequence from that reference curve as a four-dimensional vector. The representation can be used to distinguish differences between beats or to model shape changes in a subject’s T-wave over time. This model provides physically interpretable parameters characterizing T-wave shape, and is robust to the determination of the endpoint of the T-wave. Thus, this dimension reduction methodology offers the strong potential for definition of more robust and more informative biomarkers of cardiac abnormalities than the QT (or QT corrected) interval in current use.

Today’s call center managers face multiple operational decision-making tasks. One of the most common is determining the weekly staffing levels to ensure customer satisfaction and meeting their needs while minimizing service costs. An initial step for producing the weekly schedule is forecasting the future system loads which involves predicting both arrival counts and average service times.

We introduce an arrival count model which is based on a mixed Poisson process approach. The model is applied to data from an Israeli Telecom company call center. In our model, we also consider the effect of events such as billing on the arrival process and we demonstrate how to incorporate them as exogenous variables in the model.

After obtaining the forecasted system load, in large call centers, a manager can choose to apply the QED (Quality-Efficiency Driven) regime’s “square-root staffing” rule in order to balance the offered-load per server with the quality of service. Implementing this staffing rule requires that the forecasted values of the arrival counts and average service times maintain certain levels of precision. We develop different goodness of fit criteria that help determine our model’s practical performance under the QED regime. These show that during most hours of the day the model can reach desired precision levels.

The study of fingerprint individuality aims to determine to what extent a fingerprint uniquely identifies an individual. Recent court cases have highlighted the need for measures of fingerprint individuality when a person is identified based on fingerprint evidence. The main challenge in studies of fingerprint individuality is to adequately capture the variability of fingerprint features in a population. In this paper hierarchical mixture models are introduced to infer the extent of individualization. Hierarchical mixtures utilize complementary aspects of mixtures at different levels of the hierarchy. At the first (top) level, a mixture is used to represent homogeneous groups of fingerprints in the population, whereas at the second level, nested mixtures are used as flexible representations of distributions of features from each fingerprint. Inference for hierarchical mixtures is more challenging since the number of unknown mixture components arise in both the first and second levels of the hierarchy. A Bayesian approach based on reversible jump Markov chain Monte Carlo methodology is developed for the inference of all unknown parameters of hierarchical mixtures. The methodology is illustrated on fingerprint images from the NIST database and is used to make inference on fingerprint individuality estimates from this population.

We compare calcium ion signaling (Ca²⁺) between two exposures; the data are present as movies, or, more prosaically, time series of images. This paper describes novel uses of singular value decompositions (SVD) and weighted versions of them (WSVD) to extract the signals from such movies, in a way that is semi-automatic and tuned closely to the actual data and their many complexities. These complexities include the following. First, the images themselves are of no interest: all interest focuses on the behavior of individual cells across time, and thus, the cells need to be segmented in an automated manner. Second, the cells themselves have 100+ pixels, so that they form 100+ curves measured over time, so that data compression is required to extract the features of these curves. Third, some of the pixels in some of the cells are subject to image saturation due to bit depth limits, and this saturation needs to be accounted for if one is to normalize the images in a reasonably unbiased manner. Finally, the Ca²⁺ signals have oscillations or waves that vary with time and these signals need to be extracted. Thus, our aim is to show how to use multiple weighted and standard singular value decompositions to detect, extract and clarify the Ca²⁺ signals. Our signal extraction methods then lead to simple although finely focused statistical methods to compare Ca²⁺ signals across experimental conditions.

The American Community Survey (ACS) provides one-year (1y), three-year (3y) and five-year (5y) multi-year estimates (MYEs) of various demographic and economic variables for each “community,” although the 1y and 3y may not be available for communities with a small population. These survey estimates are not truly measuring the same quantities, since they each cover different time spans. Using some simplistic models, we demonstrate that comparing different period-length MYEs results in spurious conclusions about trend movements. A simple method utilizing weighted averages is presented that reduces the bias inherent in comparing trends of different MYEs. These weighted averages are nonparametric, require only a short span of data, and are designed to preserve polynomial characteristics of the time series that are relevant for trends. The basic method, which only requires polynomial algebra, is outlined and applied to ACS data. In some cases there is an improvement to comparability, although a final verdict must await additional ACS data. We draw the conclusion that MYE data is not comparable across different periods.

This study proposes a nonhomogeneous birth–death model which captures the dynamics of a directly transmitted infectious disease. Our model accounts for an important aspect of observed epidemic data in which only symptomatic infecteds are observed. The nonhomogeneous birth–death process depends on survival distributions of reproduction and removal, which jointly yield an estimate of the effective reproduction number R(t) as a function of epidemic time. We employ the Burr distribution family for the survival functions and, as special cases, proportional rate and accelerated event-time models are also employed for the parameter estimation procedure. As an example, our model is applied to an outbreak of avian influenza (H7N7) in the Netherlands, 2003, confirming that the conditional estimate of R(t) declined below unity for the first time on day 23 since the detection of the index case.

We present a flexible branching process model for cell population dynamics in synchrony/time-series experiments used to study important cellular processes. Its formulation is constructive, based on an accounting of the unique cohorts in the population as they arise and evolve over time, allowing it to be written in closed form. The model can attribute effects to subsets of the population, providing flexibility not available using the models historically applied to these populations. It provides a tool for in silico synchronization of the population and can be used to deconvolve population-level experimental measurements, such as temporal expression profiles. It also allows for the direct comparison of assay measurements made from multiple experiments. The model can be fit either to budding index or DNA content measurements, or both, and is easily adaptable to new forms of data. The ability to use DNA content data makes the model applicable to almost any organism. We describe the model and illustrate its utility and flexibility in a study of cell cycle progression in the yeast Saccharomyces cerevisiae.

There have been several studies of the genome-wide temporal transcriptional program of viruses, based on microarray experiments, which are generally useful in the construction of gene regulation network. It seems that biological interpretations in these studies are directly based on the normalized data and some crude statistics, which provide rough estimates of limited features of the profile and may incur biases. This paper introduces a hierarchical Bayesian shape restricted regression method for making inference on the time course expression of virus genes. Estimates of many salient features of the expression profile like onset time, inflection point, maximum value, time to maximum value, area under curve, etc. can be obtained immediately by this method. Applying this method to a baculovirus microarray time course expression data set, we indicate that many biological questions can be formulated quantitatively and we are able to offer insights into the baculovirus biology.

Support Vector Machines are a widely used classification technique. They are computationally efficient and provide excellent predictions even for high-dimensional data. Moreover, Support Vector Machines are very flexible due to the incorporation of kernel functions. The latter allow to model nonlinearity, but also to deal with nonnumerical data such as protein strings. However, Support Vector Machines can suffer a lot from unclean data containing, for example, outliers or mislabeled observations. Although several outlier detection schemes have been proposed in the literature, the selection of outliers versus nonoutliers is often rather ad hoc and does not provide much insight in the data. In robust multivariate statistics outlier maps are quite popular tools to assess the quality of data under consideration. They provide a visual representation of the data depicting several types of outliers. This paper proposes an outlier map designed for Support Vector Machine classification. The Stahel–Donoho outlyingness measure from multivariate statistics is extended to an arbitrary kernel space. A trimmed version of Support Vector Machines is defined trimming part of the samples with largest outlyingness. Based on this classifier, an outlier map is constructed visualizing data in any type of high-dimensional kernel space. The outlier map is illustrated on 4 biological examples showing its use in exploratory data analysis.

Reliability (survival analysis, to biostatisticians) is a key ingredient for making decisions that mitigate the risk of failure. The other key ingredient is utility. A decision theoretic framework harnesses the two, but to invoke this framework we must distinguish between chance and probability. We describe a functional form for the utility of chance that incorporates all dispositions to risk, and propose a probability of choice model for eliciting this utility. To implement the model, a subject is asked to make a series of binary choices between gambles and certainty. These choices endow a statistical character to the problem of utility elicitation. The workings of our approach are illustrated via a live example involving a military planner. The material is general because it is germane to any situation involving the valuation of chance.

The active field of Functional Data Analysis (about understanding the variation in a set of curves) has been recently extended to Object Oriented Data Analysis, which considers populations of more general objects. A particularly challenging extension of this set of ideas is to populations of tree-structured objects. We develop an analog of Principal Component Analysis for trees, based on the notion of tree-lines, and propose numerically fast (linear time) algorithms to solve the resulting problems to proven optimality. The solutions we obtain are used in the analysis of a data set of 73 individuals, where each data object is a tree of blood vessels in one person’s brain. Our analysis revealed a significant relation between the age of the individuals and their brain vessel structure.

The greatest root distribution occurs everywhere in classical multivariate analysis, but even under the null hypothesis the exact distribution has required extensive tables or special purpose software. We describe a simple approximation, based on the Tracy–Widom distribution, that in many cases can be used instead of tables or software, at least for initial screening. The quality of approximation is studied, and its use illustrated in a variety of settings.

Probe-level microarray data are usually stored in matrices, where the row and column correspond to array and probe, respectively. Scientists routinely summarize each array by a single index as the expression level of each probe set (gene). We examine the adequacy of a unidimensional summary for characterizing the data matrix of each probe set. To do so, we propose a low-rank matrix model for the probe-level intensities, and develop a useful framework for testing the adequacy of unidimensionality against targeted alternatives. This is an interesting statistical problem where inference has to be made based on one data matrix whose entries are not i.i.d. We analyze the asymptotic properties of the proposed test statistics, and use Monte Carlo simulations to assess their small sample performance. Applications of the proposed tests to GeneChip data show that evidence against a unidimensional model is often indicative of practically relevant features of a probe set.

This paper describes a framework for flexible multiple hypothesis testing of autoregressive time series. The modeling approach is Bayesian, though a blend of frequentist and Bayesian reasoning is used to evaluate procedures. Nonparametric characterizations of both the null and alternative hypotheses will be shown to be the key robustification step necessary to ensure reasonable Type-I error performance. The methodology is applied to part of a large database containing up to 50 years of corporate performance statistics on 24,157 publicly traded American companies, where the primary goal of the analysis is to flag companies whose historical performance is significantly different from that expected due to chance.

We present an applied study in cancer genomics for integrating data and inferences from laboratory experiments on cancer cell lines with observational data obtained from human breast cancer studies. The biological focus is on improving understanding of transcriptional responses of tumors to changes in the pH level of the cellular microenvironment. The statistical focus is on connecting experimentally defined biomarkers of such responses to clinical outcome in observational studies of breast cancer patients. Our analysis exemplifies a general strategy for accomplishing this kind of integration across contexts. The statistical methodologies employed here draw heavily on Bayesian sparse factor models for identifying, modularizing and correlating with clinical outcome these signatures of aggregate changes in gene expression. By projecting patterns of biological response linked to specific experimental interventions into observational studies where such responses may be evidenced via variation in gene expression across samples, we are able to define biomarkers of clinically relevant physiological states and outcomes that are rooted in the biology of the original experiment. Through this approach we identify microenvironment-related prognostic factors capable of predicting long term survival in two independent breast cancer datasets. These results suggest possible directions for future laboratory studies, as well as indicate the potential for therapeutic advances though targeted disruption of specific pathway components.

Despite rapid advances in experimental cell biology, the in vivo behavior of hematopoietic stem cells (HSC) cannot be directly observed and measured. Previously we modeled feline hematopoiesis using a two-compartment hidden Markov process that had birth and emigration events in the first compartment. Here we perform Bayesian statistical inference on models which contain two additional events in the first compartment in order to determine if HSC fate decisions are linked to cell division or occur independently. Pareto Optimal Model Assessment approach is used to cross check the estimates from Bayesian inference. Our results show that HSC must divide symmetrically (i.e., produce two HSC daughter cells) in order to maintain hematopoiesis. We then demonstrate that the augmented model that adds asymmetric division events provides a better fit to the competitive transplantation data, and we thus provide evidence that HSC fate determination in vivo occurs both in association with cell division and at a separate point in time. Last we show that assuming each cat has a unique set of parameters leads to either a significant decrease or a nonsignificant increase in model fit, suggesting that the kinetic parameters for HSC are not unique attributes of individual animals, but shared within a species.

Besides serving as prediction models, classification trees are useful for finding important predictor variables and identifying interesting subgroups in the data. These functions can be compromised by weak split selection algorithms that have variable selection biases or that fail to search beyond local main effects at each node of the tree. The resulting models may include many irrelevant variables or select too few of the important ones. Either eventuality can lead to erroneous conclusions. Four techniques to improve the precision of the models are proposed and their effectiveness compared with that of other algorithms, including tree ensembles, on real and simulated data sets.

In linear regression problems with related predictors, it is desirable to do variable selection and estimation by maintaining the hierarchical or structural relationships among predictors. In this paper we propose non-negative garrote methods that can naturally incorporate such relationships defined through effect heredity principles or marginality principles. We show that the methods are very easy to compute and enjoy nice theoretical properties. We also show that the methods can be easily extended to deal with more general regression problems such as generalized linear models. Simulations and real examples are used to illustrate the merits of the proposed methods.

Outlying curves often occur in functional or longitudinal datasets, and can be very influential on parameter estimators and very hard to detect visually. In this article we introduce estimators of the mean and the principal components that are resistant to, and then can be used for detection of, outlying sample trajectories. The estimators are based on reduced-rank t-models and are specifically aimed at sparse and irregularly sampled functional data. The outlier-resistance properties of the estimators and their relative efficiency for noncontaminated data are studied theoretically and by simulation. Applications to the analysis of Internet traffic data and glycated hemoglobin levels in diabetic children are presented.

This paper illustrates novel methods for nonstationary time series modeling along with their applications to selected problems in neuroscience. These methods are semi-parametric in that inferences are derived by combining sequential Bayesian updating with a non-parametric change-point test. As a test statistic, we propose a Kullback–Leibler (KL) divergence between posterior distributions arising from different sets of data. A closed form expression of this statistic is derived for exponential family models, whereas standard Markov chain Monte Carlo output is used to approximate its value and its critical region for more general models. The behavior of one-step ahead predictive distributions under our semi-parametric framework is described analytically for a dynamic linear time series model. Conditions under which our approach reduces to fully parametric state-space modeling are also illustrated. We apply our methods to estimating the functional dynamics of a wide range of neural data, including multi-channel electroencephalogram recordings, longitudinal behavioral experiments and in-vivo multiple spike trains recordings. The estimated dynamics are related to the presentation of visual stimuli, to the evaluation of a learning performance and to changes in the functional connections between neurons over a sequence of experiments.

Rapid developments in geographical information systems (GIS) continue to generate interest in analyzing complex spatial datasets. One area of activity is in creating smoothed disease maps to describe the geographic variation of disease and generate hypotheses for apparent differences in risk. With multiple diseases, a multivariate conditionally autoregressive (MCAR) model is often used to smooth across space while accounting for associations between the diseases. The MCAR, however, imposes complex covariance structures that are difficult to interpret and estimate. This article develops a much simpler alternative approach building upon the techniques of smoothed ANOVA (SANOVA). Instead of simply shrinking effects without any structure, here we use SANOVA to smooth spatial random effects by taking advantage of the spatial structure. We extend SANOVA to cases in which one factor is a spatial lattice, which is smoothed using a CAR model, and a second factor is, for example, type of cancer. Datasets routinely lack enough information to identify the additional structure of MCAR. SANOVA offers a simpler and more intelligible structure than the MCAR while performing as well. We demonstrate our approach with simulation studies designed to compare SANOVA with different design matrices versus MCAR with different priors. Subsequently a cancer-surveillance dataset, describing incidence of 3-cancers in Minnesota’s 87 counties, is analyzed using both approaches, showing the competitiveness of the SANOVA approach.