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}.

Modern methods for quantifying and predicting species distribution play a crucial part in biodiversity conservation. Occupancy models are a popular choice for analyzing species occurrence data as they allow to separate the observational error induced by imperfect detection, and the sources of bias affecting the occupancy process. However, the spatial and temporal variation in occupancy not accounted for by environmental covariates is often ignored or modelled through simple spatial structures as the computational costs of fitting explicit spatio-temporal models is too high. In this work, we demonstrate how INLA may be used to fit complex occupancy models and how the R-INLA package can provide a user-friendly interface to make such complex models available to users. We show how occupancy models, provided some simplification on the detection process, can be framed as latent Gaussian models and benefit from the powerful INLA machinery. A large selection of complex modelling features, and random effect modelshave already been implemented in R-INLA. These become available for occupancy models, providing the user with an efficient and flexible toolbox. We illustrate how INLA provides a computationally efficient framework for developing and fitting complex occupancy models using two case studies. Through these, we show how different spatio-temporal models that include spatial-varying trends, smooth terms, and spatio-temporal random effects can be fitted. At the cost of limiting the complexity of the detection model, INLA can incorporate a range of complex structures in the process. INLA-based occupancy models provide an alternative framework to fit complex spatiotemporal occupancy models. The need for new and more flexible computationally approaches to fit such models makes INLA an attractive option for addressing complex ecological problems, and a promising area of research.

This paper introduces the R package INLAjoint, designed as a toolbox for fitting a diverse range of regression models addressing both longitudinal and survival outcomes. INLAjoint relies on the computational efficiency of the integrated nested Laplace approximations methodology, an efficient alternative to Markov chain Monte Carlo for Bayesian inference, ensuring both speed and accuracy in parameter estimation and uncertainty quantification. The package facilitates the construction of complex joint models by treating individual regression models as building blocks, which can be assembled to address specific research questions. Joint models are relevant in biomedical studies where the collection of longitudinal markers alongside censored survival times is common. They have gained significant interest in recent literature, demonstrating the ability to rectify biases present in separate modeling approaches such as informative censoring by a survival event or confusion bias due to population heterogeneity. We provide a comprehensive overview of the joint modeling framework embedded in INLAjoint with illustrative examples. Through these examples, we demonstrate the practical utility of INLAjoint in handling complex data scenarios encountered in biomedical research.

In this work, we present a new approach for constructing models for correlation matrices with a user-defined graphical structure. The graphical structure makes correlation matrices interpretable and avoids the quadratic increase of parameters as a function of the dimension. We suggest an automatic approach to define a prior using a natural sequence of simpler models within the Penalized Complexity framework for the unknown parameters in these models. We illustrate this approach with three applications: a multivariate linear regression of four biomarkers, a multivariate disease mapping, and a multivariate longitudinal joint modelling. Each application underscores our method's intuitive appeal, signifying a substantial advancement toward a more cohesive and enlightening model that facilitates a meaningful interpretation of correlation matrices.

Performing a Bayesian inference on large spatio-temporal models requires extracting inverse elements of large sparse precision matrices for marginal variances. Although direct matrix factorizations can be used for the inversion, such methods fail to scale well for distributed problems when run on large computing clusters. On the contrary, Krylov subspace methods for the selected inversion have been gaining traction. We propose a parallel hybrid approach based on domain decomposition, which extends the Rao-Blackwellized Monte Carlo estimator for distributed precision matrices. Our approach exploits the strength of Krylov subspace methods as global solvers and efficiency of direct factorizations as base case solvers to compute the marginal variances using a divide-and-conquer strategy. By introducing subdomain overlaps, one can achieve a greater accuracy at an increased computational effort with little to no additional communication. We demonstrate the speed improvements on both simulated models and a massive US daily temperature data., Comment: 17 pages, 7 figures

Model checking is essential to evaluate the adequacy of statistical models and the validity of inferences drawn from them. Particularly, hierarchical models such as latent Gaussian models (LGMs) pose unique challenges as it is difficult to check assumptions about the distribution of the latent parameters. Discrepancy measures are often used to quantify the degree to which a model fit deviates from the observed data. We construct discrepancy measures by (a) defining an alternative model with relaxed assumptions and (b) deriving the discrepancy measure most sensitive to discrepancies induced by this alternative model. We also promote a workflow for model criticism that combines model checking with subsequent robustness analysis. As a result, we obtain a general recipe to check assumptions in LGMs and the impact of these assumptions on the results. We demonstrate the ideas by assessing the latent Gaussianity assumption, a crucial but often overlooked assumption in LGMs. We illustrate the methods via examples utilising Stan and provide functions for easy usage of the methods for general models fitted through R-INLA., Comment: 40 pages, 21 figures

This paper aims to extend the Besag model, a widely used Bayesian spatial model in disease mapping, to a non-stationary spatial model for irregular lattice-type data. The goal is to improve the model's ability to capture complex spatial dependence patterns and increase interpretability. The proposed model uses multiple precision parameters, accounting for different intensities of spatial dependence in different sub-regions. We derive a joint penalized complexity prior for the flexible local precision parameters to prevent overfitting and ensure contraction to the stationary model at a user-defined rate. The proposed methodology can be used as a basis for the development of various other non-stationary effects over other domains such as time. An accompanying R package 'fbesag' equips the reader with the necessary tools for immediate use and application. We illustrate the novelty of the proposal by modeling the risk of dengue in Brazil, where the stationary spatial assumption fails and interesting risk profiles are estimated when accounting for spatial non-stationary.

Bayesian inference tasks continue to pose a computational challenge. This especially holds for spatial-temporal modeling where high-dimensional latent parameter spaces are ubiquitous. The methodology of integrated nested Laplace approximations (INLA) provides a framework for performing Bayesian inference applicable to a large subclass of additive Bayesian hierarchical models. In combination with the stochastic partial differential equations (SPDE) approach it gives rise to an efficient method for spatial-temporal modeling. In this work we build on the INLA-SPDE approach, by putting forward a performant distributed memory variant, INLA-DIST, for large-scale applications. To perform the arising computational kernel operations, consisting of Cholesky factorizations, solving linear systems, and selected matrix inversions, we present two numerical solver options, a sparse CPU-based library and a novel blocked GPU-accelerated approach which we propose. We leverage the recurring nonzero block structure in the arising precision (inverse covariance) matrices, which allows us to employ dense subroutines within a sparse setting. Both versions of INLA-DIST are highly scalable, capable of performing inference on models with millions of latent parameters. We demonstrate their accuracy and performance on synthetic as well as real-world climate dataset applications., Comment: 22 pages, 14 figures

Bayesian methods and software for spatial data analysis are generally now well established in the scientific community. Despite the wide application of spatial models, the analysis of multivariate spatial data using R-INLA has not been widely described in the existing literature. Therefore, the main objective of this article is to demonstrate that R-INLA is a convenient toolbox to analyse different types of multivariate spatial datasets. Additionally, this will be illustrated by analysing three datasets which are publicly available. Furthermore, the details and the R code of these analyses are provided to exemplify how to adjust multivariate spatial datasets with R-INLA., Comment: Submitted to the RJournal (19 pages and 6 figures)

This tutorial shows how various Bayesian survival models can be fitted using the integrated nested Laplace approximation in a clear, legible, and comprehensible manner using the INLA and INLAjoint R-packages. Such models include accelerated failure time, proportional hazards, mixture cure, competing risks, multi-state, frailty, and joint models of longitudinal and survival data, originally presented in the article "Bayesian survival analysis with BUGS" (Alvares et al., 2021). In addition, we illustrate the implementation of a new joint model for a longitudinal semicontinuous marker, recurrent events, and a terminal event. Our proposal aims to provide the reader with syntax examples for implementing survival models using a fast and accurate approximate Bayesian inferential approach.

Latent Gaussian models (LGMs) are perhaps the most commonly used class of models in statistical applications. Nevertheless, in areas ranging from longitudinal studies in biostatistics to geostatistics, it is easy to find datasets that contain inherently non-Gaussian features, such as sudden jumps or spikes, that adversely affect the inferences and predictions made from an LGM. These datasets require more general latent non-Gaussian models (LnGMs) that can handle these non-Gaussian features automatically. However, fast implementation and easy-to-use software are lacking, which prevent LnGMs from becoming widely applicable. In this paper, we derive variational Bayes algorithms for fast and scalable inference of LnGMs. The approximation leads to an LGM that downweights extreme events in the latent process, reducing their impact and leading to more robust inferences. It can be applied to a wide range of models, such as autoregressive processes for time series, simultaneous autoregressive models for areal data, and spatial Mat\'ern models. To facilitate Bayesian inference, we introduce the ngvb package, where LGMs implemented in R-INLA can be easily extended to LnGMs by adding a single line of code., Comment: 30 pages, 10 figures

There is a growing demand for performing larger-scale Bayesian inference tasks, arising from greater data availability and higher-dimensional model parameter spaces. In this work we present parallelization strategies for the methodology of integrated nested Laplace approximations (INLA), a popular framework for performing approximate Bayesian inference on the class of Latent Gaussian models. Our approach makes use of nested OpenMP parallelism, a parallel line search procedure using robust regression in INLA's optimization phase and the state-of-the-art sparse linear solver PARDISO. We leverage mutually independent function evaluations in the algorithm as well as advanced sparse linear algebra techniques. This way we can flexibly utilize the power of today's multi-core architectures. We demonstrate the performance of our new parallelization scheme on a number of different real-world applications. The introduction of parallelism leads to speedups of a factor 10 and more for all larger models. Our work is already integrated in the current version of the open-source R-INLA package, making its improved performance conveniently available to all users., Comment: 18 pages, 7 figures

Various computational challenges arise when applying Bayesian inference approaches to complex hierarchical models. Sampling-based inference methods, such as Markov Chain Monte Carlo strategies, are renowned for providing accurate results but with high computational costs and slow or questionable convergence. On the contrary, approximate methods like the Integrated Nested Laplace Approximation (INLA) construct a deterministic approximation to the univariate posteriors through nested Laplace Approximations. This method enables fast inference performance in Latent Gaussian Models, which encode a large class of hierarchical models. R-INLA software mainly consists of three strategies to compute all the required posterior approximations depending on the accuracy requirements. The Simplified Laplace approximation (SLA) is the most attractive because of its speed performance since it is based on a Taylor expansion up to order three of a full Laplace Approximation. Here we enhance the methodology by simplifying the computations necessary for the skewness and modal configuration. Then we propose an expansion up to order four and use the Extended Skew Normal distribution as a new parametric fit. The resulting approximations to the marginal posterior densities are more accurate than those calculated with the SLA, with essentially no additional cost., Comment: 22 pages, 11 figures

Modeling longitudinal and survival data jointly offers many advantages such as addressing measurement error and missing data in the longitudinal processes, understanding and quantifying the association between the longitudinal markers and the survival events and predicting the risk of events based on the longitudinal markers. A joint model involves multiple submodels (one for each longitudinal/survival outcome) usually linked together through correlated or shared random effects. Their estimation is computationally expensive (particularly due to a multidimensional integration of the likelihood over the random effects distribution) so that inference methods become rapidly intractable, and restricts applications of joint models to a small number of longitudinal markers and/or random effects. We introduce a Bayesian approximation based on the Integrated Nested Laplace Approximation algorithm implemented in the R package R-INLA to alleviate the computational burden and allow the estimation of multivariate joint models with fewer restrictions. Our simulation studies show that R-INLA substantially reduces the computation time and the variability of the parameter estimates compared to alternative estimation strategies. We further apply the methodology to analyze 5 longitudinal markers (3 continuous, 1 count, 1 binary, and 16 random effects) and competing risks of death and transplantation in a clinical trial on primary biliary cholangitis. R-INLA provides a fast and reliable inference technique for applying joint models to the complex multivariate data encountered in health research.

The normal inverse Gaussian (NIG) and generalized asymmetric Laplace (GAL) distributions can be seen as skewed and semi-heavy-tailed extensions of the Gaussian distribution. Models driven by these more flexible noise distributions are then regarded as flexible extensions of simpler Gaussian models. Inferential procedures tend to overestimate the degree of non-Gaussianity in the data and therefore we propose controlling the flexibility of these non-Gaussian models by adding sensible priors in the inferential framework that contract the model towards Gaussianity. In our venture to derive sensible priors, we also propose a new intuitive parameterization of the non-Gaussian models and discuss how to implement them efficiently in $Stan$. The methods are derived for a generic class of non-Gaussian models that include spatial Mat\'ern fields, autoregressive models for time series, and simultaneous autoregressive models for aerial data. The results are illustrated with a simulation study and geostatistics application, where priors that penalize model complexity were shown to lead to more robust estimation and give preference to the Gaussian model, while at the same time allowing for non-Gaussianity if there is sufficient evidence in the data., Comment: 34 pages, 12 figures; Added new section about Stan implementation and a cross-validation study in the application section, removed unecessary details

This paper describes the methodology used by the team RedSea in the data competition organized for EVA 2021 conference. We develop a novel two-part model to jointly describe the wildfire count data and burnt area data provided by the competition organizers with covariates. Our proposed methodology relies on the integrated nested Laplace approximation combined with the stochastic partial differential equation (INLA-SPDE) approach. In the first part, a binary non-stationary spatio-temporal model is used to describe the underlying process that determines whether or not there is wildfire at a specific time and location. In the second part, we consider a non-stationary model that is based on log-Gaussian Cox processes for positive wildfire count data, and a non-stationary log-Gaussian model for positive burnt area data. Dependence between the positive count data and positive burnt area data is captured by a shared spatio-temporal random effect. Our two-part modeling approach performs well in terms of the prediction score criterion chosen by the data competition organizers. Moreover, our model results show that surface pressure is the most influential driver for the occurrence of a wildfire, whilst surface net solar radiation and surface pressure are the key drivers for large numbers of wildfires, and temperature and evaporation are the key drivers of large burnt areas.

Statistical analysis based on quantile regression methods is more comprehensive, flexible, and less sensitive to outliers when compared to mean regression methods. When the link between different diseases are of interest, joint disease mapping is useful for measuring directional correlation between them. Most studies study this link through multiple correlated mean regressions. In this paper we propose a joint quantile regression framework for multiple diseases where different quantile levels can be considered. We are motivated by the theorized link between the presence of Malaria and the gene deficiency G6PD, where medical scientist have anecdotally discovered a possible link between high levels of G6PD and lower than expected levels of Malaria initially pointing towards the occurrence of G6PD inhibiting the occurrence of Malaria. This link cannot be investigated with mean regressions and thus the need for flexible joint quantile regression in a disease mapping framework. Our joint quantile disease mapping model can be used for linear and non-linear effects of covariates by stochastic splines, since we define it as a latent Gaussian model. We perform Bayesian inference of this model using the INLA framework embedded in the R software package INLA. Finally, we illustrate the applicability of model by analyzing the malaria and G6PD deficiency incidences in 21 African countries using linked quantiles of different levels.

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

Bayesian disease mapping, yet if undeniably useful to describe variation in risk over time and space, comes with the hurdle of prior elicitation on hard-to-interpret random effect precision parameters. We introduce a reparametrized version of the popular spatio-temporal interaction models, based on Kronecker product intrinsic Gaussian Markov Random Fields, that we name the variance partitioning (VP) model. The VP model includes a mixing parameter that balances the contribution of the main and interaction effects to the total (generalized) variance and enhances interpretability. The use of a penalized complexity prior on the mixing parameter aids in coding prior information in a intuitive way. We illustrate the advantages of the VP model using two case studies.

The popular Bayesian meta-analysis expressed by Bayesian normal-normal hierarchical model (NNHM) synthesizes knowledge from several studies and is highly relevant in practice. Moreover, NNHM is the simplest Bayesian hierarchical model (BHM), which illustrates problems typical in more complex BHMs. Until now, it has been unclear to what extent the data determines the marginal posterior distributions of the parameters in NNHM. To address this issue we computed the second derivative of the Bhattacharyya coefficient with respect to the weighted likelihood, defined the total empirical determinacy (TED), the proportion of the empirical determinacy of location to TED (pEDL), and the proportion of the empirical determinacy of spread to TED (pEDS). We implemented this method in the R package \texttt{ed4bhm} and considered two case studies and one simulation study. We quantified TED, pEDL and pEDS under different modeling conditions such as model parametrization, the primary outcome, and the prior. This clarified to what extent the location and spread of the marginal posterior distributions of the parameters are determined by the data. Although these investigations focused on Bayesian NNHM, the method proposed is applicable more generally to complex BHMs., Comment: 22 pages, 1 figure

We show that many machine-learning algorithms are specific instances of a single algorithm called the \emph{Bayesian learning rule}. The rule, derived from Bayesian principles, yields a wide-range of algorithms from fields such as optimization, deep learning, and graphical models. This includes classical algorithms such as ridge regression, Newton's method, and Kalman filter, as well as modern deep-learning algorithms such as stochastic-gradient descent, RMSprop, and Dropout. The key idea in deriving such algorithms is to approximate the posterior using candidate distributions estimated by using natural gradients. Different candidate distributions result in different algorithms and further approximations to natural gradients give rise to variants of those algorithms. Our work not only unifies, generalizes, and improves existing algorithms, but also helps us design new ones.

The generalised extreme value (GEV) distribution is a three parameter family that describes the asymptotic behaviour of properly renormalised maxima of a sequence of independent and identically distributed random variables. If the shape parameter $\xi$ is zero, the GEV distribution has unbounded support, whereas if $\xi$ is positive, the limiting distribution is heavy-tailed with infinite upper endpoint but finite lower endpoint. In practical applications, we assume that the GEV family is a reasonable approximation for the distribution of maxima over blocks, and we fit it accordingly. This implies that GEV properties, such as finite lower endpoint in the case $\xi>0$, are inherited by the finite-sample maxima, which might not have bounded support. This is particularly problematic when predicting extreme observations based on multiple and interacting covariates. To tackle this usually overlooked issue, we propose a blended GEV distribution, which smoothly combines the left tail of a Gumbel distribution (GEV with $\xi=0$) with the right tail of a Fr\'echet distribution (GEV with $\xi>0$) and, therefore, has unbounded support. Using a Bayesian framework, we reparametrise the GEV distribution to offer a more natural interpretation of the (possibly covariate-dependent) model parameters. Independent priors over the new location and spread parameters induce a joint prior distribution for the original location and scale parameters. We introduce the concept of property-preserving penalised complexity (P$^3$C) priors and apply it to the shape parameter to preserve first and second moments. We illustrate our methods with an application to NO$_2$ pollution levels in California, which reveals the robustness of the bGEV distribution, as well as the suitability of the new parametrisation and the P$^3$C prior framework., Comment: 19 pages, 3 figures

A new method is proposed for modelling the yearly maxima of sub-daily precipitation, with the aim of producing spatial maps of return level estimates. Yearly precipitation maxima are modelled using a Bayesian hierarchical model with a latent Gaussian field, with the blended generalised extreme value (bGEV) distribution used as a substitute for the more standard generalised extreme value (GEV) distribution. Inference is made less wasteful with a novel two-step procedure that performs separate modelling of the scale parameter of the bGEV distribution using peaks over threshold data. Fast inference is performed using integrated nested Laplace approximations (INLA) together with the stochastic partial differential equation (SPDE) approach, both implemented in R-INLA. Heuristics for improving the numerical stability of R-INLA with the GEV and bGEV distributions are also presented. The model is fitted to yearly maxima of sub-daily precipitation from the south of Norway, and is able to quickly produce high-resolution return level maps with uncertainty. The proposed two-step procedure provides an improved model fit over standard inference techniques when modelling the yearly maxima of sub-daily precipitation with the bGEV distribution.

The Integrated Nested Laplace Approximation (INLA) is a deterministic approach to Bayesian inference on latent Gaussian models (LGMs) and focuses on fast and accurate approximation of posterior marginals for the parameters in the models. Recently, methods have been developed to extend this class of models to those that can be expressed as conditional LGMs by fixing some of the parameters in the models to descriptive values. These methods differ in the manner descriptive values are chosen. This paper proposes to combine importance sampling with INLA (IS-INLA), and extends this approach with the more robust adaptive multiple importance sampling algorithm combined with INLA (AMIS-INLA). This paper gives a comparison between these approaches and existing methods on a series of applications with simulated and observed datasets and evaluates their performance based on accuracy, efficiency, and robustness. The approaches are validated by exact posteriors in a simple bivariate linear model; then, they are applied to a Bayesian lasso model, a Bayesian imputation of missing covariate values, and lastly, in parametric Bayesian quantile regression. The applications show that the AMIS-INLA approach, in general, outperforms the other methods, but the IS-INLA algorithm could be considered for faster inference when good proposals are available.

Two-part joint models for a longitudinal semicontinuous biomarker and a terminal event have been recently introduced based on frequentist estimation. The biomarker distribution is decomposed into a probability of positive value and the expected value among positive values. Shared random effects can represent the association structure between the biomarker and the terminal event. The computational burden increases compared to standard joint models with a single regression model for the biomarker. In this context, the frequentist estimation implemented in the R package frailtypack can be challenging for complex models (i.e., large number of parameters and dimension of the random effects). As an alternative, we propose a Bayesian estimation of two-part joint models based on the Integrated Nested Laplace Approximation (INLA) algorithm to alleviate the computational burden and fit more complex models. Our simulation studies confirm that INLA provides accurate approximation of posterior estimates and to reduced computation time and variability of estimates compared to frailtypack in the situations considered. We contrast the Bayesian and frequentist approaches in the analysis of two randomized cancer clinical trials (GERCOR and PRIME studies), where INLA has a reduced variability for the association between the biomarker and the risk of event. Moreover, the Bayesian approach was able to characterize subgroups of patients associated with different responses to treatment in the PRIME study. Our study suggests that the Bayesian approach using INLA algorithm enables to fit complex joint models that might be of interest in a wide range of clinical applications.

The class of autoregressive (AR) processes is extensively used to model temporal dependence in observed time series. Such models are easily available and routinely fitted using freely available statistical software like R. A potential caveat in analyzing short time series is that commonly applied estimators for the coefficients of AR processes are severely biased. This paper suggests a model-based approach for bias correction of well-known estimators for the coefficients of first and second-order stationary AR processes, taking the sampling distribution of the original estimator into account. This is achieved by modeling the relationship between the true and estimated AR coefficients using weighted orthogonal polynomial regression, fitted to a huge number of simulations. The finite-sample distributions of the new estimators are approximated using transformations of skew-normal densities and their properties are demonstrated by simulations and in the analysis of a real ecological data set. The new estimators are easily available in our accompanying R-package ARbiascorrect for time series of length n = 10, 11, ... , 50, where original estimates are found using exact or conditional maximum likelihood, Burg's method or the Yule-Walker equations., Comment: 21 pages, 11 figures

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

For incremental quantile estimators the step size and possibly other tuning parameters must be carefully set. However, little attention has been given on how to set these values in an online manner. In this article we suggest two novel procedures that address this issue. The core part of the procedures is to estimate the current tracking mean squared error (MSE). The MSE is decomposed in tracking variance and bias and novel and efficient procedures to estimate these quantities are presented. It is shown that estimation bias can be tracked by associating it with the portion of observations below the quantile estimates. The first procedure runs an ensemble of $L$ quantile estimators for wide range of values of the tuning parameters and typically around $L = 100$. In each iteration an oracle selects the best estimate by the guidance of the estimated MSEs. The second method only runs an ensemble of $L = 3$ estimators and thus the values of the tuning parameters need from time to time to be adjusted for the running estimators. The procedures have a low memory foot print of $8L$ and a computational complexity of $8L$ per iteration. The experiments show that the procedures are highly efficient and track quantiles with an error close to the theoretical optimum. The Oracle approach performs best, but comes with higher computational cost. The procedures were further applied to a massive real-life data stream of tweets and proofed real world applicability of them.

The concept of depth represents methods to measure how deep an arbitrary point is positioned in a dataset and can be seen as the opposite of outlyingness. It has proved very useful and a wide range of methods have been developed based on the concept. To address the well-known computational challenges associated with the depth concept, we suggest to estimate Tukey depth contours using recently developed incremental quantile estimators. The suggested algorithm can estimate depth contours when the dataset in known in advance, but also recursively update and even track Tukey depth contours for dynamically varying data stream distributions. Tracking was demonstrated in a real-life data example where changes in human activity was detected in real-time from accelerometer observations.

The integrated nested Laplace approximation (INLA) for Bayesian inference is an efficient approach to estimate the posterior marginal distributions of the parameters and latent effects of Bayesian hierarchical models that can be expressed as latent Gaussian Markov random fields (GMRF). The representation as a GMRF allows the associated software R-INLA to estimate the posterior marginals in a fraction of the time as typical Markov chain Monte Carlo algorithms. INLA can be extended by means of Bayesian model averaging (BMA) to increase the number of models that it can fit to conditional latent GMRF. In this paper we review the use of BMA with INLA and propose a new example on spatial econometrics models., Comment: Submitted to Econometrics (MDPI)

This paper presents the development of a spatial block-Nearest Neighbor Gaussian process (block-NNGP) for location-referenced large spatial data. The key idea behind this approach is to divide the spatial domain into several blocks which are dependent under some constraints. The cross-blocks capture the large-scale spatial dependence, while each block captures the small-scale spatial dependence. The resulting block-NNGP enjoys Markov properties reflected on its sparse precision matrix. It is embedded as a prior within the class of latent Gaussian models, thus Bayesian inference is obtained using the integrated nested Laplace approximation (INLA). The performance of the block-NNGP is illustrated on simulated examples and massive real data for locations in the order of $10^4$., Comment: 60 pages, 20 figures (including the ones in the Supplementary Material), 4 tables

Assessing water quality and recognizing its associated risks to human health and the broader environment is undoubtedly essential. Groundwater is widely used to supply water for drinking, industry, and agriculture purposes. The groundwater quality measurements vary for different climates and various human behaviors, and consequently, their spatial variability can be substantial. In this paper, we aim to analyze a groundwater dataset from the Golestan province, Iran, for November 2003 to November 2013. Our target response variable to monitor the quality of groundwater is the number of counts that the quality of water is good for a drink. Hence, we are facing spatial count data. Due to the ubiquity of over or underdispersion in count data, we propose a Bayesian hierarchical modeling approach based on the renewal theory that relates nonexponential waiting times between events and the distribution of the counts, relaxing the assumption of equidispersion at the cost of an additional parameter. Particularly, we extend the methodology for the analysis of spatial count data based on the gamma distribution assumption for waiting times. The model can be formulated as a latent Gaussian model, and therefore, we can carry out the fast computation by using the integrated nested Laplace approximation method. The analysis of the groundwater dataset and a simulation study show a significant improvement over both Poisson and negative binomial models., Comment: 24 pages, 6 figures

Estimation of quantiles is one of the most fundamental real-time analysis tasks. Most real-time data streams vary dynamically with time and incremental quantile estimators document state-of-the art performance to track quantiles of such data streams. However, most are not able to make joint estimates of multiple quantiles in a consistent manner, and estimates may violate the monotone property of quantiles. In this paper we propose the general concept of *conditional quantiles* that can extend incremental estimators to jointly track multiple quantiles. We apply the concept to propose two new estimators. Extensive experimental results, on both synthetic and real-life data, show that the new estimators clearly outperform legacy state-of-the-art joint quantile tracking algorithm and achieve faster adaptivity in dynamically varying data streams., Comment: arXiv admin note: text overlap with arXiv:1901.04681

Variance parameters in additive models are typically assigned independent priors that do not account for model structure. We present a new framework for prior selection based on a hierarchical decomposition of the total variance along a tree structure to the individual model components. For each split in the tree, an analyst may be ignorant or have a sound intuition on how to attribute variance to the branches. In the former case a Dirichlet prior is appropriate to use, while in the latter case a penalised complexity (PC) prior provides robust shrinkage. A bottom-up combination of the conditional priors results in a proper joint prior. We suggest default values for the hyperparameters and offer intuitive statements for eliciting the hyperparameters based on expert knowledge. The prior framework is applicable for R packages for Bayesian inference such as INLA and RStan. Three simulation studies show that, in terms of the application-specific measures of interest, PC priors improve inference over Dirichlet priors when used to penalise different levels of complexity in splits. However, when expressing ignorance in a split, Dirichlet priors perform equally well and are preferred for their simplicity. We find that assigning current state-of-the-art default priors for each variance parameter individually is less transparent and does not perform better than using the proposed joint priors. We demonstrate practical use of the new framework by analysing spatial heterogeneity in neonatal mortality in Kenya in 2010-2014 based on complex survey data., Comment: 64 pages, 33 figures

The Exponentially Weighted Average (EWA) of observations is known to be state-of-art estimator for tracking expectations of dynamically varying data stream distributions. However, how to devise an EWA estimator to rather track quantiles of data stream distributions is not obvious. In this paper, we present a lightweight quantile estimator using a generalized form of the EWA. To the best of our knowledge, this work represents the first reported quantile estimator of this form in the literature. An appealing property of the estimator is that the update step size is adjusted online proportionally to the difference between current observation and the current quantile estimate. Thus, if the estimator is off-track compared to the data stream, large steps will be taken to promptly get the estimator back on-track. The convergence of the estimator to the true quantile is proven using the theory of stochastic learning. Extensive experimental results using both synthetic and real-life data show that our estimator clearly outperforms legacy state-of-the-art quantile tracking estimators and achieves faster adaptivity in dynamic environments. The quantile estimator was further tested on real-life data where the objective is efficient online control of indoor climate. We show that the estimator can be incorporated into a concept drift detector for efficiently decide when a machine learning model used to predict future indoor temperature should be retrained/updated.

In this paper, we extend and analyze a Bayesian hierarchical spatio-temporal model for physical systems. A novelty is to model the discrepancy between the output of a computer simulator for a physical process and the actual process values with a multivariate random walk. For computational efficiency, linear algebra for bandwidth limited matrices is utilized, and first-order emulator inference allows for the fast emulation of a numerical partial differential equation (PDE) solver. A test scenario from a physical system motivated by glaciology is used to examine the speed and accuracy of the computational methods used, in addition to the viability of modeling assumptions. We conclude by discussing how the model and associated methodology can be applied in other physical contexts besides glaciology., Comment: Revision accepted for publication by the Journal of Agricultural, Biological, and Environmental Statistics

Renewable sources of energy such as wind power have become a sustainable alternative to fossil fuel-based energy. However, the uncertainty and fluctuation of the wind speed derived from its intermittent nature bring a great threat to the wind power production stability, and to the wind turbines themselves. Lately, much work has been done on developing models to forecast average wind speed values, yet surprisingly little has focused on proposing models to accurately forecast extreme wind speeds, which can damage the turbines. In this work, we develop a flexible spliced Gamma-Generalized Pareto model to forecast extreme and non-extreme wind speeds simultaneously. Our model belongs to the class of latent Gaussian models, for which inference is conveniently performed based on the integrated nested Laplace approximation method. Considering a flexible additive regression structure, we propose two models for the latent linear predictor to capture the spatio-temporal dynamics of wind speeds. Our models are fast to fit and can describe both the bulk and the tail of the wind speed distribution while producing short-term extreme and non-extreme wind speed probabilistic forecasts., Comment: 25 pages

The main goal of disease mapping is to estimate disease risk and identify high-risk areas. Such analyses are hampered by the limited geographical resolution of the available data. Typically the available data are counts per spatial unit and the common approach is the Besag--York--Molli{\'e} (BYM) model. When precise geocodes are available, it is more natural to use Log-Gaussian Cox processes (LGCPs). In a simulation study mimicking childhood leukaemia incidence using actual residential locations of all children in the canton of Z\"urich, Switzerland, we compare the ability of these models to recover risk surfaces and identify high-risk areas. We then apply both approaches to actual data on childhood leukaemia incidence in the canton of Z\"urich during 1985-2015. We found that LGCPs outperform BYM models in almost all scenarios considered. Our findings suggest that there are important gains to be made from the use of LGCPs in spatial epidemiology., Comment: 28 pages, 4 figures, 2 Tables

Varying coefficient models are useful in applications where the effect of the covariate might depend on some other covariate such as time or location. Various applications of these models often give rise to case-specific prior distributions for the parameter(s) describing how much the coefficients vary. In this work, we introduce a unified view of varying coefficients models, arguing for a way of specifying these prior distributions that are coherent across various applications, avoid overfitting and have a coherent interpretation. We do this by considering varying coefficients models as a flexible extension of the natural simpler model and capitalising on the recently proposed framework of penalized complexity (PC) priors. We illustrate our approach in two spatial examples where varying coefficient models are relevant.