Daily precipitation has an enormous impact on human activity, and the study of how it varies over time and space, and what global indicators influence it, is of paramount importance to Australian agriculture. We analyze over 294 million daily rainfall measurements since 1876, spanning 17,606 sites across continental Australia. The data are not only large but also complex, and the topic would benefit from a common and publicly available statistical framework. We propose a Bayesian hierarchical mixture model that accommodates mixed discrete-continuous data. The observational level describes site-specific temporal and climatic variation via a mixture-of-experts model. At the next level of the hierarchy, spatial variability of the mixture weights’ parameters is modeled by a spatial Gaussian process prior. A parallel and distributed Markov chain Monte Carlo sampler is developed which scales the model to large data sets. We present examples of posterior inference on the mixture weights, monthly intensity levels, daily temporal dependence, offsite prediction of the effects of climate drivers and long-term rainfall trends across the entire continent. Computer code implementing the methods proposed in this paper is available as an R package.

In the last 25 years there has been an important increase in the amount of data collected from animal-mounted sensors (bio-probes) which are often used to study the animals’ behaviour or environment. We focus here on an example of the latter, where the interest is in sea surface temperature (SST), and measurements are taken from sensors mounted on elephant seals in the southern Indian Ocean. We show that standard geostatistical models may not be reliable for this type of data, due to the possibility that the regions visited by the animals may depend on the SST. This phenomenon is know in the literature as preferential sampling, and, if ignored, it may affect the resulting spatial predictions and parameter estimates. Research on this topic has been mostly restricted to stationary sampling locations such as monitoring sites. The main contribution of this manuscript is to extend this methodology to observations obtained by devices that move through the region of interest, as is the case with the tagged seals. More specifically, we propose a flexible framework for inference on preferentially sampled fields where the process that generates the sampling locations is stochastic and moving over time through a two-dimensional space. Our simulation studies confirm that predictions obtained from the preferential sampling model are more reliable when this phenomenon is present, and they compare very well to the standard ones when there is no preferential sampling. Finally, we note that the conclusions of our analysis of the SST data can change considerably when we incorporate preferential sampling in the model.

Computer simulation models are central to environmental science. These mathematical models are used to understand complex weather and climate patterns and to predict the climate’s response to different forcings. Climate models are of course not perfect reflections of reality, and so comparison with observed data is needed to quantify and to correct for biases and other deficiencies. We propose a new method to calibrate model output using observed data. Our approach not only matches the marginal distributions of the model output and gridded observed data, but it simultaneously postprocesses the model output to have the same spatial correlation as the observed data. This comprehensive calibration method permits realistic spatial simulations for regional impact studies. We apply the proposed method to global climate model output in North America and show that it successfully calibrates the model output for temperature and precipitation.

The conventional model for assessing insensitivity to hidden bias in paired observational studies constructs a worst-case distribution for treatment assignments subject to bounds on the maximal bias to which any given pair is subjected. In studies where rare cases of extreme hidden bias are suspected, the maximal bias may be substantially larger than the typical bias across pairs, such that a correctly specified bound on the maximal bias would yield an unduly pessimistic perception of the study’s robustness to hidden bias. We present an extended sensitivity analysis which allows researchers to simultaneously bound the maximal and typical bias perturbing the pairs under investigation while maintaining the desired Type I error rate. We motivate and illustrate our method with two sibling studies on the impact of schooling on earnings, one containing information of cognitive ability of siblings and the other not. Cognitive ability, clearly influential of both earnings and degree of schooling, is likely similar between members of most sibling pairs yet could, conceivably, vary drastically for some siblings. The method is straightforward to implement, simply requiring the solution to a quadratic program. R code is provided in the Supplementary Material.

We introduce a Bayesian multivariate hierarchical framework to estimate a space-time model for a joint series of monthly extreme temperatures and amounts of precipitation. Data are available for 360 monitoring stations over 60 years, with missing data affecting almost all series. Model components account for spatio-temporal correlation and annual cycles, dependence on covariates and between responses. Spatio-temporal dependence is modeled by the nearest neighbor Gaussian process (GP), response multivariate dependencies are represented by the linear model of coregionalization and effects of annual cycles are included by a circular representation of time. The proposed approach allows imputation of missing values and interpolation of climate surfaces at the national level. It also provides a characterization of the so called Italian ecoregions, namely broad and discrete ecologically homogeneous areas of similar potential as regards the climate, physiography, hydrography, vegetation and wildlife. To now, Italian ecoregions are hierarchically classified into 4 tiers that go from 2 Divisions to 35 Subsections and are defined by informed expert judgments. The current climatic characterization of Italian ecoregions is based on bioclimatic indices for the period 1955–2000.

With the growing cost of health care in the United States, the need to improve efficiency and efficacy has become increasingly urgent. There has been a keen interest in developing interventions to effectively coordinate the typically fragmented care of patients with many comorbidities. Evaluation of such interventions is often challenging given their long-term nature and their differential effectiveness among different patients. Furthermore, care coordination interventions are often highly resource-intensive. Hence there is pressing need to identify which patients would benefit the most from a care coordination program. In this work we introduce a subgroup identification procedure for long-term interventions whose effects are expected to change smoothly over time. We allow differential effects of an intervention to vary over time and encourage these effects to be more similar for closer time points by utilizing a fused lasso penalty. Our approach allows for flexible modeling of temporally changing intervention effects while also borrowing strength in estimation over time. We utilize our approach to construct a personalized enrollment decision rule for a complex case management intervention in a large health system and demonstrate that the enrollment decision rule results in improvement in health outcomes and care costs. The proposed methodology could have broad usage for the analysis of different types of long-term interventions or treatments including other interventions commonly implemented in health systems.

Bulk gene expression experiments relied on aggregations of thousands of cells to measure the average expression in an organism. Advances in microfluidic and droplet sequencing now permit expression profiling in single cells. This study of cell-to-cell variation reveals that individual cells lack detectable expression of transcripts that appear abundant on a population level, giving rise to zero-inflated expression patterns. To infer gene coregulatory networks from such data, we propose a multivariate Hurdle model. It is comprised of a mixture of singular Gaussian distributions. We employ neighborhood selection with the pseudo-likelihood and a group lasso penalty to select and fit undirected graphical models that capture conditional independences between genes. The proposed method is more sensitive than existing approaches in simulations, even under departures from our Hurdle model. The method is applied to data for T follicular helper cells, and a high-dimensional profile of mouse dendritic cells. It infers network structure not revealed by other methods, or in bulk data sets. A R implementation is available at https://github.com/amcdavid/HurdleNormal.

We present TreeClone, a latent feature allocation model to reconstruct tumor subclones subject to phylogenetic evolution that mimics tumor evolution. Similar to most current methods, we consider data from next-generation sequencing of tumor DNA. Unlike most methods that use information in short reads mapped to single nucleotide variants (SNVs), we consider subclone phylogeny reconstruction using pairs of two proximal SNVs that can be mapped by the same short reads. As part of the Bayesian inference model, we construct a phylogenetic tree prior. The use of the tree structure in the prior greatly strengthens inference. Only subclones that can be explained by a phylogenetic tree are assigned non-negligible probabilities. The proposed Bayesian framework implies posterior distributions on the number of subclones, their genotypes, cellular proportions and the phylogenetic tree spanned by the inferred subclones. The proposed method is validated against different sets of simulated and real-world data using single and multiple tumor samples. An open source software package is available at http://www.compgenome.org/treeclone.

The Eurovision Song Contest is a popular TV singing competition held annually among country members of the European Broadcasting Union. In this competition, each member can be both contestant and jury, as it can participate with a song and/or vote for other countries’ tunes. During the years, the voting system has repeatedly been accused of being biased by tactical voting; votes would represent strategic interests rather than actual musical preferences of the voting countries. In this work, we develop a latent space model to investigate the presence of a latent structure underlying the exchange of votes. Focusing on the period from 1998 to 2015, we represent the vote exchange as a multivariate network: each edition is a network, where countries are the nodes and two countries are linked by an edge if one voted for the other. The different networks are taken to be independent replicates of a conditional Bernoulli distribution, with success probability specified as a function of a common latent space capturing the overall relationships among the countries. Proximity denotes similarity, and countries close in the latent space are more likely to exchange votes. If the exchange of votes depends on the similarity between countries, the quality of the competing songs might not be a relevant factor in the determination of the voting preferences, and this would suggest the presence of some bias. A Bayesian hierarchical modelling approach is employed to estimate the parameters, where the probability of a connection between any two countries is a function of their distance in the latent space, network-specific parameters and edge-specific covariates. The estimated latent space is found to be relevant in the determination of edge probabilities, however, the positions of the countries in such space only partially correspond to their actual geographical positions.

There is a growing awareness of the important roles that microbial communities play in complex biological processes. Modern investigation of these often uses next generation sequencing of metagenomic samples to determine community composition. We propose a statistical technique based on clique loglinear models and Bayes model averaging to identify microbial components in a metagenomic sample at various taxonomic levels that have significant associations. We describe the model class, a stochastic search technique for model selection, and the calculation of estimates of posterior probabilities of interest. We demonstrate our approach using data from the Human Microbiome Project and from a study of the skin microbiome in chronic wound healing. Our technique also identifies significant dependencies among microbial components as evidence of possible microbial syntrophy.

The central aim in this paper is to address variable selection questions in nonlinear and nonparametric regression. Motivated by statistical genetics, where nonlinear interactions are of particular interest, we introduce a novel and interpretable way to summarize the relative importance of predictor variables. Methodologically, we develop the “RelATive cEntrality” (RATE) measure to prioritize candidate genetic variants that are not just marginally important, but whose associations also stem from significant covarying relationships with other variants in the data. We illustrate RATE through Bayesian Gaussian process regression, but the methodological innovations apply to other “black box” methods. It is known that nonlinear models often exhibit greater predictive accuracy than linear models, particularly for phenotypes generated by complex genetic architectures. With detailed simulations and two real data association mapping studies, we show that applying RATE enables an explanation for this improved performance.

We develop the hierarchical cluster coherence (HCC) method for brain signals, a procedure for characterizing connectivity in a network by clustering nodes or groups of channels that display a high level of coordination as measured by “cluster-coherence.” While the most common approach to measure dependence between clusters is through pairs of single time series, our method proposes cluster coherence which measures dependence between pairs of whole clusters rather than between single elements. Thus it takes into account both the dependence between clusters and within channels in a cluster. The identified clusters contain time series that exhibit high cross-dependence in the spectral domain. Simulation studies demonstrate that the proposed HCC method is competitive with the other feature-based clustering methods. To study clustering in a network of multichannel electroencephalograms (EEG) during an epileptic seizure, we applied the HCC method and identified connectivity on alpha $(8,12)$ Hertz and beta $(16,30)$ Hertz bands at different phases of the recording: before an epileptic seizure, during the early and middle phases of the seizure episode. To increase the potential impact of HCC in neuroscience, we also developed the HCC-Vis, an R-Shiny app (RStudio), which can be downloaded from https://carolinaeuan.shinyapps.io/hcc-vis/.

Principal component analysis (PCA) is a commonly used statistical method in a wide range of applications. However, it does not work well when the number of features is larger than the sample size. We consider the estimation of the sparse principal subspace in the high dimensional setting with missing data motivated by the analysis of single-cell RNA sequence data. We propose a two step estimation procedure, and establish the rates of convergence for estimating the principal subspace. Simulated examples with various missing mechanisms show its competitive performance compared to existing sparse PCA methods. We apply the method to single-cell data and show that the proposed method can better distinguish cell types than other PCA methods.

Exploratory analysis is an important first step for discovering latent structure and generating hypotheses in large biological data sets. However, when the number of variables is large compared to the number of samples, standard methods such as principal components analysis give results that are unstable and difficult to interpret.

Here, we present adaptive generalized principal components analysis (adaptive gPCA), a new method that solves these problems by incorporating information about the relationships among the variables. Adaptive gPCA gives a low-dimensional representation of the samples with axes that are interpretable in terms of groups of closely related variables. We show that adaptive gPCA does well at reconstructing true latent structure in simulated data and demonstrate its use on a study of the effect of antibiotics on the human gut microbiota.

Glass fragments are often compared in the course of a forensic evaluation using their chemical composition determined with technologies such as LA-ICP-MS. At present forensic scientists advocate the use of two comparison criteria based on univariate intervals around all mean elemental concentrations for fragments originating from a known piece of broken glass. The main drawback of this approach is that it does not consider the dependencies between concentrations. Further, when the elemental concentrations are more variable within panes, it becomes harder to reject the null hypothesis of no difference between fragments. In the legal context higher variance would tend to incriminate the defendant because the intervals would tend to be wider. We demonstrate that a score-based approach to assess the probative value of evidence in glass comparisons outperforms the two standard interval methods and other methods proposed in the literature, at least in terms of minimizing classification error in the glass fragment sources we analyzed. We use machine learning algorithms to construct a similarity score between pairs of glass fragments. The learning algorithms exploit the dependencies among elemental concentrations and result in an empirical class probability; so, we can report the degree of similarity between two fragments. Our group is in the process of assembling the first glass composition database with enough information within and between glass samples to permit computing well-conditioned estimates of high-dimensional covariance matrices. These data will be available to anyone who wishes to carry out research in this area.

The advent of high-throughput sequencing technologies has led to an increasing availability of large multi-tissue data sets which contain gene expression measurements across different tissues and individuals. In this setting, variation in expression levels arises due to contributions specific to genes, tissues, individuals, and interactions thereof. Classical clustering methods are ill-suited to explore these three-way interactions and struggle to fully extract the insights into transcriptome complexity contained in the data. We propose a new statistical method, called MultiCluster, based on semi-nonnegative tensor decomposition which permits the investigation of transcriptome variation across individuals and tissues simultaneously. We further develop a tensor projection procedure which detects covariate-related genes with high power, demonstrating the advantage of tensor-based methods in incorporating information across similar tissues. Through simulation and application to the GTEx RNA-seq data from 53 human tissues, we show that MultiCluster identifies three-way interactions with high accuracy and robustness.

This work is motivated by the problem of testing for differences in the mean electricity prices before and after Germany’s abrupt nuclear phaseout after the nuclear disaster in Fukushima Daiichi, Japan, in mid-March 2011. Taking into account the nature of the data and the auction design of the electricity market, we approach this problem using a Local Linear Kernel (LLK) estimator for the nonparametric mean function of sparse covariate-adjusted functional data. We build upon recent theoretical work on the LLK estimator and propose a two-sample test statistics using a finite sample correction to avoid size distortions. Our nonparametric test results on the price differences point to a Simpson’s paradox explaining an unexpected result recently reported in the literature.

An important objective in vaccine studies entails identifying an immune response which is predictive of disease risk. Nonparametric methods are developed for inference on immune response thresholds that are associated with specified levels of disease risk, including where the risk level is zero. This threshold is defined as the minimum immune response value above which disease risk is less than or equal to the desired level. The proposed nonparametric methods are compared to previously developed parametric methods in simulation studies. The methods are extended for use in studies that only measure the immune response in a subset of participants, such as case-cohort or case-control studies, and with right censored time to disease outcomes. Finally, these methods are used to estimate neutralizing antibody thresholds for virologically confirmed dengue risk using data from two recent dengue vaccine trials.

The Italian National Institute for Statistics regularly provides estimates of unemployment indicators using data from the labor force survey. However, direct estimates of unemployment incidence cannot be released for local labor market areas. These are unplanned domains defined as clusters of municipalities; many are out-of-sample areas, and the majority is characterized by a small sample size which renders direct estimates inadequate. The empirical best predictor represents an appropriate, model-based alternative. However, for non-Gaussian responses its computation and the computation of the analytic approximation to its mean squared error require the solution of (possibly) multiple integrals that, generally, have not a closed form. To solve the issue, Monte Carlo methods and parametric bootstrap are common choices, even though the computational burden is a nontrivial task. In this paper, we propose a semiparametric empirical best predictor for a (possibly) nonlinear mixed effect model by leaving the distribution of the area-specific random effects unspecified and estimating it from the observed data. This approach is known to lead to a discrete mixing distribution which helps avoid unverifiable parametric assumptions and heavy integral approximations. We also derive a second-order, bias-corrected analytic approximation to the corresponding mean squared error. Finite sample properties of the proposed approach are tested via a large scale simulation study. Furthermore, the proposal is applied to unit-level data from the 2012 Italian Labor Force Survey to estimate unemployment incidence for 611 local labor market areas using auxiliary information from administrative registers and the 2011 Census.

For zero-truncated count data, as they typically arise in capture-recapture modelling, we consider modelling under one-inflation. This is motivated by police data on drink-driving in Britain which shows high one-inflation. The data, which are used here, are from the years 2011 to 2015 and are based on DR10 endorsements. We show that inference for an arbitrary count density with one-inflation can be equivalently based upon the associated zero-one truncated count density. This simplifies inference considerably including maximum likelihood estimation and likelihood ratio testing. For the drink-driving application, we use the geometric distribution which shows a good fit. We estimate the total drink-driving as about $2{,}300{,}000$ drink-drivers in the observational period. As $227{,}578$ were observed, this means that only about 10% of the drink-driving population is observed with a bootstrap confidence interval of 9%–12%.

Tumor cells acquire different genetic alterations during the course of evolution in cancer patients. As a result of competition and selection, only a few subgroups of cells with distinct genotypes survive. These subgroups of cells are often referred to as subclones. In recent years, many statistical and computational methods have been developed to identify tumor subclones, leading to biologically significant discoveries and shedding light on tumor progression, metastasis, drug resistance and other processes. However, most existing methods are either not able to infer the phylogenetic structure among subclones, or not able to incorporate copy number variations (CNV). In this article, we propose SIFA (tumor Subclone Identification by Feature Allocation), a Bayesian model which takes into account both CNV and tumor phylogeny structure to infer tumor subclones. We compare the performance of SIFA with two other commonly used methods using simulation studies with varying sequencing depth, evolutionary tree size, and tree complexity. SIFA consistently yields better results in terms of Rand Index and cellularity estimation accuracy. The usefulness of SIFA is also demonstrated through its application to whole genome sequencing (WGS) samples from four patients in a breast cancer study.

Most causal inference studies rely on the assumption of overlap to estimate population or sample average causal effects. When data suffer from non-overlap, estimation of these estimands requires reliance on model specifications due to poor data support. All existing methods to address non-overlap, such as trimming or down-weighting data in regions of poor data support, change the estimand so that inference cannot be made on the sample or the underlying population. In environmental health research settings where study results are often intended to influence policy, population-level inference may be critical and changes in the estimand can diminish the impact of the study results, because estimates may not be representative of effects in the population of interest to policymakers. Researchers may be willing to make additional, minimal modeling assumptions in order to preserve the ability to estimate population average causal effects. We seek to make two contributions on this topic. First, we propose a flexible, data-driven definition of propensity score overlap and non-overlap regions. Second, we develop a novel Bayesian framework to estimate population average causal effects with minor model dependence and appropriately large uncertainties in the presence of non-overlap and causal effect heterogeneity. In this approach the tasks of estimating causal effects in the overlap and non-overlap regions are delegated to two distinct models suited to the degree of data support in each region. Tree ensembles are used to nonparametrically estimate individual causal effects in the overlap region, where the data can speak for themselves. In the non-overlap region where insufficient data support means reliance on model specification is necessary, individual causal effects are estimated by extrapolating trends from the overlap region via a spline model. The promising performance of our method is demonstrated in simulations. Finally, we utilize our method to perform a novel investigation of the causal effect of natural gas compressor station exposure on cancer outcomes. Code and data to implement the method and reproduce all simulations and analyses is available on Github (https://github.com/rachelnethery/overlap).

Antibodies, an essential part of our immune system, develop through an intricate process to bind a wide array of pathogens. This process involves randomly mutating DNA sequences encoding these antibodies to find variants with improved binding, though mutations are not distributed uniformly across sequence sites. Immunologists observe this nonuniformity to be consistent with “mutation motifs” which are short DNA subsequences that affect how likely a given site is to experience a mutation. Quantifying the effect of motifs on mutation rates is challenging. A large number of possible motifs makes this statistical problem high dimensional, while the unobserved history of the mutation process leads to a nontrivial missing data problem. We introduce an $\ell_{1}$-penalized proportional hazards model to infer mutation motifs and their effects. In order to estimate model parameters, our method uses a Monte Carlo EM algorithm to marginalize over the unknown ordering of mutations. We show that our method performs better on simulated data compared to current methods and leads to more parsimonious models. The application of proportional hazards to mutation processes is, to our knowledge, novel and formalizes the current methods in a statistical framework that can be easily extended to analyze the effect of other biological features on mutation rates.

Neurological diseases are due to the loss of structure or function of neurons that eventually leads to cognitive deficit, neuropsychiatric symptoms, and impaired activities of daily living. Identifying sensitive and specific biological and clinical markers for early diagnosis allows recruiting patients into a clinical trial to test therapeutic intervention. However, many biomarker studies considered a single biomarker at one time that fails to provide precise prediction for disease age at onset. In this paper, we use longitudinally collected measurements from multiple biomarkers and measurement error-corrected clinical diagnosis ages to identify which biomarkers and what features of biomarker trajectories are useful for early diagnosis. Specifically, we assume that the subject-specific biomarker trajectories depend on unobserved states of underlying latent variables with the conditional mean follows a nonlinear sigmoid shape. We show that peak degeneration age of the biomarker trajectory is useful for early diagnosis. We propose an Expectation-Maximization (EM) algorithm to obtain the maximum likelihood estimates of all parameters and conduct extensive simulation studies to examine the performance of the proposed methods. Finally, we apply our methods to studies of Alzheimer’s disease and Huntington’s disease and identify a few important biomarkers that can be used for early diagnosis.