Search

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 model. The variance partitioning 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 an intuitive way. We illustrate the advantages of the variance partitioning model using two case studies.

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., 18 pages, 7 figures

Integrated Nested Laplace Approximations (INLA) has been a successful approximate Bayesian inference framework since its proposal by Rue et al. (2009). The increased computational efficiency and accuracy when compared with sampling-based methods for Bayesian inference like MCMC methods, are some contributors to its success. Ongoing research in the INLA methodology and implementation thereof in the R package R-INLA, ensures continued relevance for practitioners and improved performance and applicability of INLA. The era of big data and some recent research developments, presents an opportunity to reformulate some aspects of the classic INLA formulation, to achieve even faster inference, improved numerical stability and scalability. The improvement is especially noticeable for data-rich models. We demonstrate the efficiency gains with various examples of data-rich models, like Cox's proportional hazards model, an item-response theory model, a spatial model including prediction, and a 3-dimensional model for fMRI data.

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

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. Supplementary materials accompanying this paper appear online.

In this paper, we extend and analyze a Bayesian hierarchical spatiotemporal 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.

In this paper we consider the problem of tracking multiple quantiles of dynamicallyvarying data stream distributions. The method is based on making incremental updates ofthe quantile estimates every time a new sample is received. The method is memory andcomputationally efficient since it only stores one value for each quantile estimate and onlyperforms one operation per quantile estimate when a new sample is received from the datastream. The estimates are realistic in the sense that the monotone property of quantiles issatisfied in every iteration. Experiments show that the method efficiently tracks multiplequantiles and outperforms state of the art methods.

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 approach, both implemented in . Heuristics for improving the numerical stability of 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. Supplementary materials accompanying this paper appear on-line.

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 article 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 article 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 Poisson mixture, a zero-inflated Poisson model and a spatial autoregressive combined model. The applications show that the AMIS-INLA approach, in general, outperforms the other methods compared, but the IS-INLA algorithm could be considered for faster inference when good proposals are available. Supplementary materials for this article are available online.

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

Reliable quantification of the global mean surface temperature (GMST) response to radiative forcing is essential for assessing the risk of dangerous anthropogenic climate change. We present the statistical foundations for an observation-based approach using a stochastic linear response model that is consistent with the long-range temporal dependence observed in global temperature variability. We have incorporated the model in a latent Gaussian modeling framework, which allows for the use of integrated nested Laplace approximations (INLAs) to perform full Bayesian analysis. As examples of applications, we estimate the GMST response to forcing from historical data and compute temperature trajectories under the Representative Concentration Pathways (RCPs) for future greenhouse gas forcing. For historic runs in the Model Intercomparison Project Phase 5 (CMIP5) ensemble, we estimate response functions and demonstrate that one can infer the transient climate response (TCR) from the instrumental temperature record. We illustrate the effect of long-range dependence by comparing the results with those obtained from one-box and two-box energy balance models. The software developed to perform the given analyses is publicly available as the R package INLA.climate.

The use of flat or weakly informative priors is popular due to the objective a priori belief in the absence of strong prior information. In the case of the Weibull model the improper uniform, equal parameter gamma and joint Jeffrey's priors for the shape parameter are popular choices. The effects and behaviors of these priors have yet to be established from a modeling viewpoint, especially their ability to reduce to the simpler exponential model. In this work we propose a new principled prior for the shape parameter of the Weibull model, originating from a prior on the distance function, and advocate this new prior as a principled choice in the absence of strong prior information. This new prior can then be used in models with a Weibull modeling component, like competing risks, joint and spatial models, to mention a few. This prior is available in the R-INLA for use, and is applied in a joint longitudinal-survival model framework using the INLA method.

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.

Fractional Gaussian noise (fGn) is a stationary time series model with long-memory properties applied in various fields like econometrics, hydrology and climatology. The computational cost in fitting an fGn model of length n using a likelihood-based approach is $${{\mathcal {O}}}(n^{2})$$ , exploiting the Toeplitz structure of the covariance matrix. In most realistic cases, we do not observe the fGn process directly but only through indirect Gaussian observations, so the Toeplitz structure is easily lost and the computational cost increases to $${{\mathcal {O}}}(n^{3})$$ . This paper presents an approximate fGn model of $${{\mathcal {O}}}(n)$$ computational cost, both with direct and indirect Gaussian observations, with or without conditioning. This is achieved by approximating fGn with a weighted sum of independent first-order autoregressive (AR) processes, fitting the parameters of the approximation to match the autocorrelation function of the fGn model. The resulting approximation is stationary despite being Markov and gives a remarkably accurate fit using only four AR components. Specifically, the given approximate fGn model is incorporated within the class of latent Gaussian models in which Bayesian inference is obtained using the methodology of integrated nested Laplace approximation. The performance of the approximate fGn model is demonstrated in simulations and two real data examples.

Models of excess mortality with random effects were used to estimate regional variation in relative or net survival of cancer patients. Statistical inference for these models based on the Markov chain Monte Carlo (MCMC) methods is computationally intensive and, therefore, not feasible for routine analyses of cancer register data. This study assessed the performance of the integrated nested Laplace approximation (INLA) in monitoring regional variation in cancer survival. Poisson regression model of excess mortality including both spatially correlated and unstructured random effects was fitted to the data of patients diagnosed with ovarian and breast cancer in Finland during 1955–2014 with follow up from 1960 through 2014 by using the period approach with five‐year calendar time windows. We estimated standard deviations associated with variation (i) between hospital districts and (ii) between municipalities within hospital districts. Posterior estimates based on the INLA approach were compared to those based on the MCMC simulation. The estimates of the variation parameters were similar between the two approaches. Variation within hospital districts dominated in the total variation between municipalities. In 2000–2014, the proportion of the average variation within hospital districts was 68% (95% posterior interval: 35–93%) and 82% (60–98%) out of the total variation in ovarian and breast cancer, respectively. In the estimation of regional variation, the INLA approach was accurate, fast, and easy to implement by using the R‐INLA package.

In a bivariate meta-analysis, the number of diagnostic studies involved is often very low so that frequentist methods may result in problems. Using Bayesian inference is particularly attractive as informative priors that add a small amount of information can stabilise the analysis without overwhelming the data. However, Bayesian analysis is often computationally demanding and the selection of the prior for the covariance matrix of the bivariate structure is crucial with little data. The integrated nested Laplace approximations method provides an efficient solution to the computational issues by avoiding any sampling, but the important question of priors remain. We explore the penalised complexity (PC) prior framework for specifying informative priors for the variance parameters and the correlation parameter. PC priors facilitate model interpretation and hyperparameter specification as expert knowledge can be incorporated intuitively. We conduct a simulation study to compare the properties and behaviour of differently defined PC priors to currently used priors in the field. The simulation study shows that the PC prior seems beneficial for the variance parameters. The use of PC priors for the correlation parameter results in more precise estimates when specified in a sensible neighbourhood around the truth. To investigate the usage of PC priors in practice, we reanalyse a meta-analysis using the telomerase marker for the diagnosis of bladder cancer and compare the results with those obtained by other commonly used modelling approaches. Copyright © 2017 John Wiley & Sons, Ltd.

We acknowledge the principal investigators who were responsible for collecting and analysing the soil maps (Jim Dallin, Robert John, Kyle Harms, Robert Stallard and Joe Yavitt), the funding sources (National Science Foundation grants DEB021104, 021115, 0212284 and 0212818 and Office of International Science and Engineering grant 0314581, the Smithsonian Tropical Research Institute soils initiative and the Center for Tropical Forest Science) and field assistants (Paolo Segre and Juan Di Trani). Hyperprior specifications for random fields in spatial point process modelling can have a major impact on the results. In fitting log-Gaussian Cox processes to rainforest tree species, we consider a reparameterised model combining a spatially structured and an unstructured random field into a single component. This component has one hyperpa- rameter accounting for marginal variance, while an additional hyperparameter governs the fraction of the variance explained by the spatially structured effect. This facilitates inter- pretation of the hyperparameters and significance of covariates is studied for a range of hyperprior specifications. Appropriate scaling makes the analysis invariant to grid resolution. Postprint

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 Portuguese National Statistical Institute is responsible for estimating and publishing quarterly labour market figures at national level for both NUTS I and NUTS II regions. Over recent years it has become increasinglyimportant to identify these figures at more disaggregated levels. However, based on the established direct estimation method, it is not possible to produce satisfactorily precise estimates at higher spatial resolutions. From the 4th quarter of 2014 onwards, all the sampling units, namely the residential buildings, of the Portuguese Labour Force Survey (PLFS) were georeferenced. To take full advantage of this information, Pereira et al. (2019) proposed applying a spatial marked point process approach to unemployment estimation, in which the estimation of the unemployment intensity becomes the focal point. There, the sampling units were assumed to be a realization of a spatial point process, specifically a log Gaussian Cox process, with the number of unemployed in each unit being their respective marks. Recently, further information on the geo-referenced locations of all units of the population, namely all residential buildings in the national territory, became available. Consequently, it is no longer necessary to model the spatial configuration of the units of the population. Thus, we propose a new point referenced model for the marks based on the sampled units and extrapolate this to all units of the population. As expected, this extra information, and as a consequence the new model itself, produce estimates with higher precision.

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érn Gaussian random fields. In this review, we discuss the large success of spatial modeling 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 nonstationary spatial models and nonseparable space–time models. Locked until 5.7.2019 due to copyright restrictions. This is the peer reviewed version of an article, which has been published in final form at [https://doi.org/10.1002/wics.1443]. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Self-Archiving.

Source at https://doi.org/10.1214/17-sts576rej. INTRODUCTION: The point of departure for our paper is that most modern statistical models are built to be flexible enough to model diverse data generating mechanisms. Good statistical practice requires us to limit this flexibility, which is typically controlled by a small number of parameters, to the amount “needed” to model the data at hand. The Bayesian framework provides a natural method for doing this although, as DD points out, this trend for penalising model complexity casts a broad shadow over all of modern statistics and data science.The PC prior framework argues for setting priors on these flexibility parameters that are specifically built to penalise a certain type of complexity and avoid overfitting. The discussants raised various points about this core idea. First, DD pointed out that while over-fitting a model is a bad thing, under-fitting is not better: we do not want Occam’s razor to slit our throat. We saw this behaviour when using a half-Normal prior on the distance, while the exponential prior does not lead to obvious attenuation of the estimates. This is confirmed experimentally by Klein and Kneib (2016).

The Northern Humboldt Current System (NHCS) is the world's most productive ecosystem in terms of fish. In particular, the Peruvian anchovy (Engraulis ringens) is the major prey of the main top predators, like seabirds, fish, humans, and other mammals. In this context, it is important to understand the dynamics of the anchovy distribution to preserve it as well as to exploit its economic capacities. Using the data collected by the "Instituto del Mar del Peru" (IMARPE) during a scientific survey in 2005, we present a statistical analysis that has as main goals: (i) to adapt to the characteristics of the sampled data, such as spatial dependence, high proportions of zeros and big size of samples; (ii) to provide important insights on the dynamics of the anchovy population; and (iii) to propose a model for estimation and prediction of anchovy biomass in the NHCS offshore from Peru. These data were analyzed in a Bayesian framework using the integrated nested Laplace approximation (INLA) method. Further, to select the best model and to study the predictive power of each model, we performed model comparisons and predictive checks, respectively. Finally, we carried out a Bayesian spatial influence diagnostic for the preferred model.

This work extends the integrated nested Laplace approximation (INLA) method to latent models outside the scope of latent Gaussian models, where independent components of the latent field can have a near-Gaussian distribution. The proposed methodology is an essential component of a bigger project that aims to extend the R package INLA in order to allow the user to add flexibility and challenge the Gaussian assumptions of some of the model components in a straightforward and intuitive way. Our approach is applied to two examples, and the results are compared with that obtained by Markov chain Monte Carlo, showing similar accuracy with only a small fraction of computational time. Implementation of the proposed extension is available in the R-INLA package.

Summary Motivated by the problem of cancer risk assessment near a nuclear power generating station, the paper describes a methodology for fitting a spatially correlated survival model to large retrospective cohort data sets. Retrospective cohorts, which can be assembled inexpensively from population-based health databases, can partially account for lags between exposures and outcome of chronic diseases such as cancer. These data sets overcome one of the principal limitations of cross-sectional spatial analyses, though performing statistical inference requires accommodating censored and truncated event times as well as spatial dependence. The use of spatial survival models for large retrospective cohorts is described, and Bayesian inference using Markov random-field approximations and integrated nested Laplace approximations is presented. The method is applied to data from individuals living near Pickering Nuclear Generating Station in Canada, showing that the effect of ambient radiation on cancer is not statistically significant.

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., 28 pages, 7 figures

In recent years, disease mapping studies have become a routine application within geographical epidemiology and are typically analysed within a Bayesian hierarchical model formulation. A variety of model formulations for the latent level have been proposed but all come with inherent issues. In the classical BYM model, the spatially structured component cannot be seen independently from the unstructured component. This makes prior definitions for the hyperparameters of the two random effects challenging. There are alternative model formulations that address this confounding, however, the issue on how to choose interpretable hyperpriors is still unsolved. Here, we discuss a recently proposed parameterisation of the BYM model that leads to improved parameter control as the hyperparameters can be seen independently from each other. Furthermore, the need for a scaled spatial component is addressed, which facilitates assignment of interpretable hyperpriors and make these transferable between spatial applications with different graph structures. We provide implementation details for the new model formulation which preserve sparsity properties, and we investigate systematically the model performance and compare it to existing parameterisations. Through a simulation study, we show that the new model performs well, both showing good learning abilities and good shrinkage behaviour. In terms of model choice criteria, the proposed model performs at least equally well as existing parameterisations, but only the new formulation offers parameters that are interpretable and hyperpriors that have a clear meaning., Comment: 24 pages, 6 figures, 5 tables