Quantities with right-skewed distributions are ubiquitous in complex social systems, including political conflict, economics and social networks, and these systems sometimes produce extremely large events. For instance, the 9/11 terrorist events produced nearly 3000 fatalities, nearly six times more than the next largest event. But, was this enormous loss of life statistically unlikely given modern terrorism’s historical record? Accurately estimating the probability of such an event is complicated by the large fluctuations in the empirical distribution’s upper tail. We present a generic statistical algorithm for making such estimates, which combines semi-parametric models of tail behavior and a nonparametric bootstrap. Applied to a global database of terrorist events, we estimate the worldwide historical probability of observing at least one 9/11-sized or larger event since 1968 to be 11–35%. These results are robust to conditioning on global variations in economic development, domestic versus international events, the type of weapon used and a truncated history that stops at 1998. We then use this procedure to make a data-driven statistical forecast of at least one similar event over the next decade.

A first step in exploring population structure in crop plants and other organisms is to define the number of subpopulations that exist for a given data set. The genetic marker data sets being generated have become increasingly large over time and commonly are of the high-dimension, low sample size (HDLSS) situation. An algorithm for deciding the number of clusters is proposed, and is validated on simulated data sets varying in both the level of structure and the number of clusters covering the range of variation observed empirically. The algorithm was then tested on six empirical data sets across three small grain species. The algorithm uses bootstrapping, three methods of clustering, and defines the optimum number of clusters based on a common criterion, the Hubert’s gamma statistic. Validation on simulated sets coupled with testing on empirical sets suggests that the algorithm can be used for a wide variety of genetic data sets.

Statistical modeling of presence-only data has attracted much recent attention in the ecological literature, leading to a proliferation of methods, including the inhomogeneous Poisson process (IPP) model, maximum entropy (Maxent) modeling of species distributions and logistic regression models. Several recent articles have shown the close relationships between these methods. We explain why the IPP intensity function is a more natural object of inference in presence-only studies than occurrence probability (which is only defined with reference to quadrat size), and why presence-only data only allows estimation of relative, and not absolute intensity of species occurrence.

All three of the above techniques amount to parametric density estimation under the same exponential family model (in the case of the IPP, the fitted density is multiplied by the number of presence records to obtain a fitted intensity). We show that IPP and Maxent give the exact same estimate for this density, but logistic regression in general yields a different estimate in finite samples. When the model is misspecified—as it practically always is—logistic regression and the IPP may have substantially different asymptotic limits with large data sets. We propose “infinitely weighted logistic regression,” which is exactly equivalent to the IPP in finite samples. Consequently, many already-implemented methods extending logistic regression can also extend the Maxent and IPP models in directly analogous ways using this technique.

Previous studies on event-related functional magnetic resonance imaging experimental designs are primarily based on linear models, in which a known shape of the hemodynamic response function (HRF) is assumed. However, the HRF shape is usually uncertain at the design stage. To address this issue, we consider a nonlinear model to accommodate a wide spectrum of feasible HRF shapes, and propose efficient approaches for obtaining maximin and maximin-efficient designs. Our approaches involve a reduction in the parameter space and a search algorithm that helps to efficiently search over a restricted class of designs for good designs. The obtained designs are compared with traditional designs widely used in practice. We also demonstrate the usefulness of our approaches via a motivating example.

A central challenge in computational modeling of dynamic biological systems is parameter inference from experimental time course measurements. However, one would not only like to infer kinetic parameters but also study their variability from cell to cell. Here we focus on the case where single-cell fluorescent protein imaging time series data are available for a population of cells. Based on van Kampen’s linear noise approximation, we derive a dynamic state space model for molecular populations which is then extended to a hierarchical model. This model has potential to address the sources of variability relevant to single-cell data, namely, intrinsic noise due to the stochastic nature of the birth and death processes involved in reactions and extrinsic noise arising from the cell-to-cell variation of kinetic parameters. In order to infer such a model from experimental data, one must also quantify the measurement process where one has to allow for nonmeasurable molecular species as well as measurement noise of unknown level and variance. The availability of multiple single-cell time series data here provides a unique testbed to fit such a model and quantify these different sources of variation from experimental data.

A common problem faced by statistical institutes is that data may be missing from collected data sets. The typical way to overcome this problem is to impute the missing data. The problem of imputing missing data is complicated by the fact that statistical data often have to satisfy certain edit rules and that values of variables across units sometimes have to sum up to known totals. For numerical data, edit rules are most often formulated as linear restrictions on the variables. For example, for data on enterprises edit rules could be that the profit and costs of an enterprise should sum up to its turnover and that the turnover should be at least zero. The totals of some variables across units may already be known from administrative data (e.g., turnover from a tax register) or estimated from other sources. Standard imputation methods for numerical data as described in the literature generally do not take such edit rules and totals into account. In this article we describe algorithms for imputing missing numerical data that take edit restrictions into account and ensure that sums are calibrated to known totals. These algorithms are based on a sequential regression approach that uses regression predictions to impute the variables one by one. To assess the performance of the imputation methods, a simulation study is carried out as well as an evaluation study based on a real data set.

In computational neuroscience, it is important to estimate well the proportion of signal variance in the total variance of neural activity measurements. This explainable variance measure helps neuroscientists assess the adequacy of predictive models that describe how images are encoded in the brain. Complicating the estimation problem are strong noise correlations, which may confound the neural responses corresponding to the stimuli. If not properly taken into account, the correlations could inflate the explainable variance estimates and suggest false possible prediction accuracies.

We propose a novel method to estimate the explainable variance in functional MRI (fMRI) brain activity measurements when there are strong correlations in the noise. Our shuffle estimator is nonparametric, unbiased, and built upon the random effect model reflecting the randomization in the fMRI data collection process. Leveraging symmetries in the measurements, our estimator is obtained by appropriately permuting the measurement vector in such a way that the noise covariance structure is intact but the explainable variance is changed after the permutation. This difference is then used to estimate the explainable variance. We validate the properties of the proposed method in simulation experiments. For the image-fMRI data, we show that the shuffle estimates can explain the variation in prediction accuracy for voxels within the primary visual cortex (V1) better than alternative parametric methods.

This study presents a Bayesian parametric model for the purpose of estimating the extreme load on a wind turbine. The extreme load is the highest stress level imposed on a turbine structure that the turbine would experience during its service lifetime. A wind turbine should be designed to resist such a high load to avoid catastrophic structural failures. To assess the extreme load, turbine structural responses are evaluated by conducting field measurement campaigns or performing aeroelastic simulation studies. In general, data obtained in either case are not sufficient to represent various loading responses under all possible weather conditions. An appropriate extrapolation is necessary to characterize the structural loads in a turbine’s service life. This study devises a Bayesian spline method for this extrapolation purpose, using load data collected in a period much shorter than a turbine’s service life. The spline method is applied to three sets of turbine’s load response data to estimate the corresponding extreme loads at the roots of the turbine blades. Compared to the current industry practice, the spline method appears to provide better extreme load assessment.

Estimates of the effects of treatment on cost from observational studies are subject to bias if there are unmeasured confounders. It is therefore advisable in practice to assess the potential magnitude of such biases. We derive a general adjustment formula for loglinear models of mean cost and explore special cases under plausible assumptions about the distribution of the unmeasured confounder. We assess the performance of the adjustment by simulation, in particular, examining robustness to a key assumption of conditional independence between the unmeasured and measured covariates given the treatment indicator. We apply our method to SEER-Medicare cost data for a stage II/III muscle-invasive bladder cancer cohort. We evaluate the costs for radical cystectomy vs. combined radiation/chemotherapy, and find that the significance of the treatment effect is sensitive to plausible unmeasured Bernoulli, Poisson and Gamma confounders.

In June of 2010 access via robotic means was obtained to a tomb adjacent to the one studied in Feuerverger [Ann. Appl. Stat. 2 (2008) 3–54]. In this update, we lay out and attempt to interpret the remarkable findings from this second tomb and comment on the statistical and scientific significance of these new data and of their possible inferential connections to the data from the first tomb. Readers are then invited to formulate their own conclusions.

The primary goal of randomized trials is to compare the effects of different interventions on some outcome of interest. In addition to the treatment assignment and outcome, data on baseline covariates, such as demographic characteristics or biomarker measurements, are typically collected. Incorporating such auxiliary covariates in the analysis of randomized trials can increase power, but questions remain about how to preserve type I error when incorporating such covariates in a flexible way, particularly when the number of randomized units is small. Using the Young Citizens study, a cluster-randomized trial of an educational intervention to promote HIV awareness, we compare several methods to evaluate intervention effects when baseline covariates are incorporated adaptively. To ascertain the validity of the methods shown in small samples, extensive simulation studies were conducted. We demonstrate that randomization inference preserves type I error under model selection while tests based on asymptotic theory may yield invalid results. We also demonstrate that covariate adjustment generally increases power, except at extremely small sample sizes using liberal selection procedures. Although shown within the context of HIV prevention research, our conclusions have important implications for maximizing efficiency and robustness in randomized trials with small samples across disciplines.

A major practical impediment when implementing adaptive dose-finding designs is that the toxicity outcome used by the decision rules may not be observed shortly after the initiation of the treatment. To address this issue, we propose the data augmentation continual reassessment method (DA-CRM) for dose finding. By naturally treating the unobserved toxicities as missing data, we show that such missing data are nonignorable in the sense that the missingness depends on the unobserved outcomes. The Bayesian data augmentation approach is used to sample both the missing data and model parameters from their posterior full conditional distributions. We evaluate the performance of the DA-CRM through extensive simulation studies and also compare it with other existing methods. The results show that the proposed design satisfactorily resolves the issues related to late-onset toxicities and possesses desirable operating characteristics: treating patients more safely and also selecting the maximum tolerated dose with a higher probability. The new DA-CRM is illustrated with two phase I cancer clinical trials.

We consider penalized estimation in hidden Markov models (HMMs) with multivariate Normal observations. In the moderate-to-large dimensional setting, estimation for HMMs remains challenging in practice, due to several concerns arising from the hidden nature of the states. We address these concerns by $\ell_{1}$-penalization of state-specific inverse covariance matrices. Penalized estimation leads to sparse inverse covariance matrices which can be interpreted as state-specific conditional independence graphs. Penalization is nontrivial in this latent variable setting; we propose a penalty that automatically adapts to the number of states $K$ and the state-specific sample sizes and can cope with scaling issues arising from the unknown states. The methodology is adaptive and very general, applying in particular to both low- and high-dimensional settings without requiring hand tuning. Furthermore, our approach facilitates exploration of the number of states $K$ by coupling estimation for successive candidate values $K$. Empirical results on simulated examples demonstrate the effectiveness of the proposed approach. In a challenging real data example from genome biology, we demonstrate the ability of our approach to yield gains in predictive power and to deliver richer estimates than existing methods.

Understanding the mechanisms that allow biological species to co-occur is of great interest to ecologists. Here we investigate the factors that influence co-occurrence of members of the genus Protea in the Cape Floristic Region of southwestern Africa, a global hot spot of biodiversity. Due to the binomial nature of our response, a critical issue is to choose appropriate link functions for the regression model. In this paper we propose a new family of flexible link functions for modeling binomial response data. By introducing a power parameter into the cumulative distribution function (c.d.f.) corresponding to a symmetric link function and its mirror reflection, greater flexibility in skewness can be achieved in both positive and negative directions. Through simulated data sets and analysis of the Protea co-occurrence data, we show that the proposed link function is quite flexible and performs better against link misspecification than standard link functions.

Time-course high-throughput gene expression data are emerging in genomic and translational medicine. Extracting interesting time-course patterns from a patient cohort can provide biological insights for further clinical research and patient treatment. We propose principal trend analysis (PTA) to extract principal trends of time-course gene expression data from a group of patients, and identify genes that make dominant contributions to the principal trends. Through simulations, we demonstrate the utility of PTA for dimension reduction, time-course signal recovery and feature selection with high-dimensional data. Moreover, PTA derives new insights in real biological and clinical research. We demonstrate the usefulness of PTA by applying it to longitudinal gene expression data of a circadian regulation system and burn patients. These applications show that PTA can extract interesting time-course trends with biological significance, which helps the understanding of biological mechanisms of circadian regulation systems as well as the recovery of burn patients. Overall, the proposed PTA approach will benefit the genomic medicine research. Our method is implemented into an R-package: PTA (Principal Trend Analysis).

Because of the decreasing cost and high digital resolution, next-generation sequencing (NGS) is expected to replace the traditional hybridization-based microarray technology. For genetics study, the first-step analysis of NGS data is often to identify genomic variants among sequenced samples. Several statistical models and tests have been developed for variant calling in NGS study. The existing approaches, however, are based on either conventional Bayesian or frequentist methods, which are unable to address the multiplicity and testing efficiency issues simultaneously. In this paper, we derive an optimal empirical Bayes testing procedure to detect variants for NGS study. We utilize the empirical Bayes technique to exploit the across-site information among many testing sites in NGS data. We prove that our testing procedure is valid and optimal in the sense of rejecting the maximum number of nonnulls while the Bayesian false discovery rate is controlled at a given nominal level. We show by both simulation studies and real data analysis that our testing efficiency can be greatly enhanced over the existing frequentist approaches that fail to pool and utilize information across the multiple testing sites.

The main objective of accelerated life tests (ALTs) is to predict fraction failings of products in the field. However, there are often discrepancies between the predicted fraction failing from the lab testing data and that from the field failure data, due to the yet unobserved heterogeneities in usage and operating conditions. Most previous research on ALT planning and data analysis ignores the discrepancies, resulting in inferior test plans and biased predictions. In this paper we model the heterogeneous environments together with their effects on the product failures as a frailty term to link the lab failure time distribution and field failure time distribution of a product. We show that in the presence of the heterogeneous operating conditions, the hazard rate function of the field failure time distribution exhibits a range of shapes. Statistical inference procedure for the frailty models is developed when both the ALT data and the field failure data are available. Based on the frailty models, optimal ALT plans aimed at predicting the field failure time distribution are obtained. The developed methods are demonstrated through a real life example.

Motivated by the increasing use of and rapid changes in array technologies, we consider the prediction problem of fitting a linear regression relating a continuous outcome $Y$ to a large number of covariates $\mathbf{X}$, for example, measurements from current, state-of-the-art technology. For most of the samples, only the outcome $Y$ and surrogate covariates, $\mathbf{W}$, are available. These surrogates may be data from prior studies using older technologies. Owing to the dimension of the problem and the large fraction of missing information, a critical issue is appropriate shrinkage of model parameters for an optimal bias-variance trade-off. We discuss a variety of fully Bayesian and Empirical Bayes algorithms which account for uncertainty in the missing data and adaptively shrink parameter estimates for superior prediction. These methods are evaluated via a comprehensive simulation study. In addition, we apply our methods to a lung cancer data set, predicting survival time ($Y$) using qRT-PCR ($\mathbf{X}$) and microarray ($\mathbf{W}$) measurements.

We propose a general and flexible capture–recapture model in continuous time. Our model incorporates time-heterogeneity, observed and unobserved individual heterogeneity, and behavioral response to capture. Behavioral response can possibly have a delayed onset and a finite-time memory. Estimation of the population size is based on the conditional likelihood after use of the EM algorithm. We develop an application to the estimation of the number of adult cannabinoid users in Italy.

Continuous-time linear birth–death-immigration (BDI) processes are frequently used in ecology and epidemiology to model stochastic dynamics of the population of interest. In clinical settings, multiple birth–death processes can describe disease trajectories of individual patients, allowing for estimation of the effects of individual covariates on the birth and death rates of the process. Such estimation is usually accomplished by analyzing patient data collected at unevenly spaced time points, referred to as panel data in the biostatistics literature. Fitting linear BDI processes to panel data is a nontrivial optimization problem because birth and death rates can be functions of many parameters related to the covariates of interest. We propose a novel expectation–maximization (EM) algorithm for fitting linear BDI models with covariates to panel data. We derive a closed-form expression for the joint generating function of some of the BDI process statistics and use this generating function to reduce the E-step of the EM algorithm, as well as calculation of the Fisher information, to one-dimensional integration. This analytical technique yields a computationally efficient and robust optimization algorithm that we implemented in an open-source R package. We apply our method to DNA fingerprinting of Mycobacterium tuberculosis, the causative agent of tuberculosis, to study intrapatient time evolution of IS6110 copy number, a genetic marker frequently used during estimation of epidemiological clusters of Mycobacterium tuberculosis infections. Our analysis reveals previously undocumented differences in IS6110 birth–death rates among three major lineages of Mycobacterium tuberculosis, which has important implications for epidemiologists that use IS6110 for DNA fingerprinting of Mycobacterium tuberculosis.

The causal effect of a randomized job training program, the JOBS II study, on trainees’ depression is evaluated. Principal stratification is used to deal with noncompliance to the assigned treatment. Due to the latent nature of the principal strata, strong structural assumptions are often invoked to identify principal causal effects. Alternatively, distributional assumptions may be invoked using a model-based approach. These often lead to weakly identified models with substantial regions of flatness in the posterior distribution of the causal effects. Information on multiple outcomes is routinely collected in practice, but is rarely used to improve inference. This article develops a Bayesian approach to exploit multivariate outcomes to sharpen inferences in weakly identified principal stratification models. We show that inference for the causal effect on depression is significantly improved by using the re-employment status as a secondary outcome in the JOBS II study. Simulation studies are also performed to illustrate the potential gains in the estimation of principal causal effects from jointly modeling more than one outcome. This approach can also be used to assess plausibility of structural assumptions and sensitivity to deviations from these structural assumptions. Two model checking procedures via posterior predictive checks are also discussed.

Multivariate data that combine binary, categorical, count and continuous outcomes are common in the social and health sciences. We propose a semiparametric Bayesian latent variable model for multivariate data of arbitrary type that does not require specification of conditional distributions. Drawing on the extended rank likelihood method by Hoff [Ann. Appl. Stat. 1 (2007) 265–283], we develop a semiparametric approach for latent variable modeling with mixed outcomes and propose associated Markov chain Monte Carlo estimation methods. Motivated by cognitive testing data, we focus on bifactor models, a special case of factor analysis. We employ our semiparametric Bayesian latent variable model to investigate the association between cognitive outcomes and MRI-measured regional brain volumes.

Signal averaging is the process that consists in computing a mean shape from a set of noisy signals. In the presence of geometric variability in time in the data, the usual Euclidean mean of the raw data yields a mean pattern that does not reflect the typical shape of the observed signals. In this setting, it is necessary to use alignment techniques for a precise synchronization of the signals, and then to average the aligned data to obtain a consistent mean shape. In this paper, we study the numerical performances of Fréchet means of curves which are extensions of the usual Euclidean mean to spaces endowed with non-Euclidean metrics. This yields a new algorithm for signal averaging and for the estimation of the time variability of a set of signals. We apply this approach to the analysis of heartbeats from ECG records.

The main focus of this work is on developing models for the activity profile of a terrorist group, detecting sudden spurts and downfalls in this profile, and, in general, tracking it over a period of time. Toward this goal, a $d$-state hidden Markov model (HMM) that captures the latent states underlying the dynamics of the group and thus its activity profile is developed. The simplest setting of $d=2$ corresponds to the case where the dynamics are coarsely quantized as Active and Inactive, respectively. A state estimation strategy that exploits the underlying HMM structure is then developed for spurt detection and tracking. This strategy is shown to track even nonpersistent changes that last only for a short duration at the cost of learning the underlying model. Case studies with real terrorism data from open-source databases are provided to illustrate the performance of the proposed methodology.

We consider the problem of quantifying temporal coordination between multiple high-dimensional responses. We introduce a family of multi-way stochastic blockmodels suited for this problem, which avoids preprocessing steps such as binning and thresholding commonly adopted for this type of data, in biology. We develop two inference procedures based on collapsed Gibbs sampling and variational methods. We provide a thorough evaluation of the proposed methods on simulated data, in terms of membership and blockmodel estimation, predictions out-of-sample and run-time. We also quantify the effects of censoring procedures such as binning and thresholding on the estimation tasks. We use these models to carry out an empirical analysis of the functional mechanisms driving the coordination between gene expression and metabolite concentrations during carbon and nitrogen starvation, in S. cerevisiae.