Semi-Implicit Approaches for Large-Scale Bayesian Spatial Interpolation

Molecular ecology regularly requires the analysis of count data that reflect the relative abundance of features of a composition (e.g., taxa in a community, gene transcripts in a tissue). The sampling process that generates these data can be modelled using the multinomial distribution. Replicate multinomial samples inform the relative abundances of features in an underlying Dirichlet distribution. These distributions together form a hierarchical model for relative abundances among replicates and sampling groups. This type of Dirichlet-multinomial modelling (DMM) has been described previously, but its benefits and limitations are largely untested. With simulated data, we quantified the ability of DMM to detect differences in proportions between treatment and control groups, and compared the efficacy of three computational methods to implement DMM-Hamiltonian Monte Carlo (HMC), variational inference (VI), and Gibbs Markov chain Monte Carlo. We report that DMM was better able to detect shifts in relative abundances than analogous analytical tools, while identifying an acceptably low number of false positives. Among methods for implementing DMM, HMC provided the most accurate estimates of relative abundances, and VI was the most computationally efficient. The sensitivity of DMM was exemplified through analysis of previously published data describing lung microbiomes. We report that DMM identified several potentially pathogenic, bacterial taxa as more abundant in the lungs of children who aspirated foreign material during swallowing; these differences went undetected with different statistical approaches. Our results suggest that DMM has strong potential as a statistical method to guide inference in molecular ecology.

MotivationThe negative binomial distribution has been shown to be a good model for counts data from both bulk and single-cell RNA-sequencing (RNA-seq). Gaussian process (GP) regression provides a useful non-parametric approach for modelling temporal or spatial changes in gene expression. However, currently available GP regression methods that implement negative binomial likelihood models do not scale to the increasingly large datasets being produced by single-cell and spatial transcriptomics.ResultsThe GPcounts package implements GP regression methods for modelling counts data using a negative binomial likelihood function. Computational efficiency is achieved through the use of variational Bayesian inference. The GP function models changes in the mean of the negative binomial likelihood through a logarithmic link function and the dispersion parameter is fitted by maximum likelihood. We validate the method on simulated time course data, showing better performance to identify changes in over-dispersed counts data than methods based on Gaussian or Poisson likelihoods. To demonstrate temporal inference, we apply GPcounts to single-cell RNA-seq datasets after pseudotime and branching inference. To demonstrate spatial inference, we apply GPcounts to data from the mouse olfactory bulb to identify spatially variable genes and compare to two published GP methods. We also provide the option of modelling additional dropout using a zero-inflated negative binomial. Our results show that GPcounts can be used to model temporal and spatial counts data in cases where simpler Gaussian and Poisson likelihoods are unrealistic.Availability and implementationGPcounts is implemented using the GPflow library in Python and is available at https://github.com/ManchesterBioinference/GPcounts along with the data, code and notebooks required to reproduce the results presented here. The version used for this paper is archived at https://doi.org/10.5281/zenodo.5027066.Supplementary informationSupplementary data are available at Bioinformatics online.

Probabilistic programming languages (PPLs) are at the interface between statistics and the theory of programming languages. PPLs formulate statistical models as stochastic programs that enable automatic inference algorithms and optimization. Pyro [1] and its sibling NumPyro [2] are PPLs built on top of the deep learning frameworks PyTorch [3] and Jax [4], respectively. Both PPLs provide simple, highly similar interfaces for inference using efficient implementations of Hamiltonian Monte Carlo (HMC), the No-U-Turn Sampler (NUTS), and Stochastic Variational Inference (SVI). They automatically generate variational distributions from a model, automatically enumerate discrete variables, and support formulating deep probabilistic models such as variational autoencoders and deep Markov models. The Sine von Mises distribution and its skewed variant are toroidal distributions relevant to protein bioinformatics. They provide a natural way to model the dihedral angles of protein structures, which is important in protein structure prediction, simulation and analysis. We present efficient implementations of the Sine von Mises distribution and its skewing in Pyro and NumPyro, and devise a simulation method that increases efficiency with several orders of magnitude when using parallel hardware (i.e., modern CPUs, GPUs, and TPUs). We demonstrate the use of the skewed Sine von Mises distribution by modeling dihedral angles of proteins using a Bayesian mixture model inferred using NUTS, exploiting NumPyro's facilities for automatic enumeration [5].

Gaussian Processes (GPs) are widely-used tools in spatial statistics and machine learning and the formulae for the mean function and covariance kernel of a GP [Formula: see text] that is the image of another GP [Formula: see text] under a linear transformation [Formula: see text] acting on the sample paths of [Formula: see text] are well known, almost to the point of being folklore. However, these formulae are often used without rigorous attention to technical details, particularly when [Formula: see text] is an unbounded operator such as a differential operator, which is common in many modern applications. This note provides a self-contained proof of the claimed formulae for the case of a closed, densely-defined operator [Formula: see text] acting on the sample paths of a square-integrable (not necessarily Gaussian) stochastic process. Our proof technique relies upon Hille’s theorem for the Bochner integral of a Banach-valued random variable.

Many techniques have been proposed to model space-varying observation processes with a nonstationary spatial covariance structure and/or anisotropy, usually on a geostatistical framework. Nevertheless, there is an increasing interest in point process applications, and methodologies that take nonstationarity into account are welcomed. In this sense, this work proposes an extension of a class of spatial Cox process using spatial deformation. The proposed method enables the deformation behavior to be data-driven, through a multivariate latent Gaussian process. Inference leads to intractable posterior distributions that are approximated via MCMC. The convergence of algorithms based on the Metropolis-Hastings steps proved to be slow, and the computational efficiency of the Bayesian updating scheme was improved by adopting Hamiltonian Monte Carlo (HMC) methods. Our proposal was also compared against an alternative anisotropic formulation. Studies based on synthetic data provided empirical evidence of the benefit brought by the adoption of nonstationarity through our anisotropic structure. A real data application was conducted on the spatial spread of the Spodoptera frugiperda pest in a corn-producing agricultural area in southern Brazil. Once again, the proposed method demonstrated its benefit over alternatives.

Variational inference has become one of the most widely used methods in latent variable modeling. In its basic form, variational inference employs a fully factorized variational distribution and minimizes its Kullback–Leibler divergence to the posterior. As the minimization can only be carried out approximately, this approximation induces a bias. In this paper, we revisit perturbation theory as a powerful way of improving the variational approximation. Perturbation theory relies on a form of Taylor expansion of the log marginal likelihood, vaguely in terms of the log ratio of the true posterior and its variational approximation. While first order terms give the classical variational bound, higher-order terms yield corrections that tighten it. However, traditional perturbation theory does not provide a lower bound, making it inapt for stochastic optimization. In this paper, we derive a similar yet alternative way of deriving corrections to the evidence lower bound that resemble perturbation theory, but that result in a valid bound. We show in experiments on Gaussian processes and variational autoencoders that the new bounds are more mass covering, and that the resulting posterior covariances are closer to the true posterior and lead to higher likelihoods on held-out data.

Gaussian processes (GPs) are a popular class of Bayesian nonparametric models, but its training can be computationally burdensome for massive training datasets. While there has been notable work on scaling up these models for big data, existing methods typically rely on a stationary GP assumption for approximation, and can thus perform poorly when the underlying response surface is nonstationary, that is, it has some regions of rapid change and other regions with little change. Such non-stationarity is, however, ubiquitous in real-world problems, including our motivating application for surrogate modeling of computer experiments. We propose a new Product of Sparse GP (ProSpar-GP) method for scalable GP modeling with massive nonstationary data. The ProSpar-GP makes use of a carefully-constructed product-of-experts formulation of sparse GP experts, where different experts are placed within local regions of non-stationarity. These GP experts are fit via a novel variational inference approach, which capitalizes on mini-batching and GPU acceleration for efficient optimization of inducing points and length-scale parameters for each expert. We further show that the ProSpar-GP is Kolmogorov-consistent, in that its generative distribution defines a valid stochastic process over the prediction space; such a property provides essential stability for variational inference, particularly in the presence of non-stationarity. We then demonstrate the improved performance of the ProSpar-GP over the state-of-the-art, in a suite of numerical experiments and an application for surrogate modeling of a satellite drag simulator. Supplemental materials for this article are available online.

It is in the present decade that quanntitative-geographical studies began to undergo a marked development based on a new viewpoint in geography. In geography, too, mathematical and statistical techniques have been utilized since a long while ago as well as in other sciences. Those techniques of late are however characterized by the quantitative-geographical viewpoint as to the ‘nomothetic’ subjects in pursuit of some general law concerning the spatial regularities of various phenomena on the earth. In this respect, it should be clearly distinct from the ‘ideographic’ viewpoint in regional geography, which is interested in unique items to put emphasis on the intensive studies of individual cases, if the mathematical and statistical techniques are employed as a mere means for description of regions.The studies from quantitative-geographical viewpoint pay further attention to spatial aspects of the phenomena on the earth, make models derived from some hypothesis, put them to practical tests, and formulate the theories on which the hypothesis is based. The models used there range from deterministic to probabilistic, and static to dynamic models.The quantitative geography is assumed to contain the following seven branches of research: (1) Point Pattern Analysis, (2) Network Analysis, (3) Trend-Surface Analysis, (4) Regionalization Analysis, (5) Spatial Interaction Analysis, (6) Spatial Diffusion Analysis, and (7) Spatial Behavior Analysis. A large number of studies have added to the literatures of quantitative geography. Among the abovelisted seven branches of quantitative geography, the spatial diffusion analysis seems to show more sufficient development in connection with a theoretical formulation.The spatial diffusion analysis has been studied in a considerable number of geographical papers, since a Swedish geographer Torsten Hägerstrand understood the diffusion of information and innovation as a spatial process. The significance of Hagerstrand's studies seem to be briefly summarized in two sentences as follows: (1) the spatial process of diffusion in a geographic space was approached through his analysis of the diffusion of information and innovation, and (2) the Monte Carlo simulation technique was first employed for the analysis of geographic space. His diffusion theory of information and innovation might be regarded as comparable to the Christaller's central-place theory.The present writer attempted an analysis of the intra-urban migration in the pursuit of a theoretical formulation of geographic space, assuming the geographic space as a system of purely spatial components. And the geographic space is assumed to organize a ‘geographic field’ as probability surface given in form of ‘Mean Migration Field’ for the intra-urdan migration. Thus the MMF as a geographic field organized rationally for a given geographic space may probabilistically control the intra-urban migration as a spatial and geographic phenomenon. An example was drawn from an analysis of intra-urban migation on the built-up area, defined as a Densely Inhabited District by the national population census, of Fukaya city (approximately 60, 000 population), Saitama prefecture. The following four stochastic models were formulated to test the hypotheses concerning geographic space considered for the intra-urban migration by means of the Hägerstrand type of Monte Carlo simulation: (1) MODEL I to simulate the migration under a complete random process on a hypothesis of the completely uniform geographic space, ‘Isotropic Plain’, (2) MODEL II for a simulation of the migration pattern within a geographic space assuming the friction effect of distance involved, (3) MODEL III by which simulation is attempted for the intra-urban migration within a geographic space characterized by the uneven distribution of population or the densities effect

Uncertainty propagation using infinite mixture of Gaussian processes and variational Bayesian inference

We investigate how well a variational Bayes (VB) method performs as an approximation for the demanding computations typically required in Gaussian process regression. We study a technique using so-called emph{inducing variables} to drastically reduce the computational load. The key question is: can we make equally good estimates with less computational effort? We show that this is indeed possible, provided that one carefully chooses emph{how many} and emph{which} inducing variables. If too few are chosen, the resulting model becomes too `smooth', leading to inaccurate estimates and an overestimation of uncertainty. However, with a sufficient number of well-chosen inducing variables, the VB method can be just as accurate as the original, much slower approach. Faster computation with retained accuracy is feasible, as long as the number of inducing variables is tailored to the complexity of the true underlying function. Besides providing a point estimate, we also want to know how certain we can be about it. A key advantage of the Bayesian approach is that it allows for emph{uncertainty quantification}. We analyse VB methods based on spectral features and examine how well these methods capture uncertainty. Surprisingly, these methods can yield reliable uncertainty estimates under certain conditions, which is not always true for other VB approaches. Thus, this method can offer not only accurate predictions but also trustworthy uncertainty estimates, provided the priors are well-matched to the underlying data. Finally, we address a practical challenge: how do you choose the right settings (hyperparameters) for a model when you don’t know how complex the real world is? We consider an emph{adaptive} VB method that adjusts itself to the complexity of the data without requiring prior knowledge. We combine a hierarchical Bayesian model with a variational approach, and demonstrate both theoretically and through experiments that this works well. The result is a fast method that automatically adapts to various types of underlying functions while still achieving optimal performance. Self-tuning VB methods are therefore possible and deliver strong results without requiring prior knowledge of the complexity of the data-generating mechanism.

Inversion of hierarchical Bayesian models using Gaussian processes

In many machine learning problems, one has to work with data of different types, including continuous, discrete, and categorical data. Further, it is often the case that many of these data are missing from the database. This paper proposes a Gaussian process framework that efficiently captures the information from mixed numerical and categorical data that effectively incorporates missing variables. First, we propose a generative model for the mixed-type data. The generative model exploits Gaussian processes with kernels constructed from the latent vectors. We also propose a method for inference of the unknowns, and in its implementation, we rely on a sparse spectrum approximation of the Gaussian processes and variational inference. We demonstrate the performance of the method for both supervised and unsupervised tasks. First, we investigate the imputation of missing variables in an unsupervised setting, and then we show the results of joint imputation and classification on IBM employee data.

The article develops a hybrid variational Bayes (VB) algorithm that combines the mean-field and stochastic linear regression fixed-form VB methods. The new estimation algorithm can be used to approximate any posterior without relying on conjugate priors. We propose a divide and recombine strategy for the analysis of large datasets, which partitions a large dataset into smaller subsets and then combines the variational distributions that have been learned in parallel on each separate subset using the hybrid VB algorithm. We also describe an efficient model selection strategy using cross-validation, which is straightforward to implement as a by-product of the parallel run. The proposed method is applied to fitting generalized linear mixed models. The computational efficiency of the parallel and hybrid VB algorithm is demonstrated on several simulated and real datasets. Supplementary material for this article is available online.

Statistical models for spatial processes play a central role in analyses of spatial data. Yet, it is the simple, interpretable, and well understood models that are routinely employed even though, as is revealed through prior and posterior predictive checks, these can poorly characterise the spatial heterogeneity in the underlying process of interest. Here, we propose a new, flexible class of spatial-process models, which we refer to as spatial Bayesian neural networks (SBNNs). An SBNN leverages the representational capacity of a Bayesian neural network; it is tailored to a spatial setting by incorporating a spatial “embedding layer” into the network and, possibly, spatially-varying network parameters. An SBNN is calibrated by matching its finite-dimensional distribution at locations on a fine gridding of space to that of a target process of interest. That process could be easy to simulate from or we may have many realisations from it. We propose several variants of SBNNs, most of which are able to match the finite-dimensional distribution of the target process at the selected grid better than conventional BNNs of similar complexity. We also show that an SBNN can be used to represent a variety of spatial processes often used in practice, such as Gaussian processes, lognormal processes, and max-stable processes. We briefly discuss the tools that could be used to make inference with SBNNs, and we conclude with a discussion of their advantages and limitations.

The Expected Value of Perfect Partial Information (EVPPI) is a decision‐theoretic measure of the ‘cost’ of parametric uncertainty in decision making used principally in health economic decision making. Despite this decision‐theoretic grounding, the uptake of EVPPI calculations in practice has been slow. This is in part due to the prohibitive computational time required to estimate the EVPPI via Monte Carlo simulations. However, recent developments have demonstrated that the EVPPI can be estimated by non‐parametric regression methods, which have significantly decreased the computation time required to approximate the EVPPI. Under certain circumstances, high‐dimensional Gaussian Process (GP) regression is suggested, but this can still be prohibitively expensive. Applying fast computation methods developed in spatial statistics using Integrated Nested Laplace Approximations (INLA) and projecting from a high‐dimensional into a low‐dimensional input space allows us to decrease the computation time for fitting these high‐dimensional GP, often substantially. We demonstrate that the EVPPI calculated using our method for GP regression is in line with the standard GP regression method and that despite the apparent methodological complexity of this new method, R functions are available in the package BCEA to implement it simply and efficiently. © 2016 The Authors. Statistics in Medicine Published by John Wiley & Sons Ltd.