Search

The integrated nested Laplace approximation (INLA) method has become a popular approach for computationally efficient approximate Bayesian computation. In particular, by leveraging sparsity in random effect precision matrices, INLA is commonly used in spatial and spatio-temporal applications. However, the speed of INLA comes at the cost of restricting the user to the family of latent Gaussian models and the likelihoods currently implemented in {INLA}, the main software implementation of the INLA methodology. {inlabru} is a software package that extends the types of models that can be fitted using INLA by allowing the latent predictor to be non-linear in its parameters, moving beyond the additive linear predictor framework to allow more complex functional relationships. For inference it uses an approximate iterative method based on the first-order Taylor expansion of the non-linear predictor, fitting the model using INLA for each linearised model configuration. {inlabru} automates much of the workflow required to fit models using {R-INLA}, simplifying the process for users to specify, fit and predict from models. There is additional support for fitting joint likelihood models by building each likelihood individually. {inlabru} also supports the direct use of spatial data structures, such as those implemented in the {sf} and {terra} packages. In this paper we outline the statistical theory, model structure and basic syntax required for users to understand and develop their own models using {inlabru}. We evaluate the approximate inference method using a Bayesian method checking approach. We provide three examples modelling simulated spatial data that demonstrate the benefits of the additional flexibility provided by {inlabru}.

This work aims to combine two primary meteorological data sources in the Philippines: data from a sparse network of weather stations and outcomes of a numerical weather prediction model. To this end, we propose a data fusion model which is primarily motivated by the problem of sparsity in the observational data and the use of a numerical prediction model as an additional data source in order to obtain better predictions for the variables of interest. The proposed data fusion model assumes that the different data sources are error-prone realizations of a common latent process. The outcomes from the weather stations follow the classical error model while the outcomes of the numerical weather prediction model involves a constant multiplicative bias parameter and an additive bias which is spatially-structured and time-varying. We use a Bayesian model averaging approach with the integrated nested Laplace approximation (INLA) for doing inference. The proposed data fusion model outperforms the stations-only model and the regression calibration approach, when assessed using leave-group-out cross-validation (LGOCV). We assess the benefits of data fusion and evaluate the accuracy of predictions and parameter estimation through a simulation study. The results show that the proposed data fusion model generally gives better predictions compared to the stations-only approach especially with sparse observational data.

In credit risk analysis, survival models with fixed and time-varying covariates are widely used to predict a borrower's time-to-event. When the time-varying drivers are endogenous, modelling jointly the evolution of the survival time and the endogenous covariates is the most appropriate approach, also known as the joint model for longitudinal and survival data. In addition to the temporal component, credit risk models can be enhanced when including borrowers' geographical information by considering spatial clustering and its variation over time. We propose the Spatio-Temporal Joint Model (STJM) to capture spatial and temporal effects and their interaction. This Bayesian hierarchical joint model reckons the survival effect of unobserved heterogeneity among borrowers located in the same region at a particular time. To estimate the STJM model for large datasets, we consider the Integrated Nested Laplace Approximation (INLA) methodology. We apply the STJM to predict the time to full prepayment on a large dataset of 57,258 US mortgage borrowers with more than 2.5 million observations. Empirical results indicate that including spatial effects consistently improves the performance of the joint model. However, the gains are less definitive when we additionally include spatio-temporal interactions.

Questions about information encoded by the brain demand statistical frameworks for inferring relationships between neural firing and features of the world. The landmark discovery of grid cells demonstrates that neurons can represent spatial information through regularly repeating firing fields. However, the influence of covariates may be masked in current statistical models of grid cell activity, which by employing approaches such as discretizing, aggregating and smoothing, are computationally inefficient and do not account for the continuous nature of the physical world. These limitations motivated us to develop likelihood-based procedures for modelling and estimating the firing activity of grid cells conditionally on biologically relevant covariates. Our approach models firing activity using Poisson point processes with latent Gaussian effects, which accommodate persistent inhomogeneous spatial-directional patterns and overdispersion. Inference is performed in a fully Bayesian manner, which allows us to quantify uncertainty. Applying these methods to experimental data, we provide evidence for temporal and local head direction effects on grid firing. Our approaches offer a novel and principled framework for analysis of neural representations of space., Comment: 44 pages, 16 figures

The Epidemic Type Aftershock Sequence (ETAS) model is widely used to model seismic sequences and underpins Operational Earthquake Forecasting (OEF). However, it remains challenging to assess the reliability of inverted ETAS parameters for a range of reasons. The most common algorithms just return point estimates with little quantification of uncertainty, and Bayesian Markov Chain Monte Carlo implementations remain slow to run, do not scale well and few have been extended to include spatial structure. Here we present a new approach to ETAS modelling using an alternative Bayesian method, the Integrated Nested Laplace Approximation (INLA). We have implemented this model in a new R-Package called ETAS.inlabru, which builds on the R packages R-INLA and inlabru . Whilst we just present the temporal component here, the model scales to a spatio-temporal model and may include a variety of spatial covariates. Using a series of synthetic case studies, we explore the robustness of our ETAS inversion method. We demonstrate that reliable estimates of the model parameters require that the catalogue data contains periods of relative quiescence as well as triggered sequences. We explore the robustness under stochastic uncertainty in the training data and show that the method is robust to a wide range of starting conditions. We show how the inclusion of historic earthquakes prior to the modelled domain affects the quality of the inversion. Finally, we show that rate dependent incompleteness after large earthquakes has a significant and detrimental effect on the ETAS posteriors. We believe that the speed of the inlabru inversion, which include a rigorous estimation of uncertainty, will enable a deeper exploration of how to use ETAS robustly for seismicity modelling and operational earthquake forecasting.

Hawkes process are very popular mathematical tools for modelling phenomena exhibiting a \textit{self-exciting} or \textit{self-correcting} behaviour. Typical examples are earthquakes occurrence, wild-fires, drought, capture-recapture, crime violence, trade exchange, and social network activity. The widespread use of Hawkes process in different fields calls for fast, reproducible, reliable, easy-to-code techniques to implement such models. We offer a technique to perform approximate Bayesian inference of Hawkes process parameters based on the use of the R-package \inlabru. The \inlabru R-package, in turn, relies on the INLA methodology to approximate the posterior of the parameters. Our Hawkes process approximation is based on a decomposition of the log-likelihood in three parts, which are linearly approximated separately. The linear approximation is performed with respect to the mode of the parameters' posterior distribution, which is determined with an iterative gradient-based method. The approximation of the posterior parameters is therefore deterministic, ensuring full reproducibility of the results. The proposed technique only requires the user to provide the functions to calculate the different parts of the decomposed likelihood, which are internally linearly approximated by the R-package \inlabru. We provide a comparison with the \bayesianETAS R-package which is based on an MCMC method. The two techniques provide similar results but our approach requires two to ten times less computational time to converge, depending on the amount of data., Comment: 2o pages, 7 figures, 5 tables

Gaussian processes and random fields have a long history, covering multiple approaches to representing spatial and spatio-temporal dependence structures, such as covariance functions, spectral representations, reproducing kernel Hilbert spaces, and graph based models. This article describes how the stochastic partial differential equation approach to generalising Mat\'ern covariance models via Hilbert space projections connects with several of these approaches, with each connection being useful in different situations. In addition to an overview of the main ideas, some important extensions, theory, applications, and other recent developments are discussed. The methods include both Markovian and non-Markovian models, non-Gaussian random fields, non-stationary fields and space-time fields on arbitrary manifolds, and practical computational considerations., Comment: 33 pages, 1 figure

Diurnal temperature range is an important variable in climate science that can provide information regarding climate variability and climate change. Changes in diurnal temperature range can have implications for hydrology, human health and ecology, among others. Yet, the statistical literature on modelling diurnal temperature range is lacking. In this paper we propose to model the distribution of diurnal temperature range using the five-parameter lambda (FPL) distribution. Additionally, in order to model diurnal temperature range with explanatory variables, we propose a distributional quantile regression model that combines quantile regression with marginal modelling using the FPL distribution. Inference is performed using the method of quantiles. The models are fitted to 30 years of daily observations of diurnal temperature range from 112 weather stations in the southern part of Norway. The flexible FPL distribution shows great promise as a model for diurnal temperature range, and performs well against competing models. The distributional quantile regression model is fitted to diurnal temperature range data using geographic, orographic and climatological explanatory variables. It performs well and captures much of the spatial variation in the distribution of diurnal temperature range in Norway., Comment: 28 pages, 9 figures; v2: revision of the introduction, more references added and minor corrections of the text

Operational earthquake forecasting for risk management and communication during seismic sequences depends on our ability to select an optimal forecasting model. To do this, we need to compare the performance of competing models with each other in prospective forecasting mode, and to rank their performance using a fair, reproducible and reliable method. The Collaboratory for the Study of Earthquake Predictability (CSEP) conducts such prospective earthquake forecasting experiments around the globe. One metric that has been proposed to rank competing models is the Parimutuel Gambling score, which has the advantage of allowing alarm-based (categorical) forecasts to be compared with probabilistic ones. Here we examine the suitability of this score for ranking competing earthquake forecasts. First, we prove analytically that this score is in general improper, meaning that, on average, it does not prefer the model that generated the data. Even in the special case where it is proper, we show it can still be used in an improper way. Then, we compare its performance with two commonly-used proper scores (the Brier and logarithmic scores), taking into account the uncertainty around the observed average score. We estimate the confidence intervals for the expected score difference which allows us to define if and when a model can be preferred. Our findings suggest the Parimutuel Gambling score should not be used to distinguishing between multiple competing forecasts. They also enable a more rigorous approach to distinguish between the predictive skills of candidate forecasts in addition to their rankings., Comment: 29 pages, 14 figures. Work presented at vEGU21 as vPico presentation

Gaussian random fields with Mat\'ern covariance functions are popular models in spatial statistics and machine learning. In this work, we develop a spatio-temporal extension of the Gaussian Mat\'ern fields formulated as solutions to a stochastic partial differential equation. The spatially stationary subset of the models have marginal spatial Mat\'ern covariances, and the model also extends to Whittle-Mat\'ern fields on curved manifolds, and to more general non-stationary fields. In addition to the parameters of the spatial dependence (variance, smoothness, and practical correlation range) it additionally has parameters controlling the practical correlation range in time, the smoothness in time, and the type of non-separability of the spatio-temporal covariance. Through the separability parameter, the model also allows for separable covariance functions. We provide a sparse representation based on a finite element approximation, that is well suited for statistical inference and which is implemented in the R-INLA software. The flexibility of the model is illustrated in an application to spatio-temporal modeling of global temperature data., Comment: 40 pages, 10 figures

This paper introduces a Laplace approximation to Bayesian inference in Dirichlet regression models, which can be used to analyze a set of variables on a simplex exhibiting skewness and heteroscedasticity, without having to transform the data. These data, which mainly consist of proportions or percentages of disjoint categories, are widely known as compositional data and are common in areas such as ecology, geology, and psychology. We provide both the theoretical foundations and a description of how Laplace approximation can be implemented in the case of Dirichlet regression. The paper also introduces the package dirinla in the R-language that extends the R-INLA package, which can not deal directly with Dirichlet likelihoods. Simulation studies are presented to validate the good behaviour of the proposed method, while a real data case-study is used to show how this approach can be applied.

Bayesian whole-brain functional magnetic resonance imaging (fMRI) analysis with three-dimensional spatial smoothing priors has been shown to produce state-of-the-art activity maps without pre-smoothing the data. The proposed inference algorithms are computationally demanding however, and the proposed spatial priors have several less appealing properties, such as being improper and having infinite spatial range. We propose a statistical inference framework for whole-brain fMRI analysis based on the class of Mat\'ern covariance functions. The framework uses the Gaussian Markov random field (GMRF) representation of possibly anisotropic spatial Mat\'ern fields via the stochastic partial differential equation (SPDE) approach of Lindgren et al. (2011). This allows for more flexible and interpretable spatial priors, while maintaining the sparsity required for fast inference in the high-dimensional whole-brain setting. We develop an accelerated stochastic gradient descent (SGD) optimization algorithm for empirical Bayes (EB) inference of the spatial hyperparameters. Conditionally on the inferred hyperparameters, we make a fully Bayesian treatment of the brain activity. The Mat\'ern prior is applied to both simulated and experimental task-fMRI data and clearly demonstrates that it is a more reasonable choice than the previously used priors, using comparisons of activity maps, prior simulation and cross-validation.

Coming up with Bayesian models for spatial data is easy, but performing inference with them can be challenging. Writing fast inference code for a complex spatial model with realistically-sized datasets from scratch is time-consuming, and if changes are made to the model, there is little guarantee that the code performs well. The key advantages of R-INLA are the ease with which complex models can be created and modified, without the need to write complex code, and the speed at which inference can be done even for spatial problems with hundreds of thousands of observations. R-INLA handles latent Gaussian models, where fixed effects, structured and unstructured Gaussian random effects are combined linearly in a linear predictor, and the elements of the linear predictor are observed through one or more likelihoods. The structured random effects can be both standard areal model such as the Besag and the BYM models, and geostatistical models from a subset of the Mat\'ern Gaussian random fields. In this review, we discuss the large success of spatial modelling with R-INLA and the types of spatial models that can be fitted, we give an overview of recent developments for areal models, and we give an overview of the stochastic partial differential equation (SPDE) approach and some of the ways it can be extended beyond the assumptions of isotropy and separability. In particular, we describe how slight changes to the SPDE approach leads to straight-forward approaches for non-stationary spatial models and non-separable space-time models., Comment: Extensive update, restructuring of sections

The Gaussian process is an indispensable tool for spatial data analysts. The onset of the "big data" era, however, has lead to the traditional Gaussian process being computationally infeasible for modern spatial data. As such, various alternatives to the full Gaussian process that are more amenable to handling big spatial data have been proposed. These modern methods often exploit low rank structures and/or multi-core and multi-threaded computing environments to facilitate computation. This study provides, first, an introductory overview of several methods for analyzing large spatial data. Second, this study describes the results of a predictive competition among the described methods as implemented by different groups with strong expertise in the methodology. Specifically, each research group was provided with two training datasets (one simulated and one observed) along with a set of prediction locations. Each group then wrote their own implementation of their method to produce predictions at the given location and each which was subsequently run on a common computing environment. The methods were then compared in terms of various predictive diagnostics. Supplementary materials regarding implementation details of the methods and code are available for this article online.

The use of sparse precision (inverse covariance) matrices has become popular because they allow for efficient algorithms for joint inference in high-dimensional models. Many applications require the computation of certain elements of the covariance matrix, such as the marginal variances, which may be non-trivial to obtain when the dimension is large. This paper introduces a fast Rao-Blackwellized Monte Carlo sampling based method for efficiently approximating selected elements of the covariance matrix. The variance and confidence bounds of the approximations can be precisely estimated without additional computational costs. Furthermore, a method that iterates over subdomains is introduced, and is shown to additionally reduce the approximation errors to practically negligible levels in an application on functional magnetic resonance imaging data. Both methods have low memory requirements, which is typically the bottleneck for competing direct methods.

The R software package excursions contains methods for calculating probabilistic excursion sets, contour credible regions, and simultaneous confidence bands for latent Gaussian stochastic processes and fields. It also contains methods for uncertainty quantification of contour maps and computation of Gaussian integrals. This article describes the theoretical and computational methods used in the package. The main functions of the package are introduced and two examples illustrate how the package can be used.

The key operation in Bayesian inference, is to compute high-dimensional integrals. An old approximate technique is the Laplace method or approximation, which dates back to Pierre- Simon Laplace (1774). This simple idea approximates the integrand with a second order Taylor expansion around the mode and computes the integral analytically. By developing a nested version of this classical idea, combined with modern numerical techniques for sparse matrices, we obtain the approach of Integrated Nested Laplace Approximations (INLA) to do approximate Bayesian inference for latent Gaussian models (LGMs). LGMs represent an important model-abstraction for Bayesian inference and include a large proportion of the statistical models used today. In this review, we will discuss the reasons for the success of the INLA-approach, the R-INLA package, why it is so accurate, why the approximations are very quick to compute and why LGMs make such a useful concept for Bayesian computing., Comment: 28 pages, 7 figures

Discussion of "Cross-Covariance Functions for Multivariate Geostatistics" by Genton and Kleiber [arXiv:1507.08017]., Comment: Published at http://dx.doi.org/10.1214/15-STS515 in the Statistical Science (http://www.imstat.org/sts/) by the Institute of Mathematical Statistics (http://www.imstat.org)

Contour maps are widely used to display estimates of spatial fields. Instead of showing the estimated field, a contour map only shows a fixed number of contour lines for different levels. However, despite the ubiquitous use of these maps, the uncertainty associated with them has been given a surprisingly small amount of attention. We derive measures of the statistical uncertainty, or quality, of contour maps, and use these to decide an appropriate number of contour lines, that relates to the uncertainty in the estimated spatial field. For practical use in geostatistics and medical imaging, computational methods are constructed, that can be applied to Gaussian Markov random fields, and in particular be used in combination with integrated nested Laplace approximations for latent Gaussian models. The methods are demonstrated on simulated data and an application to temperature estimation is presented.

Priors are important for achieving proper posteriors with physically meaningful covariance structures for Gaussian random fields (GRFs) since the likelihood typically only provides limited information about the covariance structure under in-fill asymptotics. We extend the recent Penalised Complexity prior framework and develop a principled joint prior for the range and the marginal variance of one-dimensional, two-dimensional and three-dimensional Mat\'ern GRFs with fixed smoothness. The prior is weakly informative and penalises complexity by shrinking the range towards infinity and the marginal variance towards zero. We propose guidelines for selecting the hyperparameters, and a simulation study shows that the new prior provides a principled alternative to reference priors that can leverage prior knowledge to achieve shorter credible intervals while maintaining good coverage. We extend the prior to a non-stationary GRF parametrized through local ranges and marginal standard deviations, and introduce a scheme for selecting the hyperparameters based on the coverage of the parameters when fitting simulated stationary data. The approach is applied to a dataset of annual precipitation in southern Norway and the scheme for selecting the hyperparameters leads to concervative estimates of non-stationarity and improved predictive performance over the stationary model.

Estimation of stationary dependence structure parameters using only a single realisation of the spatial process, typically leads to inaccurate estimates and poorly identified parameters. A common way to handle this is to fix some of the parameters, or within the Bayesian framework, impose prior knowledge. In many applied settings, stationary models are not flexible enough to model the process of interest, thus non-stationary spatial models are used. However, more flexible models usually means more parameters, and the identifiability problem becomes even more challenging. We investigate aspects of estimation of a Bayesian non-stationary spatial model for annual precipitation using observations from multiple years. The model contains replicates of the spatial field, which increases precision of the estimates and makes them less prior sensitive. Using R-INLA, we analyse precipitation data from southern Norway, and investigate statistical properties of the replicate model in a simulation study. The non-stationary spatial model we explore belongs to a recently introduced class of stochastic partial differential equation (SPDE) based spatial models. This model class allows for non-stationary models with explanatory variables in the dependence structure. We derive conditions to facilitate prior specification for these types of non-stationary spatial models.

A stationary spatial model is an idealization and we expect that the true dependence structures of physical phenomena are spatially varying, but how should we handle this non-stationarity in practice? We study the challenges involved in applying a flexible non-stationary model to a dataset of annual precipitation in the conterminous US, where exploratory data analysis shows strong evidence of a non-stationary covariance structure. The aim of this paper is to investigate the modelling pipeline once non-stationarity has been detected in spatial data. We show that there is a real danger of over-fitting the model and that careful modelling is necessary in order to properly account for varying second-order structure. In fact, the example shows that sometimes non-stationary Gaussian random fields are not necessary to model non-stationary spatial data., Comment: Minor change from previous version. arXiv admin note: text overlap with arXiv:1306.0408

This paper propose a novel decomposable graphical model to accommodate skew Gaussian graphical models. We encode conditional independence structure among the components of the multivariate closed skew normal random vector by means of a decomposable graph and so that the pattern of zero off-diagonal elements in the precision matrix corresponds to the missing edges of the given graph. We present conditions that guarantee the propriety of the posterior distributions under the standard noninformative priors for mean vector and precision matrix, and a proper prior for skewness parameter. The identifiability of the parameters is investigated by a simulation study. Finally, we apply our methodology to two data sets.

In this paper we propose a new approach for constructing \emph{multivariate} Gaussian random fields (GRFs) with oscillating covariance functions through systems of stochastic partial differential equations (SPDEs). We discuss how to build systems of SPDEs that introduces oscillation characteristics in the covariance functions of the multivariate GRFs. By choosing different parametrization of the equations, some GRFs can be made with oscillating covariance functions but other fields can have Mat\'ern covariance functions or close to Mat\'ern covariance functions. The multivariate GRFs constructed by solving the systems of SPDEs automatically fulfill the hard requirement of nonnegative definiteness for the covariance functions. The approximate weak solutions to the systems of SPDEs are used to represent the multivariate GRFs by multivariate Gaussian \emph{Markov} random fields (GMRFs). Since the multivariate GMRFs have sparse precision matrices (inverse of the covariance matrices), numerical algorithms for sparse matrices can be applied to the precision matrices for sampling and inference. Thus from a computational point of view, the \emph{big-n} problem can be partially solved with these types of models. Another advantage of the method is that the oscillation in the covariance function can be controlled directly by the parameters in the system of SPDEs. We show how to use this proposed approach with simulated data and real data examples., Comment: 40 pages, 22 figures

In this paper a new approach for constructing \emph{multivariate} Gaussian random fields (GRFs) using systems of stochastic partial differential equations (SPDEs) has been introduced and applied to simulated data and real data. By solving a system of SPDEs, we can construct multivariate GRFs. On the theoretical side, the notorious requirement of non-negative definiteness for the covariance matrix of the GRF is satisfied since the constructed covariance matrices with this approach are automatically symmetric positive definite. Using the approximate stochastic weak solutions to the systems of SPDEs, multivariate GRFs are represented by multivariate Gaussian \emph{Markov} random fields (GMRFs) with sparse precision matrices. Therefore, on the computational side, the sparse structures make it possible to use numerical algorithms for sparse matrices to do fast sampling from the random fields and statistical inference. Therefore, the \emph{big-n} problem can also be partially resolved for these models. These models out-preform existing multivariate GRF models on a commonly used real dataset., Comment: 47 pages, 19 figures

A non-stationary spatial Gaussian random field (GRF) is described as the solution of an inhomogeneous stochastic partial differential equation (SPDE), where the covariance structure of the GRF is controlled by the coefficients in the SPDE. This allows for a flexible way to vary the covariance structure, where intuition about the resulting structure can be gained from the local behaviour of the differential equation. Additionally, computations can be done with computationally convenient Gaussian Markov random fields which approximate the true GRFs. The model is applied to a dataset of annual precipitation in the conterminous US. The non-stationary model performs better than a stationary model measured with both CRPS and the logarithmic scoring rule.

Gaussian random fields (GRFs) constitute an important part of spatial modelling, but can be computationally infeasible for general covariance structures. An efficient approach is to specify GRFs via stochastic partial differential equations (SPDEs) and derive Gaussian Markov random field (GMRF) approximations of the solutions. We consider the construction of a class of non-stationary GRFs with varying local anisotropy, where the local anisotropy is introduced by allowing the coefficients in the SPDE to vary with position. This is done by using a form of diffusion equation driven by Gaussian white noise with a spatially varying diffusion matrix. This allows for the introduction of parameters that control the GRF by parametrizing the diffusion matrix. These parameters and the GRF may be considered to be part of a hierarchical model and the parameters estimated in a Bayesian framework. The results show that the use of an SPDE with non-constant coefficients is a promising way of creating non-stationary spatial GMRFs that allow for physical interpretability of the parameters, although there are several remaining challenges that would need to be solved before these models can be put to general practical use.

An interesting statistical problem is to find regions where some studied process exceeds a certain level. Estimating such regions so that the probability for exceeding the level in the entire set is equal to some predefined value is a difficult problem that occurs in several areas of applications ranging from brain imaging to astrophysics. In this work, a method for solving this problem, as well as the related problem of finding uncertainty regions for contour curves, for latent Gaussian models is proposed. The method is based on using a parametric family for the excursion sets in combination with a sequential importance sampling method for estimating joint probabilities. The accuracy of the method is investigated using simulated data and two environmental applications are presented. In the first application, areas where the air pollution in the Piemonte region in northern Italy exceeds the daily limit value, set by the European Union for human health protection, are estimated. In the second application, regions in the African Sahel that experienced an increase in vegetation after the drought period in the early 1980s are estimated.

The INLA approach for approximate Bayesian inference for latent Gaussian models has been shown to give fast and accurate estimates of posterior marginals and also to be a valuable tool in practice via the R-package R-INLA. In this paper we formalize new developments in the R-INLA package and show how these features greatly extend the scope of models that can be analyzed by this interface. We also discuss the current default method in R-INLA to approximate posterior marginals of the hyperparameters using only a modest number of evaluations of the joint posterior distribution of the hyperparameters, without any need for numerical integration.

The smoothing spline is one of the most popular curve-fitting methods, partly because of empirical evidence supporting its effectiveness and partly because of its elegant mathematical formulation. However, there are two obstacles that restrict the use of smoothing spline in practical statistical work. Firstly, it becomes computationally prohibitive for large data sets because the number of basis functions roughly equals the sample size. Secondly, its global smoothing parameter can only provide constant amount of smoothing, which often results in poor performances when estimating inhomogeneous functions. In this work, we introduce a class of adaptive smoothing spline models that is derived by solving certain stochastic differential equations with finite element methods. The solution extends the smoothing parameter to a continuous data-driven function, which is able to capture the change of the smoothness of underlying process. The new model is Markovian, which makes Bayesian computation fast. A simulation study and real data example are presented to demonstrate the effectiveness of our method., Comment: 26 Pages, 3 Figures

This paper introduces a new method for performing computational inference on log-Gaussian Cox processes. The likelihood is approximated directly by making novel use of a continuously specified Gaussian random field. We show that for sufficiently smooth Gaussian random field prior distributions, the approximation can converge with arbitrarily high order, while an approximation based on a counting process on a partition of the domain only achieves first-order convergence. The given results improve on the general theory of convergence of the stochastic partial differential equation models, introduced by Lindgren et al. (2011). The new method is demonstrated on a standard point pattern data set and two interesting extensions to the classical log-Gaussian Cox process framework are discussed. The first extension considers variable sampling effort throughout the observation window and implements the method of Chakraborty et al. (2011). The second extension constructs a log-Gaussian Cox process on the world's oceans. The analysis is performed using integrated nested Laplace approximation for fast approximate inference., Comment: 22 Pages, 8 figures

Gaussian Markov random fields (GMRFs) are frequently used as computationally efficient models in spatial statistics. Unfortunately, it has traditionally been difficult to link GMRFs with the more traditional Gaussian random field models as the Markov property is difficult to deploy in continuous space. Following the pioneering work of Lindgren et al. (2011), we expound on the link between Markovian Gaussian random fields and GMRFs. In particular, we discuss the theoretical and practical aspects of fast computation with continuously specified Markovian Gaussian random fields, as well as the clear advantages they offer in terms of clear, parsimonious and interpretable models of anisotropy and non-stationarity., Comment: 15 Pages, 5 Figures; 9/2011, Department of Mathematical Sciences, Norwegian University of Science and Technology (NTNU)

The Mat\'ern covariance function is a popular choice for modeling dependence in spatial environmental data. Standard Mat\'ern covariance models are, however, often computationally infeasible for large data sets. In this work, recent results for Markov approximations of Gaussian Mat\'{e}rn fields based on Hilbert space approximations are extended using wavelet basis functions. These Markov approximations are compared with two of the most popular methods for efficient covariance approximations; covariance tapering and the process convolution method. The results show that, for a given computational cost, the Markov methods have a substantial gain in accuracy compared with the other methods., Comment: Updated title and revised Section 4 to clarify the simulation setup

Latent Gaussian models are an extremely popular, flexible class of models. Bayesian inference for these models is, however, tricky and time consuming. Recently, Rue, Martino and Chopin introduced the Integrated Nested Laplace Approximation (INLA) method for deterministic fast approximate inference. In this paper, we outline the INLA approximation and its related R package. We will discuss the newer components of the r-INLA program as well as some possible extensions., Comment: 8 Pages, 2 Figures. Presented at ISI 2011

A new class of stochastic field models is constructed using nested stochastic partial differential equations (SPDEs). The model class is computationally efficient, applicable to data on general smooth manifolds, and includes both the Gaussian Mat\'{e}rn fields and a wide family of fields with oscillating covariance functions. Nonstationary covariance models are obtained by spatially varying the parameters in the SPDEs, and the model parameters are estimated using direct numerical optimization, which is more efficient than standard Markov Chain Monte Carlo procedures. The model class is used to estimate daily ozone maps using a large data set of spatially irregular global total column ozone data., Comment: Published in at http://dx.doi.org/10.1214/10-AOAS383 the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)