Latent Gaussian Process Research Articles

Iterative Methods for Vecchia-Laplace Approximations for Latent Gaussian Process Models

Latent Gaussian process (GP) models are flexible probabilistic non-parametric function models. Vecchia approximations are accurate approximations for GPs to overcome computational bottlenecks for large data, and the Laplace approximation is a fast method with asymptotic convergence guarantees to approximate marginal likelihoods and posterior predictive distributions for non-Gaussian likelihoods. Unfortunately, the computational complexity of combined Vecchia-Laplace approximations grows faster than linearly in the sample size when used in combination with direct solver methods such as the Cholesky decomposition. Computations with Vecchia-Laplace approximations can thus become prohibitively slow precisely when the approximations are usually the most accurate, i.e., on large data sets. In this article, we present iterative methods to overcome this drawback. Among other things, we introduce and analyze several preconditioners, derive new convergence results, and propose novel methods for accurately approximating predictive variances. We analyze our proposed methods theoretically and in experiments with simulated and real-world data. In particular, we obtain a speed-up of an order of magnitude compared to Cholesky-based calculations and a threefold increase in prediction accuracy in terms of the continuous ranked probability score compared to a state-of-the-art method on a large satellite data set. All methods are implemented in a free C++ software library with high-level Python and R packages. 4

Generalized Additive Models

Generalized additive models are generalized linear models in which the linear predictor includes a sum of smooth functions of covariates, where the shape of the functions is to be estimated. They have also been generalized beyond the original generalized linear model setting to distributions outside the exponential family and to situations in which multiple parameters of the response distribution may depend on sums of smooth functions of covariates. The widely used computational and inferential framework in which the smooth terms are represented as latent Gaussian processes, splines, or Gaussian random effects is reviewed, paying particular attention to the case in which computational and theoretical tractability is obtained by prior rank reduction of the model terms. An empirical Bayes approach is taken, and its relatively good frequentist performance discussed, along with some more overtly frequentist approaches to model selection. Estimation of the degree of smoothness of component functions via cross validation or marginal likelihood is covered, alongside the computational strategies required in practice, including when data and models are reasonably large. It is briefly shown how the framework extends easily to location-scale modeling, and, with more effort, to techniques such as quantile regression. Also covered are the main classes of smooths of multiple covariates that may be included in models: isotropic splines and tensor product smooth interaction terms.

Generating Independent Replicates Directly from the Posterior Distribution for a Class of Spatial Hierarchical Models

Markov chain Monte Carlo (MCMC) allows one to generate dependent replicates from a posterior distribution for effectively any Bayesian hierarchical model. However, MCMC can produce a significant computational burden. This motivates us to consider finding expressions of the posterior distribution that are computationally straightforward to obtain independent replicates from directly. We focus on a broad class of Bayesian hierarchical models for spatially dependent data, which are often modeled via a latent Gaussian process (LGP). First, we derive a new class of distributions referred to as the generalized conjugate multivariate (GCM) distribution. The GCM distribution’s theoretical development follows that of the conjugate multivariate (CM) distribution with two main differences: the GCM allows for latent Gaussian process assumptions, and the GCM explicitly accounts for hyperparameters through marginalization. The development of GCM is needed to obtain independent replicates directly from the exact posterior distribution, which has an efficient regression form. Hence, we refer to our method as Exact Posterior Regression (EPR). Simulation studies with weakly stationary spatial processes and spatial basis function expansions are provided. We provide an analysis of poverty incidence from the U.S. Census Bureau, and an analysis of high-dimensional remote sensing data. Supplementary materials for this article are available online.

Mapping socio-economic status using mixed data: a hierarchical Bayesian approach

Abstract We propose a Bayesian hierarchical model to estimate a socio-economic status (SES) index based on mixed dichotomous and continuous variables. In particular, we extend Quinn’s ([2004]. Bayesian factor analysis for mixed ordinal and continuous responses. Political Analysis, 12(4), 338–353. https://doi.org/10.1093/pan/mph022) and Schliep and Hoeting’s ([2013]. Multilevel latent Gaussian process model for mixed discrete and continuous multivariate response data. Journal of Agricultural, Biological, and Environmental Statistics, 18(4), 492–513. https://doi.org/10.1007/s13253-013-0136-z) factor analysis models for mixed dichotomous and continuous variables by allowing a spatial hierarchical structure of key parameters of the model. Unlike most SES assessment models proposed in the literature, the hierarchical nature of this model enables the use of census observations at the household level without needing to aggregate any information a priori. Therefore, it better accommodates the variability of the SES between census tracts and the number of households per area. The proposed model is used in the estimation of a socio-economic index using 10% of the 2010 Ghana census in the Greater Accra Metropolitan area. Out of the 20 observed variables, the number of people per room, access to water piping and flushable toilets differentiated high and low SES areas the best.

Bayesian Modeling for Nonstationary Spatial Point Process via Spatial Deformations.

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.

Unsupervised Anomaly Detection via Nonlinear Manifold Learning

Abstract Anomalies are samples that significantly deviate from the rest of the data and their detection plays a major role in building machine learning models that can be reliably used in applications such as data-driven design and novelty detection. The majority of existing anomaly detection methods either are exclusively developed for (semi) supervised settings, or provide poor performance in unsupervised applications where there are no training data with labeled anomalous samples. To bridge this research gap, we introduce a robust, efficient, and interpretable methodology based on nonlinear manifold learning to detect anomalies in unsupervised settings. The essence of our approach is to learn a low-dimensional and interpretable latent representation (aka manifold) for all the data points such that normal samples are automatically clustered together and hence can be easily and robustly identified. We learn this low-dimensional manifold by designing a learning algorithm that leverages either a latent map Gaussian process (LMGP) or a deep autoencoder (AE). Our LMGP-based approach, in particular, provides a probabilistic perspective on the learning task and is ideal for high-dimensional applications with scarce data. We demonstrate the superior performance of our approach over existing technologies via multiple analytic examples and real-world datasets.

Bayesian hierarchical model for bias-correcting climate models

Abstract. Climate models, derived from process understanding, are essential tools in the study of climate change and its wide-ranging impacts. Hindcast and future simulations provide comprehensive spatiotemporal estimates of climatology that are frequently employed within the environmental sciences community, although the output can be afflicted with bias that impedes direct interpretation. Post-processing bias correction approaches utilise observational data to address this challenge, although they are typically criticised for not being physically justified and not considering uncertainty in the correction. This paper proposes a novel Bayesian bias correction framework that robustly propagates uncertainty and models underlying spatial covariance patterns. Shared latent Gaussian processes are assumed between the in situ observations and climate model output, with the aim of partially preserving the covariance structure from the climate model after bias correction, which is based on well-established physical laws. Results demonstrate added value in modelling shared generating processes under several simulated scenarios, with the most value added for the case of sparse in situ observations and smooth underlying bias. Additionally, the propagation of uncertainty to a simulated final bias-corrected time series is illustrated, which is of key importance to a range of stakeholders, such as climate scientists engaged in impact studies, decision-makers trying to understand the likelihood of particular scenarios and individuals involved in climate change adaption strategies where accurate risk assessment is required for optimal resource allocation. This paper focuses on one-dimensional simulated examples for clarity, although the code implementation is developed to also work on multi-dimensional input data, encouraging follow-on real-world application studies that will further validate performance and remaining limitations. The Bayesian framework supports uncertainty propagation under model adaptations required for specific applications, providing a flexible approach that increases the scope of data assimilation tasks more generally.

A geospatial analysis of local intermediate snail host distributions provides insight into schistosomiasis risk within under-sampled areas of southern Lake Malawi

BackgroundAlong the southern shoreline of Lake Malawi, the incidence of schistosomiasis is increasing with snails of the genera Bulinus and Biomphalaria transmitting urogenital and intestinal schistosomiasis, respectively. Since the underlying distribution of snails is partially known, often being focal, developing pragmatic spatial models that interpolate snail information across under-sampled regions is required to understand and assess current and future risk of schistosomiasis.MethodsA secondary geospatial analysis of recently collected malacological and environmental survey data was undertaken. Using a Bayesian Poisson latent Gaussian process model, abundance data were fitted for Bulinus and Biomphalaria. Interpolating the abundance of snails along the shoreline (given their relative distance along the shoreline) was achieved by smoothing, using extracted environmental rainfall, land surface temperature (LST), evapotranspiration, normalised difference vegetation index (NDVI) and soil type covariate data for all predicted locations. Our adopted model used a combination of two-dimensional (2D) and one dimensional (1D) mapping.ResultsA significant association between normalised difference vegetation index (NDVI) and abundance of Bulinus spp. was detected (log risk ratio − 0.83, 95% CrI − 1.57, − 0.09). A qualitatively similar association was found between NDVI and Biomphalaria sp. but was not statistically significant (log risk ratio − 1.42, 95% CrI − 3.09, 0.10). Analyses of all other environmental data were considered non-significant.ConclusionsThe spatial range in which interpolation of snail distributions is possible appears < 10km owing to fine-scale biotic and abiotic heterogeneities. The forthcoming challenge is to refine geospatial sampling frameworks with future opportunities to map schistosomiasis within actual or predicted snail distributions. In so doing, this would better reveal local environmental transmission possibilities.Graphical

A Bayesian survival treed hazards model using latent Gaussian processes.

Survival models are used to analyze time-to-event data in a variety of disciplines. Proportional hazard models provide interpretable parameter estimates, but proportional hazard assumptions are not always appropriate. Non-parametric models are more flexible but often lack a clear inferential framework. We propose a Bayesian treed hazards partition model that is both flexible and inferential. Inference is obtained through the posterior tree structure and flexibility is preserved by modeling the log-hazard function in each partition using a latent Gaussian process. An efficient reversible jump Markov chain Monte Carlo algorithm is accomplished by marginalizing the parameters in each partition element via a Laplace approximation. Consistency properties for the estimator are established. The method can be used to help determine subgroups as well as prognostic and/or predictive biomarkers in time-to-event data. The method is compared with some existing methods on simulated data and a liver cirrhosis dataset.

A spatio-temporal model of multi-marker antimalarial resistance.

The emergence and spread of drug-resistant Plasmodium falciparum parasites have hindered efforts to eliminate malaria. Monitoring the spread of drug resistance is vital, as drug resistance can lead to widespread treatment failure. We develop a Bayesian model to produce spatio-temporal maps that depict the spread of drug resistance, and apply our methods for the antimalarial sulfadoxine-pyrimethamine. We infer from genetic count data the prevalences over space and time of various malaria parasite haplotypes associated with drug resistance. Previous work has focused on inferring the prevalence of individual molecular markers. In reality, combinations of mutations at multiple markers confer varying degrees of drug resistance to the parasite, indicating that multiple markers should be modelled together. However, the reporting of genetic count data is often inconsistent as some studies report haplotype counts, whereas some studies report mutation counts of individual markers separately. In response, we introduce a latent multinomial Gaussian process model to handle partially reported spatio-temporal count data. As drug-resistant mutations are often used as a proxy for treatment efficacy, point estimates from our spatio-temporal maps can help inform antimalarial drug policies, whereas the uncertainties from our maps can help with optimizing sampling strategies for future monitoring of drug resistance.

Latent Gaussian Processes Based Graph Learning for Urban Traffic Prediction

Latent Gaussian Processes Based Graph Learning for Urban Traffic Prediction

Mixed-Variable Global Sensitivity Analysis for Knowledge Discovery and Efficient Combinatorial Materials Design

Abstract Global Sensitivity Analysis (GSA) is the study of the influence of any given input on the outputs of a model. In the context of engineering design, GSA has been widely used to understand both individual and collective contributions of design variables on the design objectives. So far, global sensitivity studies have often been limited to design spaces with only quantitative (numerical) design variables. However, many engineering systems also contain, if not only, qualitative (categorical) design variables in addition to quantitative design variables. In this paper, we integrate Latent Variable Gaussian Process (LVGP) with Sobol’ analysis to develop the first metamodel-based mixed-variable GSA method. Through numerical case studies, we validate and demonstrate the effectiveness of our proposed method for mixed-variable problems. Furthermore, while the proposed GSA method is general enough to benefit various engineering design applications, we integrate it with multi-objective Bayesian optimization (BO) to create a sensitivity-aware design framework in accelerating the Pareto front design exploration for metal-organic framework (MOF) materials with many-level combinatorial design spaces. Although MOFs are constructed only from qualitative variables that are notoriously difficult to design, our method can utilize sensitivity analysis to navigate the optimization in the many-level large combinatorial design space, greatly expediting the exploration of novel MOF candidates.

Delineating important killer whale foraging areas using a spatiotemporal logistic model

Conservation management planning for highly mobile species requires an understanding of the distribution of areas that are biologically important to the species of concern. Collecting data on the locations of animal behaviors linked to biological characteristics, such as foraging, can be used to spatially describe biological important areas. However, spatial modeling of free-ranging animal behavior can be challenging, as behavioral observations of animals are often clustered, and sampling is commonly conducted at a higher frequency than changes in behavioral states, resulting in data that are usually highly autocorrelated in space and time. Here, we fit latent Gaussian process models to observational behavioral data to generate spatially-explicit predictions of foraging behavior within the critical habitat of an endangered population of fish-eating killer whales (Orcinus orca) in southern British Columbia, Canada, and northern Washington State, USA. We compare spatial models treating temporal autocorrelation in behavior in three ways: (1) ignoring temporal autocorrelation entirely; (2) traditional data-thinning to remove temporal autocorrelation, and; (3) using a temporal Gaussian process to account for temporal autocorrelation. Comparisons of autocorrelative structures for each model and visual comparison of broad spatial patterns demonstrate that our third approach yields more accurate results than when ignoring temporal autocorrelation entirely and higher precision results than when applying data-thinning methods. Within the identified areas of critical habitat, our models indicate two primary regions of intense killer whale foraging activity, and we delineate areas wherein the probability of foraging was particularly high as candidate locations for conservation management actions. This study underscores the value of refining our understanding of high-use areas for at-risk species by incorporating animal behavior data to inform area-based conservation measures.

Nonparametric quantile scalar-on-image regression

A quantile scalar-on-image regression model is developed to comprehensively study the relationship between cognitive decline and various clinical covariates and imaging factors. As a motivating example, the high-dimensional brain imaging data from the research on Alzheimer's disease are considered predictors of patients' cognitive decline. A Bayesian nonparametric model is proposed to handle the complex spatially distributed imaging data, where the coefficient function is assumed to be a latent Gaussian process. A soft-thresholding operator is introduced to capture the sparse structure of the regression coefficients. Utilizing kernel basis functions to approximate the latent Gaussian process facilitates easy-to-implement computation and consistent estimation. Inference is performed within the Bayesian framework, using an efficient Markov chain Monte Carlo algorithm. The proposed method is compared with the functional principal component analysis method in simulations and applied to a study of Alzheimer's disease.

Rapid design of top-performing metal-organic frameworks with qualitative representations of building blocks

Data-driven materials design often encounters challenges where systems possess qualitative (categorical) information. Specifically, representing Metal-organic frameworks (MOFs) through different building blocks poses a challenge for designers to incorporate qualitative information into design optimization, and leads to a combinatorial challenge, with large number of MOFs that could be explored. In this work, we integrated Latent Variable Gaussian Process (LVGP) and Multi-Objective Batch-Bayesian Optimization (MOBBO) to identify top-performing MOFs adaptively, autonomously, and efficiently. We showcased that our method (i) requires no specific physical descriptors and only uses building blocks that construct the MOFs for global optimization through qualitative representations, (ii) is application and property independent, and (iii) provides an interpretable model of building blocks with physical justification. By searching only ~1% of the design space, LVGP-MOBBO identified all MOFs on the Pareto front and 97% of the 50 top-performing designs for the CO2 working capacity and CO2/N2 selectivity properties.

Estimation and inference in a high-dimensional semiparametric Gaussian copula vector autoregressive model

This paper develops simple, robust estimation and inference methods for the transition matrix of a high-dimensional semiparametric Gaussian copula vector autoregressive process. Our estimator is based on rank estimators of the large variance and auto-covariance matrices of a transformed latent high-dimensional Gaussian process. We derive rates of convergence of our estimator based on which we develop de-biased inference for Granger causality. Numerical results demonstrate the efficacy of the proposed methods. Although our focus is on the observable process, by the nature of rank estimators, all the methods developed directly apply to the transformed latent process. In technical terms, our analysis relies heavily on newly developed exponential inequalities for (degenerate) U-statistics under α-mixing condition.

High-resolution global precipitation downscaling with latent Gaussian models and non-stationary stochastic partial differential equation structure

Abstract Obtaining high-resolution maps of precipitation data can provide key insights to stakeholders to assess a sustainable access to water resources at urban scale. Mapping a non-stationary, sparse process such as precipitation at very high spatial resolution requires the interpolation of global datasets at the location where ground stations are available with statistical models able to capture complex non-Gaussian global space–time dependence structures. In this work, we propose a new approach based on capturing the spatially varying anisotropy of a latent Gaussian process via a locally deformed stochastic partial differential equation (SPDE) with a buffer allowing for a different spatial structure across land and sea. The finite volume approximation of the SPDE, coupled with integrated nested Laplace approximation ensures feasible Bayesian inference for tens of millions of observations. The simulation studies showcase the improved predictability of the proposed approach against stationary and no-buffer alternatives. The proposed approach is then used to yield high-resolution simulations of daily precipitation across the United States.

Robust optimization of 3D printing process parameters considering process stability and production efficiency

With the rapid evolution and application of 3D printing technology, production capacity can be easily expanded by increasing the number of 3D printing machines or implementing distributed production to meet customer demand. According to our best knowledge, no paper has been found to address the process stability and production efficiency of 3D printers used in mass production or distributed manufacturing. This paper takes the material extrusion (MEX) 3D printing process as an example, aiming to minimize the difference in the quality of parts fabricated by different printers or nozzles at a lower printing cost by selecting the appropriate process parameter settings. First, combine the latent variable Gaussian process (LVGP) model with the split-plot experimental design to investigate the mechanical properties and costs of printing parts in relation to critical process parameters. Second, under consideration of the in-process uncertainty, the interval number theory is used to measure the split-plot effect of the mechanical properties of printing parts. Finally, a multi-objective optimization model integrating the split-plot effect and printing cost is established to get the Pareto optimal solution. The grey relational analysis (GRA) is used to select the optimal compromise solutions from the Pareto solution sets. The results of the validation experiments show that a more efficient and stable printing process can be achieved under the optimal parameter settings obtained by the proposed approach. This study will provide theoretical references and technical support for the standardization and industrialization of 3D printing technology.

A latent Gaussian process model for the spatial distribution of liquefaction manifestation

This paper presents a model for distributing zones of liquefaction and nonliquefaction for use in regional liquefaction risk analysis. There are two broad methodologies that have been used to evaluate liquefaction risk on the regional scale: (a) application of site-specific procedures using soil properties inferred from geology, or (b) application of geospatial proxies for liquefaction. The first approach will tend to predict similar liquefaction probabilities across broad areas with similar geology, water table depths, and shaking intensities. The second approach yields the probability of liquefaction, which can be interpreted as the portion of the area affected by liquefaction ([Formula: see text]). Neither approach, however, gives an informed prediction of the spatial distribution of liquefaction and the resulting displacements, which are particularly important for assessments of seismic risk for spatially distributed infrastructure systems. We propose a methodology for incorporating spatial correlation into a geospatial proxy for liquefaction to create maps of liquefaction and nonliquefaction for a given earthquake scenario. First, we describe a latent Gaussian process that is assumed to govern the spatial distribution of liquefaction. Next, a database of empirical observations of liquefaction is used to obtain the coefficients that describe that latent Gaussian process. The proposed model yields random realizations of maps of liquefaction and nonliquefaction conditioned on a map of [Formula: see text]. Such maps can be used to constrain the area over which displacements are estimated using soil properties inferred from geology and are therefore a critical component in reducing bias in assessments of liquefaction risk at the regional scale.

Feature-based Bathymetric Matching of Autonomous Underwater Vehicle Transects using Robust Gaussian Processes

Autonomous underwater vehicles have become essential tools for the collection of high-resolution bathymetric and geophysical data. An inertial navigation system, aided by a Doppler velocity log and a surface-mounted ultrashort baseline acoustic positioning system, is generally capable of providing sufficient accuracy at conventional survey scales. However, the accuracy decreases with depth, and the resulting relative positioning across different transects may not be sufficient to resolve features of interest at fine scales. This work presents a method for accurate coregistration of the position of adjacent transects. The approach is based on detecting and matching local features in overlapping multibeam echosounder swaths. The navigational errors for the transects are taken to be described by latent Gaussian processes, observed through these matches. The hyperparameters of the Gaussian process are learned from the data themselves and do neither require tuning of filter parameters nor intimate knowledge of the autonomous system or its sensor configuration. The proposed method is robust to outliers by considering a non-Gaussian observation model. The approach is demonstrated on a data set collected at the Arctic Mid-Ocean Ridge (AMOR). The method can be used to construct high-resolution bathymetric models from repeated passes over the same area or to accurately coregister other sensors such as cameras, subbottom profilers, and magnetometers. The primary contribution of this work is the application of feature-based matching of bathymetry with a robust Gaussian process.