Ultra-high energy cosmic rays (UHECRs) are atomic nuclei with energies over ten million times energies accessible to human-made particle accelerators. Evidence suggests that they originate from relatively nearby extragalactic sources, but the nature of the sources is unknown. We develop a multilevel Bayesian framework for assessing association of UHECRs and candidate source populations, and Markov chain Monte Carlo algorithms for estimating model parameters and comparing models by computing, via Chib’s method, marginal likelihoods and Bayes factors. We demonstrate the framework by analyzing measurements of 69 UHECRs observed by the Pierre Auger Observatory (PAO) from 2004–2009, using a volume-complete catalog of 17 local active galactic nuclei (AGN) out to 15 megaparsecs as candidate sources. An early portion of the data (“period 1,” with 14 events) was used by PAO to set an energy cut maximizing the anisotropy in period 1; the 69 measurements include this “tuned” subset, and subsequent “untuned” events with energies above the same cutoff. Also, measurement errors are approximately summarized. These factors are problematic for independent analyses of PAO data. Within the context of “standard candle” source models (i.e., with a common isotropic emission rate), and considering only the 55 untuned events, there is no significant evidence favoring association of UHECRs with local AGN vs. an isotropic background. The highest-probability associations are with the two nearest, adjacent AGN, Centaurus A and NGC 4945. If the association model is adopted, the fraction of UHECRs that may be associated is likely nonzero but is well below 50%. Our framework enables estimation of the angular scale for deflection of cosmic rays by cosmic magnetic fields; relatively modest scales of $\approx\!3^{\circ}$ to $30^{\circ}$ are favored. Models that assign a large fraction of UHECRs to a single nearby source (e.g., Centaurus A) are ruled out unless very large deflection scales are specified a priori, and even then they are disfavored. However, including the period 1 data alters the conclusions significantly, and a simulation study supports the idea that the period 1 data are anomalous, presumably due to the tuning. Accurate and optimal analysis of future data will likely require more complete disclosure of the data.

Geomagnetic storms play a critical role in space weather physics with the potential for far reaching economic impacts including power grid outages, air traffic rerouting, satellite damage and GPS disruption. The LFM–MIX is a state-of-the-art coupled magnetospheric–ionospheric model capable of simulating geomagnetic storms. Imbedded in this model are physical equations for turning the magnetohydrodynamic state parameters into energy and flux of electrons entering the ionosphere, involving a set of input parameters. The exact values of these input parameters in the model are unknown, and we seek to quantify the uncertainty about these parameters when model output is compared to observations. The model is available at different fidelities: a lower fidelity which is faster to run, and a higher fidelity but more computationally intense version. Model output and observational data are large spatiotemporal systems; the traditional design and analysis of computer experiments is unable to cope with such large data sets that involve multiple fidelities of model output. We develop an approach to this inverse problem for large spatiotemporal data sets that incorporates two different versions of the physical model. After an initial design, we propose a sequential design based on expected improvement. For the LFM–MIX, the additional run suggested by expected improvement diminishes posterior uncertainty by ruling out a posterior mode and shrinking the width of the posterior distribution. We also illustrate our approach using the Lorenz ‘96 system of equations for a simplified atmosphere, using known input parameters. For the Lorenz ‘96 system, after performing sequential runs based on expected improvement, the posterior mode converges to the true value and the posterior variability is reduced.

Biological structure and function depend on complex regulatory interactions between many genes. A wealth of gene expression data is available from high-throughput genome-wide measurement technologies, but effective gene regulatory network inference methods are still needed. Model-based methods founded on quantitative descriptions of gene regulation are among the most promising, but many such methods rely on simple, local models or on ad hoc inference approaches lacking experimental interpretability. We propose an experimental design and develop an associated statistical method for inferring a gene network by learning a standard quantitative, interpretable, predictive, biophysics-based ordinary differential equation model of gene regulation. We fit the model parameters using gene expression measurements from perturbed steady-states of the system, like those following overexpression or knockdown experiments. Although the original model is nonlinear, our design allows us to transform it into a convex optimization problem by restricting attention to steady-states and using the lasso for parameter selection. Here, we describe the model and inference algorithm and apply them to a synthetic six-gene system, demonstrating that the model is detailed and flexible enough to account for activation and repression as well as synergistic and self-regulation, and the algorithm can efficiently and accurately recover the parameters used to generate the data.

To identify novel dynamic patterns of gene expression, we develop a statistical method to cluster noisy measurements of gene expression collected from multiple replicates at multiple time points, with an unknown number of clusters. We propose a random-effects mixture model coupled with a Dirichlet-process prior for clustering. The mixture model formulation allows for probabilistic cluster assignments. The random-effects formulation allows for attributing the total variability in the data to the sources that are consistent with the experimental design, particularly when the noise level is high and the temporal dependence is not strong. The Dirichlet-process prior induces a prior distribution on partitions and helps to estimate the number of clusters (or mixture components) from the data. We further tackle two challenges associated with Dirichlet-process prior-based methods. One is efficient sampling. We develop a novel Metropolis–Hastings Markov Chain Monte Carlo (MCMC) procedure to sample the partitions. The other is efficient use of the MCMC samples in forming clusters. We propose a two-step procedure for posterior inference, which involves resampling and relabeling, to estimate the posterior allocation probability matrix. This matrix can be directly used in cluster assignments, while describing the uncertainty in clustering. We demonstrate the effectiveness of our model and sampling procedure through simulated data. Applying our method to a real data set collected from Drosophila adult muscle cells after five-minute Notch activation, we identify 14 clusters of different transcriptional responses among 163 differentially expressed genes, which provides novel insights into underlying transcriptional mechanisms in the Notch signaling pathway. The algorithm developed here is implemented in the R package DIRECT, available on CRAN.

Many investigations have used panel methods to study the relationships between fluctuations in economic activity and mortality. A broad consensus has emerged on the overall procyclical nature of mortality: perhaps counter-intuitively, mortality typically rises above its trend during expansions. This consensus has been tarnished by inconsistent reports on the specific age groups and mortality causes involved. We show that these inconsistencies result, in part, from the trend specifications used in previous panel models. Standard econometric panel analysis involves fitting regression models using ordinary least squares, employing standard errors which are robust to temporal autocorrelation. The model specifications include a fixed effect, and possibly a linear trend, for each time series in the panel. We propose alternative methodology based on nonlinear detrending. Applying our methodology on data for the 50 US states from 1980 to 2006, we obtain more precise and consistent results than previous studies. We find procyclical mortality in all age groups. We find clear procyclical mortality due to respiratory disease and traffic injuries. Predominantly procyclical cardiovascular disease mortality and countercyclical suicide are subject to substantial state-to-state variation. Neither cancer nor homicide have significant macroeconomic association.

Causal inference in observational studies typically requires making comparisons between groups that are dissimilar. For instance, researchers investigating the role of a prolonged duration of breastfeeding on child outcomes may be forced to make comparisons between women with substantially different characteristics on average. In the extreme there may exist neighborhoods of the covariate space where there are not sufficient numbers of both groups of women (those who breastfed for prolonged periods and those who did not) to make inferences about those women. This is referred to as lack of common support. Problems can arise when we try to estimate causal effects for units that lack common support, thus we may want to avoid inference for such units. If ignorability is satisfied with respect to a set of potential confounders, then identifying whether, or for which units, the common support assumption holds is an empirical question. However, in the high-dimensional covariate space often required to satisfy ignorability such identification may not be trivial. Existing methods used to address this problem often require reliance on parametric assumptions and most, if not all, ignore the information embedded in the response variable. We distinguish between the concepts of “common support” and “common causal support.” We propose a new approach for identifying common causal support that addresses some of the shortcomings of existing methods. We motivate and illustrate the approach using data from the National Longitudinal Survey of Youth to estimate the effect of breastfeeding at least nine months on reading and math achievement scores at age five or six. We also evaluate the comparative performance of this method in hypothetical examples and simulations where the true treatment effect is known.

An exceedance region is the set of locations in a spatial domain where a process exceeds some threshold. Examples of exceedance regions include areas where ozone concentrations exceed safety standards, there is high risk for tornadoes or floods, or heavy-metal levels are dangerously high. Identifying these regions in a spatial or spatio-temporal setting is an important responsibility in environmental monitoring. Exceedance regions are often estimated by finding the areas where predictions from a statistical model exceed some threshold. Even when estimation error is quantifiable at individual locations, the overall estimation error of the estimated exceedance region is still unknown. A method is presented for constructing a confidence region containing the true exceedance region of a spatio-temporal process at a certain time. The underlying latent process and any measurement error are assumed to be Gaussian. Conventional techniques are used to model the spatio-temporal data, and then conditional simulation is combined with hypothesis testing to create the desired confidence region. A simulation study is used to validate the approach for several levels of spatial and temporal dependence. The methodology is used to identify regions of Oregon having high precipitation levels and also used in comparing climate models and assessing climate change using climate models from the North American Regional Climate Change Assessment Program.

Managing risk at the aggregate level is crucial for banks and financial institutions as required by the Basel III framework. In this paper, we introduce discrete time Bayesian state space models with Poisson measurements to model aggregate mortgage default rate. We discuss parameter updating, filtering, smoothing, forecasting and estimation using Markov chain Monte Carlo methods. In addition, we investigate the dynamic behavior of the default rate and the effects of macroeconomic variables. We illustrate the use of the proposed models using actual U.S. residential mortgage data and discuss insights gained from Bayesian analysis.

When plaintiffs prevail in a discrimination case, a major component of the calculation of economic loss is the length of time they would have been in the higher position had they been treated fairly during the period in which the employer practiced discrimination. This problem is complicated by the fact that one’s eligibility for promotion is subject to termination by retirement and both the promotion and retirement processes may be affected by discriminatory practices. This semi-competing risk setup is decomposed into a retirement process and a promotion process among the employees. Predictions for the purpose of compensation are made by utilizing the expected promotion and retirement probabilities of similarly qualified members of the nondiscriminated group. The restricted mean durations of three periods are estimated—the time an employee would be at the lower position, at the higher level and in retirement. The asymptotic properties of the estimators are presented and examined through simulation studies. The proposed restricted mean job duration estimators are shown to be robust in the presence of an independent frailty term. Data from the reverse discrimination case, Alexander v. Milwaukee, where White-male lieutenants were discriminated in promotion to captain are reanalyzed. While the appellate court upheld liability, it reversed the original damage calculations, which heavily depended on the time a plaintiff would have been in each position. The results obtained by the proposed method are compared to those made at the first trial. Substantial differences in both directions are observed.

Evaluating the overall ability of players in the National Hockey League (NHL) is a difficult task. Existing methods such as the famous “plus/minus” statistic have many shortcomings. Standard linear regression methods work well when player substitutions are relatively uncommon and scoring events are relatively common, such as in basketball, but as neither of these conditions exists for hockey, we use an approach that embraces the unique characteristics of the sport. We model the scoring rate for each team as its own semi-Markov process, with hazard functions for each process that depend on the players on the ice. This method yields offensive and defensive player ability ratings which take into account quality of teammates and opponents, the game situation, and other desired factors, that themselves have a meaningful interpretation in terms of game outcomes. Additionally, since the number of parameters in this model can be quite large, we make use of two different shrinkage methods depending on the question of interest: full Bayesian hierarchical models that partially pool parameters according to player position, and penalized maximum likelihood estimation to select a smaller number of parameters that stand out as being substantially different from average. We apply the model to all five-on-five (full-strength) situations for games in five NHL seasons.

While the presence of clustering in crime and security event data is well established, the mechanism(s) by which clustering arises is not fully understood. Both contagion models and history independent correlation models are applied, but not simultaneously. In an attempt to disentangle contagion from other types of correlation, we consider a Hawkes process with background rate driven by a log Gaussian Cox process. Our inference methodology is an efficient Metropolis adjusted Langevin algorithm for filtering of the intensity and estimation of the model parameters. We apply the methodology to property and violent crime data from Chicago, terrorist attack data from Northern Ireland and Israel, and civilian casualty data from Iraq. For each data set we quantify the uncertainty in the levels of contagion vs. history independent correlation.

We develop a robust regularized singular value decomposition (RobRSVD) method for analyzing two-way functional data. The research is motivated by the application of modeling human mortality as a smooth two-way function of age group and year. The RobRSVD is formulated as a penalized loss minimization problem where a robust loss function is used to measure the reconstruction error of a low-rank matrix approximation of the data, and an appropriately defined two-way roughness penalty function is used to ensure smoothness along each of the two functional domains. By viewing the minimization problem as two conditional regularized robust regressions, we develop a fast iterative reweighted least squares algorithm to implement the method. Our implementation naturally incorporates missing values. Furthermore, our formulation allows rigorous derivation of leave-one-row/column-out cross-validation and generalized cross-validation criteria, which enable computationally efficient data-driven penalty parameter selection. The advantages of the new robust method over nonrobust ones are shown via extensive simulation studies and the mortality rate application.

Classical time series models have serious difficulties in modeling and forecasting the enormous fluctuations of electricity spot prices. Markov regime switch models belong to the most often used models in the electricity literature. These models try to capture the fluctuations of electricity spot prices by using different regimes, each with its own mean and covariance structure. Usually one regime is dedicated to moderate prices and another is dedicated to high prices. However, these models show poor performance and there is no theoretical justification for this kind of classification. The merit order model, the most important micro-economic pricing model for electricity spot prices, however, suggests a continuum of mean levels with a functional dependence on electricity demand.

We propose a new statistical perspective on modeling and forecasting electricity spot prices that accounts for the merit order model. In a first step, the functional relation between electricity spot prices and electricity demand is modeled by daily price-demand functions. In a second step, we parameterize the series of daily price-demand functions using a functional factor model. The power of this new perspective is demonstrated by a forecast study that compares our functional factor model with two established classical time series models as well as two alternative functional data models.

Global climate models aim to reproduce physical processes on a global scale and predict quantities such as temperature given some forcing inputs. We consider climate ensembles made of collections of such runs with different initial conditions and forcing scenarios. The purpose of this work is to show how the simulated temperatures in the ensemble can be reproduced (emulated) with a global space/time statistical model that addresses the issue of capturing nonstationarities in latitude more effectively than current alternatives in the literature. The model we propose leads to a computationally efficient estimation procedure and, by exploiting the gridded geometry of the data, we can fit massive data sets with millions of simulated data within a few hours. Given a training set of runs, the model efficiently emulates temperature for very different scenarios and therefore is an appealing tool for impact assessment.

In political campaigning substantial resources are spent on voter mobilization, that is, on identifying and influencing as many people as possible to vote. Campaigns use statistical tools for deciding whom to target (“microtargeting”). In this paper we describe a nonpartisan campaign that aims at increasing overall turnout using the example of the 2004 US presidential election. Based on a real data set of 19,634 eligible voters from Ohio, we introduce a modern statistical framework well suited for carrying out the main tasks of voter targeting in a single sweep: predicting an individual’s turnout (or support) likelihood for a particular cause, party or candidate as well as data-driven voter segmentation. Our framework, which we refer to as LORET (for LOgistic REgression Trees), contains standard methods such as logistic regression and classification trees as special cases and allows for a synthesis of both techniques. For our case study, we explore various LORET models with different regressors in the logistic model components and different partitioning variables in the tree components; we analyze them in terms of their predictive accuracy and compare the effect of using the full set of available variables against using only a limited amount of information. We find that augmenting a standard set of variables (such as age and voting history) with additional predictor variables (such as the household composition in terms of party affiliation) clearly improves predictive accuracy. We also find that LORET models based on tree induction beat the unpartitioned models. Furthermore, we illustrate how voter segmentation arises from our framework and discuss the resulting profiles from a targeting point of view.

Two different approaches to analysis of data from diagnostic biomarker studies are commonly employed. Logistic regression is used to fit models for probability of disease given marker values, while ROC curves and risk distributions are used to evaluate classification performance. In this paper we present a method that simultaneously accomplishes both tasks. The key step is to standardize markers relative to the nondiseased population before including them in the logistic regression model. Among the advantages of this method are the following: (i) ensuring that results from regression and performance assessments are consistent with each other; (ii) allowing covariate adjustment and covariate effects on ROC curves to be handled in a familiar way, and (iii) providing a mechanism to incorporate important assumptions about structure in the ROC curve into the fitted risk model. We develop the method in detail for the problem of combining biomarker data sets derived from multiple studies, populations or biomarker measurement platforms, when ROC curves are similar across data sources. The methods are applicable to both cohort and case–control sampling designs. The data set motivating this application concerns Prostate Cancer Antigen 3 (PCA3) for diagnosis of prostate cancer in patients with or without previous negative biopsy where the ROC curves for PCA3 are found to be the same in the two populations. The estimated constrained maximum likelihood and empirical likelihood estimators are derived. The estimators are compared in simulation studies and the methods are illustrated with the PCA3 data set.

Multivariate time series (MTS) data such as time course gene expression data in genomics are often collected to study the dynamic nature of the systems. These data provide important information about the causal dependency among a set of random variables. In this paper, we introduce a computationally efficient algorithm to learn directed acyclic graphs (DAGs) based on MTS data, focusing on learning the local structure of a given target variable. Our algorithm is based on learning all parents (P), all children (C) and some descendants (D) (PCD) iteratively, utilizing the time order of the variables to orient the edges. This time series PCD-PCD algorithm (tsPCD-PCD) extends the previous PCD-PCD algorithm to dependent observations and utilizes composite likelihood ratio tests (CLRTs) for testing the conditional independence. We present the asymptotic distribution of the CLRT statistic and show that the tsPCD-PCD is guaranteed to recover the true DAG structure when the faithfulness condition holds and the tests correctly reject the null hypotheses. Simulation studies show that the CLRTs are valid and perform well even when the sample sizes are small. In addition, the tsPCD-PCD algorithm outperforms the PCD-PCD algorithm in recovering the local graph structures. We illustrate the algorithm by analyzing a time course gene expression data related to mouse T-cell activation.

The Atmospheric Radiation Measurement program is a U.S. Department of Energy project that collects meteorological observations at several locations around the world in order to study how weather processes affect global climate change. As one of its initiatives, it operates a set of fixed but irregularly-spaced monitoring facilities in the Southern Great Plains region of the U.S. We describe methods for interpolating temperature records from these fixed facilities to locations at which no observations were made, which can be useful when values are required on a spatial grid. We interpolate by conditionally simulating from a fitted nonstationary Gaussian process model that accounts for the time-varying statistical characteristics of the temperatures, as well as the dependence on solar radiation. The model is fit by maximizing an approximate likelihood, and the conditional simulations result in well-calibrated confidence intervals for the predicted temperatures. We also describe methods for handling spatial–temporal jumps in the data to interpolate a slow-moving cold front.

We consider mark–recapture–recovery (MRR) data of animals where the model parameters are a function of individual time-varying continuous covariates. For such covariates, the covariate value is unobserved if the corresponding individual is unobserved, in which case the survival probability cannot be evaluated. For continuous-valued covariates, the corresponding likelihood can only be expressed in the form of an integral that is analytically intractable and, to date, no maximum likelihood approach that uses all the information in the data has been developed. Assuming a first-order Markov process for the covariate values, we accomplish this task by formulating the MRR setting in a state-space framework and considering an approximate likelihood approach which essentially discretizes the range of covariate values, reducing the integral to a summation. The likelihood can then be efficiently calculated and maximized using standard techniques for hidden Markov models. We initially assess the approach using simulated data before applying to real data relating to Soay sheep, specifying the survival probability as a function of body mass. Models that have previously been suggested for the corresponding covariate process are typically of the form of diffusive random walks. We consider an alternative nondiffusive AR(1)-type model which appears to provide a significantly better fit to the Soay sheep data.

Simulating sample correlation matrices is important in many areas of statistics. Approaches such as generating Gaussian data and finding their sample correlation matrix or generating random uniform $[-1,1]$ deviates as pairwise correlations both have drawbacks. We develop an algorithm for adding noise, in a highly controlled manner, to general correlation matrices. In many instances, our method yields results which are superior to those obtained by simply simulating Gaussian data. Moreover, we demonstrate how our general algorithm can be tailored to a number of different correlation models. Using our results with a few different applications, we show that simulating correlation matrices can help assess statistical methodology.

Data used to assess acute health effects from air pollution typically have good temporal but poor spatial resolution or the opposite. A modified longitudinal model was developed that sought to improve resolution in both domains by bringing together data from three sources to estimate daily levels of nitrogen dioxide ($\mathrm{NO} _{2}$) at a geographic location. Monthly $\mathrm{NO} _{2}$ measurements at 316 sites were made available by the Study of Traffic, Air quality and Respiratory health (STAR). Four US Environmental Protection Agency monitoring stations have hourly measurements of $\mathrm{NO} _{2}$. Finally, the Connecticut Department of Transportation provides data on traffic density on major roadways, a primary contributor to $\mathrm{NO} _{2}$ pollution. Inclusion of a traffic variable improved performance of the model, and it provides a method for estimating exposure at points that do not have direct measurements of the outcome. This approach can be used to estimate daily variation in levels of $\mathrm{NO} _{2}$ over a region.

Random effects are implemented for aster models using two approximations taken from Breslow and Clayton [J. Amer. Statist. Assoc. 88 (1993) 9–25]. Random effects are analytically integrated out of the Laplace approximation to the complete data log likelihood, giving a closed-form expression for an approximate missing data log likelihood. Third and higher derivatives of the complete data log likelihood with respect to the random effects are ignored, giving a closed-form expression for second derivatives of the approximate missing data log likelihood, hence approximate observed Fisher information. This method is applicable to any exponential family random effects model. It is implemented in the CRAN package aster (R Core Team [R: A Language and Environment for Statistical Computing (2012) R Foundation for Statistical Computing], Geyer [R package aster (2012) http://cran.r-project.org/package=aster]). Applications are analyses of local adaptation in the invasive California wild radish (Raphanus sativus) and the slender wild oat (Avena barbata) and of additive genetic variance for fitness in the partridge pea (Chamaecrista fasciculata).

It is generally believed that bias is minimized in well-controlled randomized clinical trials. However, bias can arise in active controlled noninferiority trials because the inference relies on a previously estimated effect size obtained from a historical trial that may have been conducted for a different population. By implementing a likelihood reweighting method through propensity scoring, a study designed to estimate a treatment effect in one trial population can be used to estimate the treatment effect size in a different target population. We illustrate this method in active controlled noninferiority trials, although it can also be used in other types of studies, such as historically controlled trials, meta-analyses, and comparative effectiveness analyses.

This paper is concerned with statistical methods for the segmental classification of linear sequence data where the task is to segment and classify the data according to an underlying hidden discrete state sequence. Such analysis is commonplace in the empirical sciences including genomics, finance and speech processing. In particular, we are interested in answering the following question: given data $y$ and a statistical model $\pi(x,y)$ of the hidden states $x$, what should we report as the prediction $\hat{x}$ under the posterior distribution $\pi(x|y)$? That is, how should you make a prediction of the underlying states? We demonstrate that traditional approaches such as reporting the most probable state sequence or most probable set of marginal predictions can give undesirable classification artefacts and offer limited control over the properties of the prediction. We propose a decision theoretic approach using a novel class of Markov loss functions and report $\hat{x}$ via the principle of minimum expected loss (maximum expected utility). We demonstrate that the sequence of minimum expected loss under the Markov loss function can be enumerated exactly using dynamic programming methods and that it offers flexibility and performance improvements over existing techniques. The result is generic and applicable to any probabilistic model on a sequence, such as Hidden Markov models, change point or product partition models.