We introduce a new estimator for the vector of coefficients $\beta$ in the linear model $y=X\beta+z$, where $X$ has dimensions $n\times p$ with $p$ possibly larger than $n$. SLOPE, short for Sorted L-One Penalized Estimation, is the solution to

\[\min_{b\in\mathbb{R}^{p}}\frac{1}{2}\Vert y-Xb\Vert_{\ell_{2}}^{2}+\lambda_{1}\vert b\vert_{(1)}+\lambda_{2}\vert b\vert_{(2)}+\cdots+\lambda_{p}\vert b\vert_{(p)},\] where $\lambda_{1}\ge\lambda_{2}\ge\cdots\ge\lambda_{p}\ge0$ and $\vert b\vert_{(1)}\ge\vert b\vert_{(2)}\ge\cdots\ge\vert b\vert_{(p)}$ are the decreasing absolute values of the entries of $b$. This is a convex program and we demonstrate a solution algorithm whose computational complexity is roughly comparable to that of classical $\ell_{1}$ procedures such as the Lasso. Here, the regularizer is a sorted $\ell_{1}$ norm, which penalizes the regression coefficients according to their rank: the higher the rank—that is, stronger the signal—the larger the penalty. This is similar to the Benjamini and Hochberg [J. Roy. Statist. Soc. Ser. B 57 (1995) 289–300] procedure (BH) which compares more significant $p$-values with more stringent thresholds. One notable choice of the sequence $\{\lambda_{i}\}$ is given by the BH critical values $\lambda_{\mathrm{BH}}(i)=z(1-i\cdot q/2p)$, where $q\in(0,1)$ and $z(\alpha)$ is the quantile of a standard normal distribution. SLOPE aims to provide finite sample guarantees on the selected model; of special interest is the false discovery rate (FDR), defined as the expected proportion of irrelevant regressors among all selected predictors. Under orthogonal designs, SLOPE with $\lambda_{\mathrm{BH}}$ provably controls FDR at level $q$. Moreover, it also appears to have appreciable inferential properties under more general designs $X$ while having substantial power, as demonstrated in a series of experiments running on both simulated and real data.

We consider adapting a canonical computer model calibration apparatus, involving coupled Gaussian process (GP) emulators, to a computer experiment simulating radiative shock hydrodynamics that is orders of magnitude larger than what can typically be accommodated. The conventional approach calls for thousands of large matrix inverses to evaluate the likelihood in an MCMC scheme. Our approach replaces that costly ideal with a thrifty take on essential ingredients, synergizing three modern ideas in emulation, calibration and optimization: local approximate GP regression, modularization, and mesh adaptive direct search. The new methodology is motivated both by necessity—considering our particular application—and by recent trends in the supercomputer simulation literature. A synthetic data application allows us to explore the merits of several variations in a controlled environment and, together with results on our motivating real-data experiment, lead to noteworthy insights into the dynamics of radiative shocks as well as the limitations of the calibration enterprise generally.

A fundamental aspect of relational data, such as from a social network, is the possibility of dependence among the relations. In particular, the relations between members of one pair of nodes may have an effect on the relations between members of another pair. This article develops a type of regression model to estimate such effects in the context of longitudinal and multivariate relational data, or other data that can be represented in the form of a tensor. The model is based on a general multilinear tensor regression model, a special case of which is a tensor autoregression model in which the tensor of relations at one time point are parsimoniously regressed on relations from previous time points. This is done via a separable, or Kronecker-structured, regression parameter along with a separable covariance model. In the context of an analysis of longitudinal multivariate relational data, it is shown how the multilinear tensor regression model can represent patterns that often appear in relational and network data, such as reciprocity and transitivity.

Holocene (the last 12,000 years) temperature variation, including the transition out of the last Ice Age to a warmer climate, is reconstructed at multiple locations in southern Finland, Sweden and Estonia based on pollen fossil data from lake sediment cores. A novel Bayesian statistical approach is proposed that allows the reconstructed temperature histories to interact through shared environmental response parameters and spatial dependence. The prior distribution for past temperatures is partially based on numerical climate simulation. The features in the reconstructions are consistent with the quantitative climate reconstructions based on more commonly used reconstruction techniques. The results suggest that the novel spatio-temporal approach can provide quantitative reconstructions that are smoother, less uncertain and generally more realistic than the site-specific individual reconstructions.

Cardiometabolic diseases have substantially increased in China in the past 20 years and blood pressure is a primary modifiable risk factor. Using data from the China Health and Nutrition Survey, we examine blood pressure trends in China from 1991 to 2009, with a concentration on age cohorts and urbanicity. Very large values of blood pressure are of interest, so we model the conditional quantile functions of systolic and diastolic blood pressure. This allows the covariate effects in the middle of the distribution to vary from those in the upper tail, the focal point of our analysis. We join the distributions of systolic and diastolic blood pressure using a copula, which permits the relationships between the covariates and the two responses to share information and enables probabilistic statements about systolic and diastolic blood pressure jointly. Our copula maintains the marginal distributions of the group quantile effects while accounting for within-subject dependence, enabling inference at the population and subject levels. Our population-level regression effects change across quantile level, year and blood pressure type, providing a rich environment for inference. To our knowledge, this is the first quantile function model to explicitly model within-subject autocorrelation and is the first quantile function approach that simultaneously models multivariate conditional response. We find that the association between high blood pressure and living in an urban area has evolved from positive to negative, with the strongest changes occurring in the upper tail. The increase in urbanization over the last twenty years coupled with the transition from the positive association between urbanization and blood pressure in earlier years to a more uniform association with urbanization suggests increasing blood pressure over time throughout China, even in less urbanized areas. Our methods are available in the R package BSquare.

We develop methods for estimating the size of hard-to-reach populations from data collected using network-based questions on standard surveys. Such data arise by asking respondents how many people they know in a specific group (e.g., people named Michael, intravenous drug users). The Network Scale up Method (NSUM) is a tool for producing population size estimates using these indirect measures of respondents’ networks. Killworth et al. [Soc. Netw. 20 (1998a) 23–50, Evaluation Review 22 (1998b) 289–308] proposed maximum likelihood estimators of population size for a fixed effects model in which respondents’ degrees or personal network sizes are treated as fixed. We extend this by treating personal network sizes as random effects, yielding principled statements of uncertainty. This allows us to generalize the model to account for variation in people’s propensity to know people in particular subgroups (barrier effects), such as their tendency to know people like themselves, as well as their lack of awareness of or reluctance to acknowledge their contacts’ group memberships (transmission bias). NSUM estimates also suffer from recall bias, in which respondents tend to underestimate the number of members of larger groups that they know, and conversely for smaller groups. We propose a data-driven adjustment method to deal with this. Our methods perform well in simulation studies, generating improved estimates and calibrated uncertainty intervals, as well as in back estimates of real sample data. We apply them to data from a study of HIV/AIDS prevalence in Curitiba, Brazil. Our results show that when transmission bias is present, external information about its likely extent can greatly improve the estimates. The methods are implemented in the NSUM R package.

Statistical analysis and inference on spike trains is one of the central topics in the neural coding. It is of great interest to understand the underlying structure of given neural data. Based on the metric distances between spike trains, recent investigations have introduced the notion of an average or prototype spike train to characterize the template pattern in the neural activity. However, as those metrics lack certain Euclidean properties, the defined averages are nonunique, and do not share the conventional properties of a mean. In this article, we propose a new framework to define the mean spike train where we adopt a Euclidean-like metric from an $L^{p}$ family. We demonstrate that this new mean spike train properly represents the average pattern in the conventional fashion, and can be effectively computed using a theoretically-proven convergent procedure. We compare this mean with other spike train averages and demonstrate its superiority. Furthermore, we apply the new framework in a recording from rodent geniculate ganglion, where background firing activity is a common issue for neural coding. We show that the proposed mean spike train can be utilized to remove the background noise and improve decoding performance.

Dust storms in the earth’s major desert regions significantly influence microphysical weather processes, the CO$_{2}$-cycle and the global climate in general. Recent increases in the spatio-temporal resolution of remote sensing instruments have created new opportunities to understand these phenomena. However, the scale of the data collected and the inherent stochasticity of the underlying process pose significant challenges, requiring a careful combination of image processing and statistical techniques. Using satellite imagery data, we develop a statistical model of atmospheric transport that relies on a latent Gaussian Markov random field (GMRF) for inference. In doing so, we make a link between the optical flow method of Horn and Schunck and the formulation of the transport process as a latent field in a generalized linear model. We critically extend this framework to satisfy the integrated continuity equation, thereby incorporating a flow field with nonzero divergence, and show that such an approach dramatically improves performance while remaining computationally feasible. Effects such as air compressibility and satellite column projection hence become intrinsic parts of this model. We conclude with a study of the dynamics of dust storms formed over Saharan Africa and show that our methodology is able to accurately and coherently track storm movement, a critical problem in this field.

Probabilistic forecasts of wind speed are important for a wide range of applications, ranging from operational decision making in connection with wind power generation to storm warnings, ship routing and aviation. We present a statistical method that provides locally calibrated, probabilistic wind speed forecasts at any desired place within the forecast domain based on the output of a numerical weather prediction (NWP) model. Three approaches for wind speed post-processing are proposed, which use either truncated normal, gamma or truncated logistic distributions to make probabilistic predictions about future observations conditional on the forecasts of an ensemble prediction system (EPS). In order to provide probabilistic forecasts on a grid, predictive distributions that were calibrated with local wind speed observations need to be interpolated. We study several interpolation schemes that combine geostatistical methods with local information on annual mean wind speeds, and evaluate the proposed methodology with surface wind speed forecasts over Germany from the COSMO-DE (Consortium for Small-scale Modelling) ensemble prediction system.

We aim to produce predictive models that are not only accurate, but are also interpretable to human experts. Our models are decision lists, which consist of a series of if…then…statements (e.g., if high blood pressure, then stroke) that discretize a high-dimensional, multivariate feature space into a series of simple, readily interpretable decision statements. We introduce a generative model called Bayesian Rule Lists that yields a posterior distribution over possible decision lists. It employs a novel prior structure to encourage sparsity. Our experiments show that Bayesian Rule Lists has predictive accuracy on par with the current top algorithms for prediction in machine learning. Our method is motivated by recent developments in personalized medicine, and can be used to produce highly accurate and interpretable medical scoring systems. We demonstrate this by producing an alternative to the CHADS$_{2}$ score, actively used in clinical practice for estimating the risk of stroke in patients that have atrial fibrillation. Our model is as interpretable as CHADS$_{2}$, but more accurate.

Ecologists often interpret variation in the spatial distribution of populations in terms of responses to environmental features, but disentangling the effects of individual variables can be difficult if latent effects and spatial and temporal correlations are not accounted for properly. Here, we use hierarchical models based on a Poisson-lognormal mixture to understand the spatial variation in relative abundance (counts per standardized unit of effort) of yellow perch, Perca flavescens, the most abundant fish species in Lake Saint Pierre, Quebec, Canada. The mixture incorporates spatially varying environmental covariates that represent local habitat characteristics, and random temporal and spatial effects that capture the effects of unobserved ecological processes. The sampling design covers the margins but not the central region of the lake. We fit spatial generalized linear mixed models based on three different prior covariance structures for the local latent effects: a single Gaussian process (GP) over the lake, a GP over a circle, and independent GP for each shore. The models allow for independence, isotropy, or nonstationary spatial effects. Nonstationarity is dealt with using two different approaches, geometric anisotropy and the inclusion of covariates in the correlation structure of the latent spatial process. The proposed approaches for specification of spatial domain and choice of Gaussian process priors may prove useful in other applications that involve spatial correlation along an irregular contour or in discontinuous spatial domains.

This article introduces a model-based approach to distributed computing for multinomial logistic (softmax) regression. We treat counts for each response category as independent Poisson regressions via plug-in estimates for fixed effects shared across categories. The work is driven by the high-dimensional-response multinomial models that are used in analysis of a large number of random counts. Our motivating applications are in text analysis, where documents are tokenized and the token counts are modeled as arising from a multinomial dependent upon document attributes. We estimate such models for a publicly available data set of reviews from Yelp, with text regressed onto a large set of explanatory variables (user, business, and rating information). The fitted models serve as a basis for exploring the connection between words and variables of interest, for reducing dimension into supervised factor scores, and for prediction. We argue that the approach herein provides an attractive option for social scientists and other text analysts who wish to bring familiar regression tools to bear on text data.

Shifting from traditional hazard-based food safety management toward risk-based management requires statistical methods for evaluating intermediate targets in food production, such as microbiological criteria (MC), in terms of their effects on human risk of illness. A fully risk-based evaluation of MC involves several uncertainties that are related to both the underlying Quantitative Microbiological Risk Assessment (QMRA) model and the production-specific sample data on the prevalence and concentrations of microbes in production batches. We used Bayesian modeling for statistical inference and evidence synthesis of two sample data sets. Thus, parameter uncertainty was represented by a joint posterior distribution, which we then used to predict the risk and to evaluate the criteria for acceptance of production batches. We also applied the Bayesian model to compare alternative criteria, accounting for the statistical uncertainty of parameters, conditional on the data sets. Comparison of the posterior mean relative risk, $E(\mathit{RR}|\mathrm{data})=E(P(\mathrm{illness}|\mathrm{criterion\ is\ met})/P(\mathrm{illness})|\mathrm{data})$, and relative posterior risk, $\mathit{RPR}=P(\mathrm{illness}|\mathrm{data,\ criterion\ is\ met})/P(\mathrm{illness}|\mathrm{data})$, showed very similar results, but computing is more efficient for RPR. Based on the sample data, together with the QMRA model, one could achieve a relative risk of 0.4 by insisting that the default criterion be fulfilled for acceptance of each batch.

For assessing in real time the short-term trend of major economic indicators, official statistical agencies generally rely on asymmetric filters that were developed by Musgrave in 1964. However, the use of the latter introduces revisions as new observations are added to the series and, from a policy-making viewpoint, they are too slow in detecting true turning points. In this paper, we use a reproducing kernel methodology to derive asymmetric filters that converge quickly and monotonically to the corresponding symmetric one. We show theoretically that proposed criteria for time-varying bandwidth selection produce real-time trend-cycle filters to be preferred to the Musgrave filters from the viewpoint of revisions and time delay to detect true turning points. We use a set of leading, coincident and lagging indicators of the US economy to illustrate the potential gains statistical agencies could have by also using our methods in their practice.

The $p$-beauty contest is a multi-player number guessing game that is widely used to study strategic behavior. Using new data from a specially-designed web experiment, we examine the evidence in favor of a popular class of behavioral economic models called $k$-step thinking models. After fitting a custom Bayesian spline model to the experimental data, we estimate that the proportion of players who could be using a $k$-step thinking strategy is approximately 25%.

Ratings of teachers’ instructional practices using standardized classroom observation instruments are increasingly being used for both research and teacher accountability. There are multiple instruments in use, each attempting to evaluate many dimensions of teaching and classroom activities, and little is known about what underlying teaching quality attributes are being measured. We use data from multiple instruments collected from 458 middle school mathematics and English language arts teachers to inform research and practice on teacher performance measurement by modeling latent constructs of high-quality teaching. We make inferences about these constructs using a novel approach to Bayesian exploratory factor analysis (EFA) that, unlike commonly used approaches for identifying factor loadings in Bayesian EFA, is invariant to how the data dimensions are ordered. Applying this approach to ratings of lessons reveals two distinct teaching constructs in both mathematics and English language arts: (1) quality of instructional practices; and (2) quality of teacher management of classrooms. We demonstrate the relationships of these constructs to other indicators of teaching quality, including teacher content knowledge and student performance on standardized tests.

Suppose that a graph is realized from a stochastic block model where one of the blocks is of interest, but many or all of the vertices’ block labels are unobserved. The task is to order the vertices with unobserved block labels into a “nomination list” such that, with high probability, vertices from the interesting block are concentrated near the list’s beginning. We propose several vertex nomination schemes. Our basic—but principled—setting and development yields a best nomination scheme (which is a Bayes–Optimal analogue), and also a likelihood maximization nomination scheme that is practical to implement when there are a thousand vertices, and which is empirically near-optimal when the number of vertices is small enough to allow comparison to the best nomination scheme. We then illustrate the robustness of the likelihood maximization nomination scheme to the modeling challenges inherent in real data, using examples which include a social network involving human trafficking, the Enron Graph, a worm brain connectome and a political blog network.

A major challenge for cancer pathologists is to determine whether a new tumor in a patient with cancer is a metastasis or an independent occurrence of the disease. In recent years numerous studies have evaluated pairs of tumor specimens to examine the similarity of the somatic characteristics of the tumors and to test for clonal relatedness. As the landscape of mutation testing has evolved, a number of statistical methods for determining clonality have developed, notably for comparing losses of heterozygosity at candidate markers, and for comparing copy number profiles. Increasingly tumors are being evaluated for point mutations in panels of candidate genes using gene sequencing technologies. Comparison of the mutational profiles of pairs of tumors presents unusual methodological challenges: mutations at some loci are much more common than others; knowledge of the marginal mutation probabilities is scanty for most loci at which mutations might occur; the sample space of potential mutational profiles is vast. We examine this problem and propose a test for clonal relatedness of a pair of tumors from a single patient. Using simulations, its properties are shown to be promising. The method is illustrated using several examples from the literature.

Conventional approaches to statistical inference preclude structures that facilitate incorporation of supplemental information acquired from similar circumstances. For example, the analysis of data obtained using perfusion computed tomography to characterize functional imaging biomarkers in cancerous regions of the liver can benefit from partially informative data collected concurrently in noncancerous regions. This paper presents a hierarchical model structure that leverages all available information about a curve, using penalized splines, while accommodating important between-source features. Our proposed methods flexibly borrow strength from the supplemental data to a degree that reflects the commensurability of the supplemental curve with the primary curve. We investigate our method’s properties for nonparametric regression via simulation, and apply it to a set of liver cancer data. We also apply our method for a semiparametric hazard model to data from a clinical trial that compares time to disease progression for three colorectal cancer treatments, while supplementing inference with information from a previous trial that tested the current standard of care.

While studies show that autism is highly heritable, the nature of the genetic basis of this disorder remains illusive. Based on the idea that highly correlated genes are functionally interrelated and more likely to affect risk, we develop a novel statistical tool to find more potentially autism risk genes by combining the genetic association scores with gene co-expression in specific brain regions and periods of development. The gene dependence network is estimated using a novel partial neighborhood selection (PNS) algorithm, where node specific properties are incorporated into network estimation for improved statistical and computational efficiency. Then we adopt a hidden Markov random field (HMRF) model to combine the estimated network and the genetic association scores in a systematic manner. The proposed modeling framework can be naturally extended to incorporate additional structural information concerning the dependence between genes. Using currently available genetic association data from whole exome sequencing studies and brain gene expression levels, the proposed algorithm successfully identified 333 genes that plausibly affect autism risk.

Ambulatory cardiovascular (CV) measurements provide valuable insights into individuals’ health conditions in “real-life,” everyday settings. Current methods of modeling ambulatory CV data do not consider the dynamic characteristics of the full data set and their relationships with covariates such as caffeine use and stress. We propose a stochastic differential equation (SDE) in the form of a dual nonlinear Ornstein–Uhlenbeck (OU) model with person-specific covariates to capture the morning surge and nighttime dipping dynamics of ambulatory CV data. To circumvent the data analytic constraint that empirical measurements are typically collected at irregular and much larger time intervals than those evaluated in simulation studies of SDEs, we adopt a Bayesian approach with a regularized Brownian Bridge sampler (RBBS) and an efficient multiresolution (MR) algorithm to fit the proposed SDE. The MR algorithm can produce more efficient MCMC samples that is crucial for valid parameter estimation and inference. Using this model and algorithm to data from the Duke Behavioral Investigation of Hypertension Study, results indicate that age, caffeine intake, gender and race have effects on distinct dynamic characteristics of the participants’ CV trajectories.

An important endpoint variable in a cocaine rehabilitation study is the time to first relapse of a patient after the treatment. We propose a joint modeling approach based on functional data analysis to study the relationship between the baseline longitudinal cocaine-use pattern and the interval censored time to first relapse. For the baseline cocaine-use pattern, we consider both self-reported cocaine-use amount trajectories and dichotomized use trajectories. Variations within the generalized longitudinal trajectories are modeled through a latent Gaussian process, which is characterized by a few leading functional principal components. The association between the baseline longitudinal trajectories and the time to first relapse is built upon the latent principal component scores. The mean and the eigenfunctions of the latent Gaussian process as well as the hazard function of time to first relapse are modeled nonparametrically using penalized splines, and the parameters in the joint model are estimated by a Monte Carlo EM algorithm based on Metropolis–Hastings steps. An Akaike information criterion (AIC) based on effective degrees of freedom is proposed to choose the tuning parameters, and a modified empirical information is proposed to estimate the variance–covariance matrix of the estimators.

We propose and develop a Bayesian plaid model for biclustering that accounts for the prior dependency between genes (and/or conditions) through a stochastic relational graph. This work is motivated by the need for improved understanding of the molecular mechanisms of human diseases for which effective drugs are lacking, and based on the extensive raw data available through gene expression profiling. We model the prior dependency information from biological knowledge gathered from gene ontologies. Our model, the Gibbs-plaid model, assumes that the relational graph is governed by a Gibbs random field. To estimate the posterior distribution of the bicluster membership labels, we develop a stochastic algorithm that is partly based on the Wang–Landau flat-histogram algorithm. We apply our method to a gene expression database created from the study of retinal detachment, with the aim of confirming known or finding novel subnetworks of proteins associated with this disorder.

We consider the high energy physics unfolding problem where the goal is to estimate the spectrum of elementary particles given observations distorted by the limited resolution of a particle detector. This important statistical inverse problem arising in data analysis at the Large Hadron Collider at CERN consists in estimating the intensity function of an indirectly observed Poisson point process. Unfolding typically proceeds in two steps: one first produces a regularized point estimate of the unknown intensity and then uses the variability of this estimator to form frequentist confidence intervals that quantify the uncertainty of the solution. In this paper, we propose forming the point estimate using empirical Bayes estimation which enables a data-driven choice of the regularization strength through marginal maximum likelihood estimation. Observing that neither Bayesian credible intervals nor standard bootstrap confidence intervals succeed in achieving good frequentist coverage in this problem due to the inherent bias of the regularized point estimate, we introduce an iteratively bias-corrected bootstrap technique for constructing improved confidence intervals. We show using simulations that this enables us to achieve nearly nominal frequentist coverage with only a modest increase in interval length. The proposed methodology is applied to unfolding the $Z$ boson invariant mass spectrum as measured in the CMS experiment at the Large Hadron Collider.