As the first successful technique used to detect exoplanets orbiting distant stars, the radial velocity method aims to detect a periodic Doppler shift in a stellar spectrum due to the star’s motion along the line sight. We introduce a new, mathematically rigorous approach to detect such a signal that accounts for the smooth functional relationship of neighboring wavelengths in the spectrum, minimizes the role of wavelength interpolation, accounts for heteroskedastic noise and easily allows for accurate calculation of the estimated radial velocity standard error. Using Hermite–Gaussian functions, we show that the problem of detecting a Doppler shift in the spectrum can be reduced to linear regression in many settings. A simulation study demonstrates that the proposed method is able to accurately estimate an individual spectrum’s radial velocity with precision below $0.3{\text{m s}}^{-1}$, corresponding to a Doppler shift much smaller than the size of a spectral pixel. Furthermore, the new method outperforms the traditional cross-correlation function approach for estimating the radial velocity by reducing the root mean squared error up to $15\text{cm}{\text{s}}^{-1}$. The proposed method is also demonstrated on a new set of observations from the EXtreme PREcision Spectrometer (EXPRES) for the host star 51 Pegasi, and successfully recovers estimates of the planetary companion’s parameters that agree well with previous studies. The method is implemented in the R package rvmethod, and supplemental Python code is also available.

Snow density estimates below the surface, used with airplane-acquired ice-penetrating radar measurements, give a site-specific history of snow water accumulation. Because it is infeasible to drill snow cores across all of Antarctica to measure snow density and because it is critical to understand how climatic changes are affecting the world’s largest freshwater reservoir, we develop methods that enable snow density estimation with uncertainty in regions where snow cores have not been drilled.

In inland West Antarctica, snow density increases monotonically as a function of depth, except for possible microscale variability or measurement error, and it cannot exceed the density of ice. We present a novel class of integrated spatial process models that allow interpolation of monotone snow density curves. For computational feasibility we construct the space-depth process through kernel convolutions of log-Gaussian spatial processes. We discuss model comparison, model fitting and prediction. Using this model, we extend estimates of snow density beyond the depth of the original core and estimate snow density curves where snow cores have not been drilled. Along flight lines with ice-penetrating radar, we use interpolated snow density curves to estimate recent water accumulation and find predominantly decreasing water accumulation over recent decades.

In this work, we estimate extreme sea surface temperature (SST) hotspots, that is, high threshold exceedance regions, for the Red Sea, a vital region of high biodiversity. We analyze high-resolution satellite-derived SST data comprising daily measurements at 16,703 grid cells across the Red Sea over the period 1985–2015. We propose a semiparametric Bayesian spatial mixed-effects linear model with a flexible mean structure to capture spatially-varying trend and seasonality, while the residual spatial variability is modeled through a Dirichlet process mixture (DPM) of low-rank spatial Student’s t processes (LTPs). By specifying cluster-specific parameters for each LTP mixture component, the bulk of the SST residuals influence tail inference and hotspot estimation only moderately. Our proposed model has a nonstationary mean, covariance, and tail dependence, and posterior inference can be drawn efficiently through Gibbs sampling. In our application, we show that the proposed method outperforms some natural parametric and semiparametric alternatives. Moreover, we show how hotspots can be identified, and we estimate extreme SST hotspots for the whole Red Sea, projected until the year 2100, based on the Representative Concentration Pathways 4.5 and 8.5. The estimated 95% credible region, for joint high threshold exceedances include large areas covering major endangered coral reefs in the southern Red Sea.

Improving and optimizing oceanographic sampling is a crucial task for marine science and maritime resource management. Faced with limited resources in understanding processes in the water column, the combination of statistics and autonomous systems provides new opportunities for experimental design. In this work we develop efficient spatial sampling methods for characterizing regions, defined by simultaneous exceedances above prescribed thresholds of several responses, with an application focus on mapping coastal ocean phenomena based on temperature and salinity measurements. Specifically, we define a design criterion based on uncertainty in the excursions of vector-valued Gaussian random fields and derive tractable expressions for the expected integrated Bernoulli variance reduction in such a framework. We demonstrate how this criterion can be used to prioritize sampling efforts at locations that are ambiguous, making exploration more effective. We use simulations to study and compare properties of the considered approaches, followed by results from field deployments with an autonomous underwater vehicle as part of a study mapping the boundary of a river plume. The results demonstrate the potential of combining statistical methods and robotic platforms to effectively inform and execute data-driven environmental sampling.

The classification of shapes is of great interest in diverse areas ranging from medical imaging to computer vision and beyond. While many statistical frameworks have been developed for the classification problem, most are strongly tied to early formulations of the problem with an object to be classified described as a vector in a relatively low-dimensional Euclidean space. Statistical shape data have two main properties that suggest a need for a novel approach: (i) shapes are inherently infinite-dimensional with strong dependence among the positions of nearby points, and (ii) shape space is not Euclidean but is fundamentally curved. To accommodate these features of the data, we work with the square-root velocity function of the curves to provide a useful formal description of the shape, pass to tangent spaces of the manifold of shapes at projection points (which effectively separate shapes for pairwise classification in the training data) and use principal components within these tangent spaces to reduce dimensionality. We illustrate the impact of the projection point and choice of subspace on the misclassification rate with a novel method of combining pairwise classifiers.

The recent curation of large-scale databases with 3D surface scans of shapes has motivated the development of tools that better detect global patterns in morphological variation. Studies, which focus on identifying differences between shapes, have been limited to simple pairwise comparisons and rely on prespecified landmarks (that are often known). We present SINATRA, the first statistical pipeline for analyzing collections of shapes without requiring any correspondences. Our novel algorithm takes in two classes of shapes and highlights the physical features that best describe the variation between them. We use a rigorous simulation framework to assess our approach. Lastly, as a case study we use SINATRA to analyze mandibular molars from four different suborders of primates and demonstrate its ability recover known morphometric variation across phylogenies.

The period-luminosity relation (PLR) of Mira variable stars is an important tool to determine astronomical distances. The common approach of estimating the PLR is a two-step procedure that first estimates the Mira periods and then runs a linear regression of magnitude on log period. When the light curves are sparse and noisy, the accuracy of period estimation decreases and can suffer from aliasing effects. Some methods improve accuracy by incorporating complex model structures at the expense of significant computational costs. Another drawback of existing methods is that they only provide point estimation without proper estimation of uncertainty. To overcome these challenges, we develop a hierarchical Bayesian model that simultaneously models the quasi-periodic variations for a collection of Mira light curves while estimating their common PLR. By borrowing strengths through the PLR, our method automatically reduces the aliasing effect, improves the accuracy of period estimation and is capable of characterizing the estimation uncertainty. We develop a scalable stochastic variational inference algorithm for computation that can effectively deal with the multimodal posterior of period. The effectiveness of the proposed method is demonstrated through simulations and an application to observations of Miras in the Local Group galaxy M33. Without using ad hoc period correction tricks, our method achieves a distance estimate of M33 that is consistent with published work. Our method also shows superior robustness to downsampling of the light curves.

Nitrogen dioxide (${\mathrm{NO}}_2$) is a primary constituent of traffic-related air pollution and has well-established harmful environmental and human-health impacts. Knowledge of the spatiotemporal distribution of ${\mathrm{NO}}_2$ is critical for exposure and risk assessment. A common approach for assessing air pollution exposure is linear regression involving spatially referenced covariates, known as land-use regression (LUR). We develop a scalable approach for simultaneous variable selection and estimation of LUR models with spatiotemporally correlated errors, by combining a general-Vecchia Gaussian-process approximation with a penalty on the LUR coefficients. In comparison to existing methods using simulated data, our approach resulted in higher model-selection specificity and sensitivity and in better prediction in terms of calibration and sharpness, for a wide range of relevant settings. In our spatiotemporal analysis of daily, US-wide, ground-level ${\mathrm{NO}}_2$ data, our approach was more accurate, and produced a sparser and more interpretable model. Our daily predictions elucidate spatiotemporal patterns of ${\mathrm{NO}}_2$ concentrations across the United States, including significant variations between cities and intra-urban variation. Thus, our predictions will be useful for epidemiological and risk-assessment studies seeking daily, national-scale predictions, and they can be used in acute-outcome health-risk assessments.

Sea ice, or frozen ocean water, freezes and melts every year in the Arctic. Forecasts of where sea ice will be located weeks to months in advance have become more important as the amount of sea ice declines due to climate change, for maritime planning and other uses. Typical sea ice forecasts are made with ensemble models, physics-based models of sea ice and the surrounding ocean and atmosphere. This paper introduces Mixture Contour Forecasting, a method to forecast sea ice probabilistically using a mixture of two distributions, one based on postprocessed output from ensembles and the other on observed sea ice patterns in recent years. At short lead times, these forecasts are better calibrated than unadjusted dynamic ensemble forecasts and other statistical reference forecasts. To produce these forecasts, a statistical technique is introduced that directly models the sea ice edge contour, the boundary around the region that is ice-covered. Mixture Contour Forecasting and reference methods are evaluated for monthly sea ice forecasts for 2008–2016 at lead times ranging from 0.5–6.5 months using one of the European Centre for Medium-Range Weather Forecasts ensembles.

Future grid management systems will coordinate distributed production and storage resources to manage, in a cost effective fashion, the increased load and variability brought by the electrification of transportation and by a higher share of weather dependent production. Electricity demand forecasts at a low level of aggregation will be key inputs for such systems. We focus on forecasting demand at the individual household level, which is more challenging than forecasting aggregate demand, due to the lower signal-to-noise ratio and to the heterogeneity of consumption patterns across households. We propose a new ensemble method for probabilistic forecasting which borrows strength across the households while accommodating their individual idiosyncrasies. In particular, we develop a set of models or “experts” which capture different demand dynamics, and we fit each of them to the data from each household. Then, we construct an aggregation of experts where the ensemble weights are estimated on the whole data set, the main innovation being that we let the weights vary with the covariates by adopting an additive model structure. In particular, the proposed aggregation method is an extension of regression stacking where the mixture weights are modelled using linear combinations of parametric, smooth or random effects. The methods for building and fitting additive stacking models are implemented by the gamFactory R package, available at https://github.com/mfasiolo/gamFactory.

Causal mediation analysis seeks to investigate how the treatment effect of an exposure on outcomes is mediated through intermediate variables. Although many applications involve longitudinal data, the existing methods are not directly applicable to settings where the mediator and outcome are measured on sparse and irregular time grids. We extend the existing causal mediation framework from a functional data analysis perspective, viewing the sparse and irregular longitudinal data as realizations of underlying smooth stochastic processes. We define causal estimands of direct and indirect effects accordingly and provide corresponding identification assumptions. For estimation and inference, we employ a functional principal component analysis approach for dimension reduction and use the first few functional principal components instead of the whole trajectories in the structural equation models. We adopt the Bayesian paradigm to accurately quantify the uncertainties. The operating characteristics of the proposed methods are examined via simulations. We apply the proposed methods to a longitudinal data set from a wild baboon population in Kenya to investigate the causal relationships between early adversity, strength of social bonds between animals and adult glucocorticoid hormone concentrations. We find that early adversity has a significant direct effect (a 9–14% increase) on females’ glucocorticoid concentrations across adulthood but find little evidence that these effects were mediated by weak social bonds.

Composed of a large number of subwavelength unit cells with designable geometries, metamaterials have been widely studied to achieve extraordinary advantageous and unusual optical properties. However, ordinary computer simulator requires a time-consuming fine-tuning to find a proper design of metamaterial for a specific optical property, making the design stage a critical bottleneck in large scale applications of metamaterials. This paper investigates the metamaterial design under the framework of computer experiments, with emphasis on dealing with the challenge of designing numerous unit cells with functional responses, simultaneously, which is not common in traditional computer experiments. We formulate the multiple related design targets as a multitarget design problem. Leveraging on the similarity between different designs, we propose an efficient Bayesian optimization strategy with a parsimonious surrogate model and an integrated acquisition function to design multiple unit cells with very few function evaluations. A wide range of simulations confirm the effectiveness and superiority of the proposed approach compared to the naive strategies where the multiple unit cells are dealt with separately or sequentially. Such a rapid design strategy has the potential to greatly promote large scale applications of metamaterials in practice.

We consider an extension to discrete-space, continuous-time models for animal movement that have previously been presented in the literature. The extension from a continuous-time Markov formulation to a continuous-time semi-Markov formulation allows for the inclusion of temporally dynamic habitat conditions as well as temporally changing movement responses by animals to that environment. We show that, with only a little additional consideration, the Poisson likelihood calculation for the Markov version can still be used within the multiple imputation framework commonly employed for analysis of telemetry data. In addition, we consider a Bayesian model selection methodology within the imputation framework. The model selection method uses a Laplace approximation to the posterior model probability to provide a computationally feasible approach. The full methodology is then used to analyze movements of 15 weaned northern fur seal (Callorhinus ursinus) pups with respect to surface winds, geostrophic currents and sea surface temperature. The highest posterior model probabilities belonged to those models containing only winds and current; SST was not a significant factor for modeling their movement.

We propose a Bayesian hierarchical Jolly–Seber model that can accommodate a semiparametric functional relationship between external covariates and capture probabilities, individual heterogeneity in departure due to an internal time-varying covariate and the dependence of arrival time on external covariates. Modelwise, we consider a stochastic process to characterize the evolution of the partially observable internal covariate that is linked to departure probabilities. Computationally, we develop a well-tailored Markov chain Monte Carlo algorithm that is free of tuning through data augmentation. Inferentially, our model allows us to make inference about stopover duration and population sizes, the impacts of various covariates on departure and arrival time and to identify flexible yet data-driven functional relationships between external covariates and capture probabilities. We demonstrate the effectiveness of our model through a motivating dataset collected for studying the migration of mallards (Anas platyrhynchos) in Sweden.

Recent advances in next-generation sequencing technology have yielded huge amounts of transcriptomic data. The discreteness and the high dimensions of RNA-seq data have posed great challenges in biological network analysis. Although estimation theories for high-dimensional modified Poisson-type graphical models have been proposed for the network analysis of count-valued data, the statistical inference of these models is still largely unknown. We herein propose a two-step procedure in both edgewise and global statistical inference of these modified Poisson-type graphical models using a cutting-edge generalized low-dimensional projection approach for bias correction. Extensive simulations and a real example with ground truth illustrate asymptotic normality of edgewise inference and more accurate inferential results in multiple testing compared to the sole estimation and the inferential method under normal assumption. Furthermore, the application of our method to novel RNA-seq data of childhood atopic asthma in Puerto Ricans demonstrates more biologically meaningful results compared to the sole estimation and the inferential methods based on Gaussian and nonparanormal graphical models.

We propose a novel way of modelling time-varying networks by inducing two-way sparsity on local models of node connectivity. This two-way sparsity separately promotes sparsity across time and sparsity across variables (within time). Separation of these two types of sparsity is achieved through a novel prior structure which draws on ideas from the Bayesian lasso and from copula modelling. We provide an efficient implementation of the proposed model via a Gibbs sampler, and we apply the model to data from neural development. In doing so, we demonstrate that the proposed model is able to identify changes in genomic network structure that match current biological knowledge. Such changes in genomic network structure can then be used by neurobiologists to identify potential targets for further experimental investigation.

On the problem of scoring genes for evidence of changes in the distribution of single-cell expression, we introduce an empirical Bayesian mixture approach and evaluate its operating characteristics in a range of numerical experiments. The proposed approach leverages cell-subtype structure revealed in cluster analysis in order to boost gene-level information on expression changes. Cell clustering informs gene-level analysis through a specially-constructed prior distribution over pairs of multinomial probability vectors; this prior meshes with available model-based tools that score patterns of differential expression over multiple subtypes. We derive an explicit formula for the posterior probability that a gene has the same distribution in two cellular conditions, allowing for a gene-specific mixture over subtypes in each condition. Advantage is gained by the compositional structure of the model not only in which a host of gene-specific mixture components are allowed but also in which the mixing proportions are constrained at the whole cell level. This structure leads to a novel form of information sharing through which the cell-clustering results support gene-level scoring of differential distribution. The result, according to our numerical experiments, is improved sensitivity compared to several standard approaches for detecting distributional expression changes.

Measuring the dependence of $\mathit{k}\ge 3$ random variables and drawing inference from such higher-order dependences are scientifically important yet challenging. Motivated here by protein coevolution with multivariate categorical features, we consider an information theoretic measure of higher-order dependence. The proposed collective dependence is a symmetrization of differential interaction information which generalizes the mutual information of a pair of random variables. We show that the collective dependence can be easily estimated and facilitates a test on the dependence of $\mathit{k}\ge 3$ random variables. Upon carefully exploring the null space of collective dependence, we devise a Classification-Assisted Large scaLe inference procedure to DEtect significant k-COllective DEpendence among $\mathit{d}\ge \mathit{k}$ random variables, with the false discovery rate controlled. Finite sample performance of our method is examined via simulations. We apply this method to the multiple protein sequence alignment data to study the residue or position coevolution for two protein families, the elongation factor P family and the zinc knuckle family. We identify novel functional triplets of amino acid residues, whose contributions to the protein function are further investigated. These confirm that the collective dependence does yield additional information important for understanding the protein coevolution compared to the pairwise measures.

There are distinguishing features or “hallmarks” of cancer that are found across tumors, individuals and types of cancer, and these hallmarks can be driven by specific genetic mutations. Yet within a single tumor there is often extensive genetic heterogeneity as evidenced by single-cell and bulk DNA sequencing data. The goal of this work is to jointly infer the underlying genotypes of tumor subpopulations and the distribution of those subpopulations in individual tumors by integrating single-cell and bulk sequencing data. Understanding the genetic composition of the tumor at the time of treatment is important in the personalized design of targeted therapeutic combinations and monitoring for possible recurrence after treatment.

We propose a hierarchical Dirichlet process mixture model that incorporates the correlation structure induced by a structured sampling arrangement, and we show that this model improves the quality of inference. We develop a representation of the hierarchical Dirichlet process prior as a Gamma–Poisson hierarchy, and we use this representation to derive a fast Gibbs sampling inference algorithm using the augment-and-marginalize method. Experiments with simulation data show that our model outperforms standard numerical and statistical methods for decomposing admixed count data. Analyses of real acute lymphoblastic leukemia cancer sequencing dataset shows that our model improves upon state-of-the-art bioinformatic methods. An interpretation of the results of our model on this real dataset reveals comutated loci across samples.

In this paper an appropriate and interpretable diagnosis statistical model is proposed to predict Nonalcoholic Steatohepatitis (NASH) from near infrared spectrometry data. In this disease, unknown patients’ profiles are expected to lead to a different diagnosis. The model has then to take into account the heterogeneity of the data and the dimension of the spectrometric data.

To this end, we propose to fit a mixture model on the joint distribution of the diagnostic binary variable and the covariates selected in the spectra. The penalized maximum likelihood estimator is considered. In practice, a twofold penalty on both regression coefficients and covariance parameters is imposed. Automatic selection criteria, such as the AIC and BIC, are used to select the amount of shrinkage and the number of clusters. The performance of the overall procedure is evaluated by a simulation study, and its application on the NASH data set is analyzed. The model leads to better prediction performance than competitive methods and provides highly interpretable results.

Phylogenetic comparative methods correct for shared evolutionary history among a set of nonindependent organisms by modeling sample traits as arising from a diffusion process along the branches of a possibly unknown history. To incorporate such uncertainty, we present a scalable Bayesian inference framework under a general Gaussian trait evolution model that exploits Hamiltonian Monte Carlo (HMC). HMC enables efficient sampling of the constrained model parameters and takes advantage of the tree structure for fast likelihood and gradient computations, yielding algorithmic complexity linear in the number of observations. This approach encompasses a wide family of stochastic processes, including the general Ornstein–Uhlenbeck (OU) process, with possible missing data and measurement errors. We implement inference tools for a biologically relevant subset of all these models into the BEAST phylogenetic software package and develop model comparison through marginal likelihood estimation. We apply our approach to study the morphological evolution in the superfamily of Musteloidea (including weasels and allies) as well as the heritability of HIV virulence. This second problem furnishes a new measure of evolutionary heritability that demonstrates its utility through a targeted simulation study.

Transcriptome-wide association studies based on genetically predicted gene expression have the potential to identify novel regions associated with various complex traits. It has been shown that incorporating expression quantitative trait loci (eQTLs) corresponding to multiple tissue types can improve power for association studies involving complex etiology. In this article we propose a new multivariate response linear regression model and method for predicting gene expression in multiple tissues simultaneously. Unlike existing methods for multitissue joint eQTL mapping, our approach incorporates tissue-tissue expression correlation which allows us to more efficiently handle missing expression measurements and to more accurately predict gene expression using a weighted summation of eQTL genotypes. We show through simulation studies that our approach performs better than the existing methods in many scenarios. We use our method to estimate eQTL weights for 29 tissues collected by GTEx, and show that our approach significantly improves expression prediction accuracy compared to competitors. Using our eQTL weights, we perform a multitissue-based S-MultiXcan (PLoS Genet. 15 (2019) e1007889) transcriptome-wide association study and show that our method leads to more discoveries in novel regions and more discoveries overall than the existing methods. Estimated eQTL weights and code for implementing the method are available for download online at github.com/ajmolstad/MTeQTLResults.

There has been a long-standing interest in the analysis of university professor salary data. The vast majority of the publications on the topic employ linear regression models in an attempt to predict individual salaries. Indeed, the administration of any academic institution is interested in adequately compensating the faculty to attract and keep the best specialists available on the market. However, higher administration and legislators are not concerned with the matter of individual compensation and need to have a bigger picture for developing university strategies and policies. This paper is the first attempt to model university compensation data at the institutional level. The analysis of university salary patterns is a challenging problem due to the heterogeneous, skewed, multiway and temporal nature of the data. This paper aims at addressing all the above-mentioned issues by proposing a novel tensor regression mixture model and applying it to the data set obtained from the American Association of University Professors. The utility of the developed model is illustrated on addressing several important questions related to gender equity and peer institution comparison.

Statistical models of community-level homeless rates typically assume a linear relationship to covariates. This linear model assumption precludes the possibility of inflection points in homeless rates—thresholds in quantifiable metrics of a community that, once breached, are associated with large increases in homelessness. In this paper we identify points of structural change in the relationship between homeless rates and community-level measures of housing affordability and extreme poverty. We utilize the Ewens–Pitman attraction (EPA) distribution to develop a Bayesian nonparametric regression model in which clusters of communities with similar covariates share common patterns of variation in homeless rates. A main finding of the study is that the expected homeless rate in a community begins to quickly increase once median rental costs exceed 30% of median income, providing a statistical link between homelessness and the U.S. government’s definition of a housing cost burden. Our analysis also identifies clusters of communities that exhibit distinct geographic patterns and yields insight into the homelessness and housing affordability crisis unfolding on both coasts of the United States.

Cross-sectional length-biased data arise from questions on the at-risk time for an event of interest from those who are at risk but have yet to experience the event. For example, in the National Survey on Family Growth (NSFG) women, who were currently attempting to become pregnant, were asked how long they had been attempting pregnancy. Cross-sectional survival analysis methods use the observed at-risk times to make inference on the distribution of the unobserved time-to-failure. However, methodological gaps in these methods remain such as how to handle semicompeting risks. For example, if the women attempting pregnancy had undergone fertility treatment during their current pregnancy attempt. In this paper we develop statistical methods that extend cross-sectional survival analysis methods to incorporate semicompeting risks. They can be used to estimate the distribution of the length of natural pregnancy attempts (i.e., without fertility treatment) while correctly accounting for women that sought fertility treatment prior to being sampled using cross-sectional data. We demonstrate our approach based on simulated data and an analysis of data from the NSFG. The proposed method results in separate survival curves for time-to-natural-pregnancy, time-to-fertility treatment and time-to-pregnancy after fertility treatment.

Surveys designed to collect data on similar variables using samples representing the same population may still result in different estimates due, for example, to differences in sample designs or modes of data collection. Considered in this paper is the case where two surveys were conducted concurrently, with one using the same methodology as used in prior rounds of the survey and the other using an updated methodology, resulting in substantial differences in several key estimates. Due to differences in sample size, only the latter survey was detailed enough for disaggregated-level estimates of publishable quality. We propose a hierarchical model to account for discrepancies in the estimates from the two surveys and a Bayesian approach for producing reliable estimates at various levels of aggregation. The model relies on a common latent structure at the disaggregated level to allow “bridging” between the two surveys. The methodology is applied to the 2016 National Survey of Fishing, Hunting and Wildlife-Associated Recreation and the 2016 50-State Surveys of Fishing, Hunting and Wildlife-Related Recreation. Aligning these two surveys is critical to extend the series of related statistics that have been published since 1955, allowing for meaningful comparisons over time despite the change in survey methodology.