Bayesian Adaptive Regression Splines Research Articles

Bayesian optimization with adaptive surrogate models for automated experimental design

Bayesian optimization (BO) is an indispensable tool to optimize objective functions that either do not have known functional forms or are expensive to evaluate. Currently, optimal experimental design is always conducted within the workflow of BO leading to more efficient exploration of the design space compared to traditional strategies. This can have a significant impact on modern scientific discovery, in particular autonomous materials discovery, which can be viewed as an optimization problem aimed at looking for the maximum (or minimum) point for the desired materials properties. The performance of BO-based experimental design depends not only on the adopted acquisition function but also on the surrogate models that help to approximate underlying objective functions. In this paper, we propose a fully autonomous experimental design framework that uses more adaptive and flexible Bayesian surrogate models in a BO procedure, namely Bayesian multivariate adaptive regression splines and Bayesian additive regression trees. They can overcome the weaknesses of widely used Gaussian process-based methods when faced with relatively high-dimensional design space or non-smooth patterns of objective functions. Both simulation studies and real-world materials science case studies demonstrate their enhanced search efficiency and robustness.

Flexible covariate representations for extremes

AbstractEnvironmental extremes often show systematic variation with covariates. Three different nonparametric descriptions (penalized B‐splines, Bayesian adaptive regression splines, and Voronoi partition) for the dependence of extreme value model parameters on covariates are considered. These descriptions take the generic form of a linear combination of basis functions on the covariate domain, but differ (i) in the way that basis functions are constructed and possibly modified, and potentially (ii) by additional penalization of the variability (e.g., variance or roughness) of basis coefficients, for a given sample, to improve inference. The three representations are used to characterize variation of parameters in a nonstationary generalized Pareto model for the magnitude of threshold exceedances with respect to covariates. Computationally efficient schemes for Bayesian inference are used, including Riemann manifold Metropolis‐adjusted Langevin algorithm and reversible jump. A simulation study assesses relative performance of the three descriptions in estimating the distribution of the T‐year maximum event (for arbitrary T greater than the period of the sample) from a peaks over threshold extreme value analysis with respect to a single periodic covariate. The three descriptions are also used to estimate a directional tail model for peaks over threshold of storm peak significant wave height at a location in the northern North Sea.

BASS: An R Package for Fitting and Performing Sensitivity Analysis of Bayesian Adaptive Spline Surfaces

We present the R package BASS as a tool for nonparametric regression. The primary focus of the package is fitting fully Bayesian adaptive spline surface (BASS) models and performing global sensitivity analyses of these models. The BASS framework is similar to that of Bayesian multivariate adaptive regression splines (BMARS) from Denison, Mallick, and Smith (1998), but with many added features. The software is built to efficiently handle significant amounts of data with many continuous or categorical predictors and with functional response. Under our Bayesian framework, most priors are automatic but these can be modified by the user to focus on parsimony and the avoidance of overfitting. If directed to do so, the software uses parallel tempering to improve the reversible jump Markov chain Monte Carlo (RJMCMC) methods used to perform inference. We discuss the implementation of these features and present the performance of BASS in a number of analyses of simulated and real data.

Estimation of neuronal firing rate using Bayesian Adaptive Kernel Smoother (BAKS).

Neurons use sequences of action potentials (spikes) to convey information across neuronal networks. In neurophysiology experiments, information about external stimuli or behavioral tasks has been frequently characterized in term of neuronal firing rate. The firing rate is conventionally estimated by averaging spiking responses across multiple similar experiments (or trials). However, there exist a number of applications in neuroscience research that require firing rate to be estimated on a single trial basis. Estimating firing rate from a single trial is a challenging problem and current state-of-the-art methods do not perform well. To address this issue, we develop a new method for estimating firing rate based on a kernel smoothing technique that considers the bandwidth as a random variable with prior distribution that is adaptively updated under an empirical Bayesian framework. By carefully selecting the prior distribution together with Gaussian kernel function, an analytical expression can be achieved for the kernel bandwidth. We refer to the proposed method as Bayesian Adaptive Kernel Smoother (BAKS). We evaluate the performance of BAKS using synthetic spike train data generated by biologically plausible models: inhomogeneous Gamma (IG) and inhomogeneous inverse Gaussian (IIG). We also apply BAKS to real spike train data from non-human primate (NHP) motor and visual cortex. We benchmark the proposed method against established and previously reported methods. These include: optimized kernel smoother (OKS), variable kernel smoother (VKS), local polynomial fit (Locfit), and Bayesian adaptive regression splines (BARS). Results using both synthetic and real data demonstrate that the proposed method achieves better performance compared to competing methods. This suggests that the proposed method could be useful for understanding the encoding mechanism of neurons in cognitive-related tasks. The proposed method could also potentially improve the performance of brain-machine interface (BMI) decoder that relies on estimated firing rate as the input.

Two-Stage Metropolis-Hastings for Tall Data.

This paper discusses the challenges presented by tall data problems associated with Bayesian classification (specifically binary classification) and the existing methods to handle them. Current methods include parallelizing the likelihood, subsampling, and consensus Monte Carlo. A new method based on the two-stage Metropolis-Hastings algorithm is also proposed. The purpose of this algorithm is to reduce the exact likelihood computational cost in the tall data situation. In the first stage, a new proposal is tested by the approximate likelihood based model. The full likelihood based posterior computation will be conducted only if the proposal passes the first stage screening. Furthermore, this method can be adopted into the consensus Monte Carlo framework. The two-stage method is applied to logistic regression, hierarchical logistic regression, and Bayesian multivariate adaptive regression splines.

Task Performance Changes the Amplitude and Timing of the BOLD Signal.

Translational studies comparing imaging data of animals and humans have gained increasing scientific interests. With this upcoming translational approach, however, identifying harmonized statistical analysis as well as shared data acquisition protocols and/or combined statistical approaches is necessary. Following this idea, we applied Bayesian Adaptive Regression Splines (BARS), which have until now mainly been used to model neural responses of electrophysiological recordings from rodent data, on human hemodynamic responses as measured via fMRI. Forty-seven healthy subjects were investigated while performing the Attention Network Task in the MRI scanner. Fluctuations in the amplitude and timing of the BOLD response were determined and validated externally with brain activation using GLM and also ecologically with the influence of task performance (i.e. good vs. bad performers). In terms of brain activation, bad performers presented reduced activation bilaterally in the parietal lobules, right prefrontal cortex (PFC) and striatum. This was accompanied by an enhanced left PFC recruitment. With regard to the amplitude of the BOLD-signal, bad performers showed enhanced values in the left PFC. In addition, in the regions of reduced activation such as the parietal and striatal regions, the temporal dynamics were higher in bad performers. Based on the relation between BOLD response and neural firing with the amplitude of the BOLD signal reflecting gamma power and timing dynamics beta power, we argue that in bad performers, an enhanced left PFC recruitment hints towards an enhanced functioning of gamma-band activity in a compensatory manner. This was accompanied by reduced parieto-striatal activity, associated with increased and potentially conflicting beta-band activity.

Quantifying Strength of Chaos in the Population Firing Rate of Neurons.

In this study, I considered quantifying the strength of chaos in the population firing rate of a pulse-coupled neural network. In particular, I considered the dynamics where the population firing rate is chaotic and the firing of each neuron is stochastic. I calculated a time histogram of firings to show the variation in the population firing rate over time. To smooth this histogram, I used Bayesian adaptive regression splines and a gaussian filter. The nonlinear prediction method, based on reconstruction, was applied to a sequence of interpeak intervals in the smoothed time histogram of firings. I propose the use of the sum of nonlinearity as a quantifier of the strength of chaos. When applying this method to the firings of a pulse-coupled neural network, the sum of nonlinearity was seen to satisfy three properties for quantifying the strength of chaos. First, it can be calculated from spiking data alone. Second, it takes large values when applied to firings that are confirmed, theoretically or numerically, to be chaotic. Third, it reflects the strength of chaos of the original dynamics.

Emulation-Based Calibration for Parameters in Parameterized Phonon Spectrum of ZrHx in TRIGA Reactor Simulations

We investigate the calibration of the uncertainties of thermal scattering of ZrHx in the fuel material in TRIGA reactor simulations. Thermal scattering cross sections of ZrHx are heavily affected by the solid-state frequency distributions, also called phonon spectra. In previous work, we proposed parameterized phonon spectrum models and explored the effects on quantities of interest (QoIs) of changing spectra with such models by varying the parameters. In this work, we establish a more general calibration framework for the phonon spectrum of ZrHx. To accomplish this calibration, we introduce two emulators, Gaussian process regression and Bayesian multivariate adaptive regression splines, to create a map from the input parameters to the QoIs into the calibration framework. Using these emulators, we perform calibrations using the emulation results with the same QoIs at 600 K. Test simulations using data generated with calibrated parameters show that uncertainties of the QoIs shrink over 50%. Moreover, we extend the test to the reactivity at a different temperature, 293.6 K, as an extrapolated test of the calibration, and obtained results close to those of the surrogate experiment. The efficacy and efficiency of implementing emulators in the calibration framework are demonstrated.

Mixtures of Self-Modelling Regressions

A shape invariant model for functions f1,…,fn specifies that each individual function fi can be related to a common shape function g through the relation fi(x)=aig(cix + di) + bi. We consider a flexible mixture model that allows multiple shape functions g1,…,gK, where each fi is a shape invariant transformation of one of those gk. We derive an MCMC algorithm for fitting the model using Bayesian Adaptive Regression Splines (BARS), propose a strategy to improve its mixing properties and utilize existing model selection criteria in semiparametric mixtures to select the number of distinct shape functions. We discuss some of the computational difficulties that arise. The method is illustrated using synaptic transmission data, where the groups of functions may indicate different active zones in a synapse.

Multiscale analysis of neural spike trains

This paper studies the multiscale analysis of neural spike trains, through both graphical and Poisson process approaches. We introduce the interspike interval plot, which simultaneously visualizes characteristics of neural spiking activity at different time scales. Using an inhomogeneous Poisson process framework, we discuss multiscale estimates of the intensity functions of spike trains. We also introduce the windowing effect for two multiscale methods. Using quasi-likelihood, we develop bootstrap confidence intervals for the multiscale intensity function. We provide a cross-validation scheme, to choose the tuning parameters, and study its unbiasedness. Studying the relationship between the spike rate and the stimulus signal, we observe that adjusting for the first spike latency is important in cross-validation. We show, through examples, that the correlation between spike trains and spike count variability can be multiscale phenomena. Furthermore, we address the modeling of the periodicity of the spike trains caused by a stimulus signal or by brain rhythms. Within the multiscale framework, we introduce intensity functions for spike trains with multiplicative and additive periodic components. Analyzing a dataset from the retinogeniculate synapse, we compare the fit of these models with the Bayesian adaptive regression splines method and discuss the limitations of the methodology. Computational efficiency, which is usually a challenge in the analysis of spike trains, is one of the highlights of these new models. In an example, we show that the reconstruction quality of a complex intensity function demonstrates the ability of the multiscale methodology to crack the neural code.

A calibration and data assimilation method using the Bayesian MARS emulator

We present a method for calibrating the uncertain inputs to a computer model using available experimental data. The goal of the procedure is to estimate the posterior distribution of the uncertain inputs such that when samples from the posterior are used as inputs to future model runs, the model is more likely to replicate (or predict) the experimental response. The calibration is performed by sampling the space of the uncertain inputs, using the computer model (or, more likely, an emulator for the computer model) to assign weights to the samples, and applying the weights to produce the posterior distributions and generate predictions of new experiments with confidence bounds. The method is similar to Metropolis–Hastings calibration methods with independently sampled updates, except that we generate samples beforehand and replace the candidate acceptance routine with a weighting scheme.We apply our method to the calibration of a Hyades 2D model of laser energy deposition in beryllium. We employ a Bayesian Multivariate Adaptive Regression Splines (BMARS) emulator as a surrogate for Hyades 2D. We treat a range of uncertainties in our application, including uncertainties in the experimental inputs, experimental measurement error, and systematic experimental timing errors. The resulting posterior distributions agree with our existing intuition, and we validate the results by performing a series of leave-one-out predictions. We find that the calibrated predictions are considerably more accurate and less uncertain than blind sampling of the forward model alone.

The use of Bayesian nonlinear regression techniques for the modelling of the retention behaviour of volatile components of Artemisia species

The main aim of this work was to assess the ability of Bayesian multivariate adaptive regression splines (BMARS) and Bayesian radial basis function (BRBF) techniques for modelling the gas chromatographic retention indices of volatile components of Artemisia species. A diverse set of molecular descriptors was calculated and used as descriptor pool for modelling the retention indices. The ability of BMARS and BRBF techniques was explored for the selection of the most relevant descriptors and proper basis functions for modelling. The results revealed that BRBF technique is more reproducible than BMARS for modelling the retention indices and can be used as a method for variable selection and modelling in quantitative structure–property relationship (QSPR) studies. It is also concluded that the Markov chain Monte Carlo (MCMC) search engine, implemented in BRBF algorithm, is a suitable method for selecting the most important features from a vast number of them. The values of correlation between the calculated retention indices and the experimental ones for the training and prediction sets (0.935 and 0.902, respectively) revealed the prediction power of the BRBF model in estimating the retention index of volatile components of Artemisia species.

Comparison of two populations of curves with an application in neuronal data analysis

Often in neurophysiological studies, scientists are interested in testing hypotheses regarding the equality of the overall intensity functions of a group of neurons when recorded under two different experimental conditions. In this paper, we consider such a hypothesis testing problem. We propose two test statistics: a parametric test similar to the modified Hotelling's T2 statistic of Behseta and Kass (Statist. Med. 2005; 24:3523–3534), as well as a nonparametric one similar to the spatial signed-rank test statistic of Möttönen and Oja (J. Nonparametric Statist. 1995; 5:201–213). We implement these tests on smooth curves obtained via fitting Bayesian Adaptive Regression Splines (BARS) to the intensity functions of neuronal Peri-Stimulus Time Histograms. Through simulation, we show that the powers of our proposed tests are extremely high even when the number of sampled neurons and the number of trials per neuron are small. Finally, we apply our methods on a group of motor cortex neurons recorded during a reaching task.

Signal Processing for Neural Spike Trains

Signal Processing for Neural Spike Trains

State-space algorithms for estimating spike rate functions.

The accurate characterization of spike firing rates including the determination of when changes in activity occur is a fundamental issue in the analysis of neurophysiological data. Here we describe a state-space model for estimating the spike rate function that provides a maximum likelihood estimate of the spike rate, model goodness-of-fit assessments, as well as confidence intervals for the spike rate function and any other associated quantities of interest. Using simulated spike data, we first compare the performance of the state-space approach with that of Bayesian adaptive regression splines (BARS) and a simple cubic spline smoothing algorithm. We show that the state-space model is computationally efficient and comparable with other spline approaches. Our results suggest both a theoretically sound and practical approach for estimating spike rate functions that is applicable to a wide range of neurophysiological data.

Bayesian Adaptive Regression Splines for Hierarchical Data

This article considers methodology for hierarchical functional data analysis, motivated by studies of reproductive hormone profiles in the menstrual cycle. Current methods standardize the cycle lengths and ignore the timing of ovulation within the cycle, both of which are biologically informative. Methods are needed that avoid standardization, while flexibly incorporating information on covariates and the timing of reference events, such as ovulation and onset of menses. In addition, it is necessary to account for within-woman dependency when data are collected for multiple cycles. We propose an approach based on a hierarchical generalization of Bayesian multivariate adaptive regression splines. Our formulation allows for an unknown set of basis functions characterizing the population-averaged and woman-specific trajectories in relation to covariates. A reversible jump Markov chain Monte Carlo algorithm is developed for posterior computation. Applying the methods to data from the North Carolina Early Pregnancy Study, we investigate differences in urinary progesterone profiles between conception and nonconception cycles.

Misinformation in the conjugate prior for the linear model with implications for free-knot spline modelling

In the conjugate prior for the normal linear model, the prior variance for the coefficients is a multiple of the error variance parameter. However, if the prior mean for the coefficients is poorly chosen, the posterior distribution of the model can be seriously distorted because of prior dependence between the coefficients and error variance. In particular, the error variance will be overestimated, as will the posterior variance of the coefficients. This occurs because the prior mean, which can be thought of as a weighted pseudo-observation, is an outlier with respect to the real observations. While this situation will be easily noticed and avoided in simple models, in more complicated models, the effect can be easily overlooked. The issue arises in the unit information (UI) prior, a conjugate prior in which the prior contributes information equal to that in one observation. In particular, a successful Bayesian nonparametric regression model - Bayesian Adaptive Regression Splines (BARS) - that relies on the UI prior for its model selection step suffers from this problem, and addressing the problem within the Bayesian paradigm alters the penalty on model dimensionality.

Testing equality of two functions using BARS

This article presents two methods of testing the hypothesis of equality of two functions H(0):f(1)(t)=f(2)(t) for all t, in a generalized non-parametric regression framework using a recently developed generalized non-parametric regression method called Bayesian adaptive regression splines (BARS). Of particular interest is the special case of testing equality of two Poisson process intensity functions lambda(1) (t)=lambda(2) (t), which arises frequently in neurophysiological applications. The first method uses Bayes factors, and the second method uses a modified Hotelling T(2) test. Both methods are applied to the analysis of 347 motor cortical neurons and, for certain choices of test criteria, the two methods lead to the same conclusions for all but 7 neurons. A small simulation study of power indicates that the Bayes factor can be somewhat more powerful in small samples. The T(2)-type test should be useful in screening large number of neurons for condition-related activity, while the Bayes factor will be especially helpful in assessing evidence in favour of H(0).

Spline-based non-parametric regression for periodic functions and its application to directional tuning of neurons

The activity of neurons in the brain often varies systematically with some quantitative feature of a stimulus or action. A well-known example is the tendency of the firing rates of neurons in the primary motor cortex to vary with the direction of a subject's arm or wrist movement. When this movement is constrained to vary in only two dimensions, the direction of movement may be characterized by an angle, and the neuronal firing rate can be written as a function of this angle. The firing rate function has traditionally been fit with a cosine, but recent evidence suggests that departures from cosine tuning occur frequently. We report here a new non-parametric regression method for fitting periodic functions and demonstrate its application to the fitting of neuronal data. The method is an extension of Bayesian adaptive regression splines (BARS) and applies both to normal and non-normal data, including Poisson data, which commonly arise in neuronal applications. We compare the new method to a periodic version of smoothing splines and some parametric alternatives and find the new method to be especially valuable when the smoothness of the periodic function varies unevenly across its domain.

Automatic correction of ocular artifacts in the EEG: a comparison of regression-based and component-based methods

A variety of procedures have been proposed to correct ocular artifacts in the electroencephalogram (EEG), including methods based on regression, principal components analysis (PCA) and independent component analysis (ICA). The current study compared these three methods, and it evaluated a modified regression approach using Bayesian adaptive regression splines to filter the electrooculogram (EOG) before computing correction factors. We applied each artifact correction procedure to real and simulated EEG data of varying epoch lengths and then quantified the impact of correction on spectral parameters of the EEG. We found that the adaptive filter improved regression-based artifact correction. An automated PCA method effectively reduced ocular artifacts and resulted in minimal spectral distortion, whereas ICA correction appeared to distort power between 5 and 20 Hz. In general, reducing the epoch length improved the accuracy of estimating spectral power in the alpha (7.5–12.5 Hz) and beta (12.5–19.5 Hz) bands, but it worsened the accuracy for power in the theta (3.5–7.5 Hz) band and distorted time domain features. Results supported the use of regression-based and PCA-based ocular artifact correction and suggested a need for further studies examining possible spectral distortion from ICA-based correction procedures.