Search

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

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.

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.

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.

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.

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

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

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

Summary Clinical data on the location of residence at the time of diagnosis of new lupus cases in Toronto, Canada, for the 40 years to 2007 are modelled with the aim of finding areas of abnormally high risk. Inference is complicated by numerous irregular changes in the census regions on which population is reported. A model is introduced consisting of a continuous random spatial surface and fixed effects for time and ages of individuals. The process is modelled on a fine grid and Bayesian inference performed by using integrated nested Laplace approximations. Predicted risk surfaces and posterior probabilities of exceedance are produced for lupus and, for comparison, psoriatic arthritis data from the same clinic. Simulations studies are also carried out to understand better the performance of the model proposed as well as to compare with existing methods.

Summary Suicide has become one of the leading causes of death of Swiss males aged between 15 and 44 years, whose age-standardized rates are about three times higher than those for females. We compared age-stratified suicide rates of Swiss men and women aged 15–79 years and analysed gender-specific differences from 1950 to 2007. Furthermore, we explored whether changes in measures of family integration can explain changes in suicide trends. The use of multivariate age–period–cohort models avoids age aggregation and allows the exploration of heterogeneous time trends across age, period and birth cohort. In addition, explanatory variables can be included. We found strong gender-specific differences in suicide mortality. Although the same risk factors may act on age and overdispersion, there was no significant correlation between gender-specific cohort effects. Family integration had an effect on Swiss suicide risk but only partially explained the underlying trends over time.

Deterministic Bayesian inference for latent Gaussian models has recently become available using integrated nested Laplace approximations (INLA). Applying the INLA-methodology, marginal estimates for elements of the latent field can be computed efficiently, providing relevant summary statistics like posterior means, variances and pointwise credible intervals. In this article, we extend the use of INLA to joint inference and present an algorithm to derive analytical simultaneous credible bands for subsets of the latent field. The algorithm is based on approximating the joint distribution of the subsets by multivariate Gaussian mixtures. Additionally, we present a saddlepoint approximation to compute Bayesian contour probabilities, representing the posterior support of fixed parameter vectors of interest. We perform a simulation study and apply the given methods to two real examples.

The data cloning method is a new computational tool for computing maximum likelihood estimates in complex statistical models such as mixed models. This method is synthesized with integrated nested Laplace approximation to compute maximum likelihood estimates efficiently via a fast implementation in generalized linear mixed models. Asymptotic behavior of the hybrid data cloning method is discussed. The performance of the proposed method is illustrated through a simulation study and real examples. It is shown that the proposed method performs well and rightly justifies the theory. Supplemental materials for this article are available online.

The article considers Bayesian analysis of hierarchical models for count, binomial and multinomial data using efficient MCMC sampling procedures. To this end, an improved method of auxiliary mixture sampling is proposed. In contrast to previously proposed samplers the method uses a bounded number of latent variables per observation, independent of the intensity of the underlying Poisson process in the case of count data, or of the number of experiments in the case of binomial and multinomial data. The bounded number of latent variables results in a more general error distribution, which is a negative log-Gamma distribution with arbitrary integer shape parameter. The required approximations of these distributions by Gaussian mixtures have been computed. Overall, the improvement leads to a substantial increase in efficiency of auxiliary mixture sampling for highly structured models. The method is illustrated for finite mixtures of generalized linear models and an epidemiological case study.

The second-order random walk (RW2) model is commonly used for smoothing data and for modelling response functions. It is computationally efficient due to the Markov properties of the joint (intrinsic) Gaussian density. For evenly spaced locations the RW2 model is well established, whereas for irregularly spaced locations there is no well established construction in the literature. By considering the RW2 model as the solution of a stochastic differential equation (SDE), a dis- cretely observed integrated Wiener process, it is possible to derive the density preserving the Markov properties by augmenting the state-space with the velocities. Here, we derive a computationally more efficient RW2 model for irregular locations using a Galerkin approximation to the solution of the SDE without the need of augmenting the state-space. Numerical comparison with the exact solution demonstrates that the error in the Galerkin approximation is small and negligible in applications.

Many commonly used models in statistics can be formulated as (Bayesian) hierarchical Gaussian Markov random field models. These are characterised by assuming a (often large) Gaussian Markov random field (GMRF) as the second stage in the hierarchical structure and a few hyperparameters at the third stage. Markov chain Monte Carlo is the common approach for Bayesian inference in such models. The variance of the Monte Carlo estimates is Op(M 1/2 ) where M is the number of samples in the chain so, in order to obtain precise estimates of marginal densities, say, we need M to be very large. Inspired by the fact that often one-block and independence samplers can be constructed for hierarchical GMRF models, we will in this work investigate whether MCMC is really needed to estimate marginal densities, which often is the goal of the analysis. By making use of GMRFapproximations, we show by typical examples that marginal densities can indeed be very precisely estimated by deterministic schemes. The methodological and practical consequence of these findings are indeed positive. We conjecture that for many hierarchical GMRF-models there is really no need for MCMC based inference to estimate marginal densities. Further, by making use of numerical methods for sparse matrices the computational costs of these deterministic schemes are nearly instant compared to the MCMC alternative. In particular, we discuss in detail the issue of computing marginal variances for GMRFs.

We illustrate how the recursive algorithm of Reeves & Pettitt (2004) for general factorizable models can be extended to allow exact sampling, maximization of distributions and computation of marginal distributions. All of the methods we describe apply to discrete-valued Markov random fields with nearest neighbour integrations defined on regular lattices; in particular we illustrate that exact inference can be performed for hidden autologistic models defined on moderately sized lattices. In this context we offer an extension of this methodology which allows approximate inference to be carried out for larger lattices without resorting to simulation techniques such as Markov chain Monte Carlo. In particular our work offers the basis for an automatic inference machine for such models. Copyright 2007, Oxford University Press.

Bayesian penalized splines (P-splines) assume an intrinsic Gaussian Markov random field prior on the spline coefficients, conditional on a precision hyper-parameter [Formula: see text]. Prior elicitation of [Formula: see text] is difficult. To overcome this issue, we aim to building priors on an interpretable property of the model, indicating the complexity of the smooth function to be estimated. Following this idea, we propose penalized complexity (PC) priors for the number of effective degrees of freedom. We present the general ideas behind the construction of these new PC priors, describe their properties and show how to implement them in P-splines for Gaussian data.

The principles behind the interface to continuous domain spatial models in the R- INLA software package for R are described. The integrated nested Laplace approximation (INLA) approach proposed by Rue, Martino, and Chopin (2009) is a computationally effective alternative to MCMC for Bayesian inference. INLA is designed for latent Gaussian models, a very wide and flexible class of models ranging from (generalized) linear mixed to spatial and spatio-temporal models. Combined with the stochastic partial differential equation approach (SPDE, Lindgren, Rue, and Lindström 2011), one can accommodate all kinds of geographically referenced data, including areal and geostatistical ones, as well as spatial point process data. The implementation interface covers stationary spatial mod- els, non-stationary spatial models, and also spatio-temporal models, and is applicable in epidemiology, ecology, environmental risk assessment, as well as general geostatistics.

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)

Weintroduceanewcopula-basedcorrectionforgeneralizedlinearmixedmodels(GLMMs)within the integrated nested Laplace approximation (INLA) approach for approximateBayesian inference for latent Gaussian models. GLMMs for Binomial and Poisson datawith many zeroes and low counts have been a somewhat difficult case for INLA, and thecase of binary data has been particularly problematic. Our new correction has been im-plemented as part of the R-INLA package, and adds only negligible computational cost.Empiricalevaluationsonbothrealandsimulateddataindicatethatthemethodworkswell.Keywords: Bayesian computation; copulas; generalized linear mixed models; integratednestedLaplaceapproximation;latentGaussianmodels 1 Introduction IntegratedNestedLaplaceApproximations(INLA)wereintroducedbyRue,MartinoCMcCulloch,SearleN2;:::;n,letProb(y

http://www.jstatsoft.org/v63/i20 The integrated nested Laplace approximation (INLA) provides an interesting way of approximating the posterior marginals of a wide range of Bayesian hierarchical models. This approximation is based on conducting a Laplace approximation of certain functions and numerical integration is extensively used to integrate some of the models parameters out. The R - INLA package o ers an interface to INLA, providing a suitable framework for data analysis. Although the INLA methodology can deal with a large number of models, only the most relevant have been implemented within R - INLA . However, many other important models are not available for R - INLA yet. In this paper we show how to fit a number of spatial models with R - INLA , including its interaction with other R packages for data analysis. Secondly, we describe a novel method to extend the number of latent models available for the model parameters. Our approach is based on conditioning on one or several model parameters and fit these conditioned models with R-INLA. Then these models are combined using Bayesian model averaging to provide the final approximations to the posterior marginals of the model. Finally, we show some examples of the application of this technique in spatial statistics. It is worth noting that our approach can be extended to a number of other fields, and not only spatial statistics

Summary To account for measurement error (ME) in explanatory variables, Bayesian approaches provide a flexible framework, as expert knowledge can be incorporated in the prior distributions. Recently, integrated nested Laplace approximations have been proven to be a computationally convenient alternative to sampling approaches for Bayesian inference in latent Gaussian models. We show how the most common approaches to adjust for ME, the classical and the Berkson ME, fit into this framework. This is achieved through a reformulation with augmented pseudo-observations and a suitable extension of the latent Gaussian field. Two specific classes are described, which allow for a particularly simple implementation using integrated nested Laplace approximations. We present three applications within the framework of generalized linear (mixed) models with ME. To illustrate the practical feasibility, R code is provided in on-line supplementary material.

Local linear M-smoothing is proposed as a method for noise reduction in one-dimensional signals. It is more appropriate than conventional local linear smoothing, because it does not introduce blurring of jumps in the signal. It improves local constant M-smoothing, by avoiding boundary effects at edges and jumps. While the idea of local linear M-smoothing is straightforward, numerical issues are challenging, because of the local minima aspect that is crucial to good performance. We give an algorithm which is both fast and robust together with the theoretical properties of the local linear M-smoother. The new M-smoother gives a large improvement for some data sets compared to the local constant M-smoother and demonstrates elsewhere good performance on various artificial and magnetic resonance data sets.

In this paper, we introduce a new concept for constructing prior distributions. We exploit the natural nested structure inherent to many model components, which defines the model component to be a flexible extension of a base model. Proper priors are defined to penalise the complexity induced by deviating from the simpler base model and are formulated after the input of a user-defined scaling parameter for that model component, both in the univariate and the multivariate case. These priors are invariant to reparameterisations, have a natural connection to Jeffreys' priors, are designed to support Occam's razor and seem to have excellent robustness properties, all which are highly desirable and allow us to use this approach to define default prior distributions. Through examples and theoretical results, we demonstrate the appropriateness of this approach and how it can be applied in various situations., Major revision of previous version. Includes a beefed up literature review and new desiderata for hierarchical priors. Removes (for space) the Cox proportional hazard model and the section on hyperparameters for Gaussian random fields

In recent years, a number of statistical models have been proposed for the purposes of high-level image analysis tasks such as object recognition. However, in general, these models remain hard to use in practice, partly as a result of their complexity, partly through lack of software. In this paper we concentrate on a particular deformable template model which has proved potentially useful for locating and labelling cells in microscope slides Rue and Hurn (1999). This model requires the specification of a number of rather non-intuitive parameters which control the shape variability of the deformed templates. Our goal is to arrange the estimation of these parameters in such a way that the microscope user's expertise is exploited to provide the necessary training data graphically by identifying a number of cells displayed on a computer screen, but that no additional statistical input is required. In this paper we use maximum likelihood estimation incorporating the error structure in the generation of our training data.

This paper addresses the task of locating and identifying an unknown number of objects of different types in an image. Baddeley & Van Lieshout (1993) advocate marked point processes as object priors, whereas Grenander & Miller (1994) use deformable template models. In this paper elements of both approaches are combined to handle scenes containing variable numbers of objects of different types, using reversible jump Markov chain Monte Carlo methods for inference (Green, 1995). The naive application of these methods here leads to slow mixing and we adapt the model and algorithm in tandem in proposing three strategies to deal with this. The first two expand the model space by introducing an additional 'unknown' object type and the idea of a variable resolution template. The third strategy, utilising the first two, augments the algorithm with classes of updates which provide intuitive transitions between realisations containing different numbers of cells by splitting or merging nearby objects.

Markov chain Monte Carlo methods are widely used to study highly structured stochastic systems. However, when the system is subject to constraints, it is difficult to find irreducible proposal distributions. We suggest a “block-wise” approach for constrained sampling and optimisation.

In this article we adopt Baddeley's delta metric as a loss function in Bayesian image restoration and classification. We develop a new algorithm that allows us to approximate the corresponding optimal Bayesian estimator. With this algorithm good practical estimates can be obtained at approximately the same computational cost as traditional estimators like the marginal posterior mode (MPM). A comparison of our proposed classification with MPM shows significant advantages, especially with respect to fine structures.

Common loss functions used for the restoration of grey scale images include the zero-one loss and the sum of squared errors. The corresponding estimators, the posterior mode and the posterior marginal mean, are optimal Bayes estimators with respect to their way of measuring the loss for different error configurations. However, both these loss functions have a fundamental weakness: the loss does not depend on the spatial structure of the errors. This is important because a systematic structure in the errors can lead to misinterpretation of the estimated image. We propose a new loss function that also penalizes strong local sample covariance in the error and we discuss how the optimal Bayes estimator can be estimated using a two-step Markov chain Monte Carlo and simulated annealing algorithm. We present simulation results for some artificial data which show improvement with respect to small structures in the image.

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.

This paper proposes 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 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.

This paper develops methodology that provides a toolbox for routinely fitting complex models to realistic spatial point pattern data. We consider models that are based on log-Gaussian Cox processes and include local interaction in these by considering constructed covariates. This enables us to use integrated nested Laplace approximation and to considerably speed up the inferential task. In addition, methods for model comparison and model assessment facilitate the modelling process. The performance of the approach is assessed in a simulation study. To demonstrate the versatility of the approach, models are fitted to two rather different examples, a large rainforest data set with covariates and a point pattern with multiple marks., Comment: Published in at http://dx.doi.org/10.1214/11-AOAS530 the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)

Unlike the development of more accurate prior distributions for use in Bayesian imaging, the design of more sensible estimators through loss functions has been neglected in the literature. We discuss the design of loss functions with a local structure that depend only on a binary misclassification vector. The proposed approach is similar to modeling with a Markov random field. The Bayes estimate is calculated in a two-step algorithm using Markov chain Monte Carlo and simulated annealing algorithms. We present simulation experiments with the Ising model, where the observations are corrupted with Gaussian and flip noise.

Next generation sequencing is quickly replacing microarrays as a technique to probe different molecular levels of the cell, such as DNA or RNA. The technology provides higher resolution, while reducing bias. RNA sequencing results in counts of RNA strands. This type of data imposes new statistical challenges. We present a novel, generic approach to model and analyze such data. Our approach aims at large flexibility of the likelihood (count) model and the regression model alike. Hence, a variety of count models is supported, such as the popular NB model, which accounts for overdispersion. In addition, complex, non-balanced designs and random effects are accommodated. Like some other methods, our method provides shrinkage of dispersion-related parameters. However, we extend it by enabling joint shrinkage of parameters, including those for which inference is desired. We argue that this is essential for Bayesian multiplicity correction. Shrinkage is effectuated by empirically estimating priors. We discuss several parametric (mixture) and non-parametric priors and develop procedures to estimate (parameters of) those. Inference is provided by means of local and Bayesian false discovery rates. We illustrate our method on several simulations and two data sets, also to compare it with other methods. Model-and data-based simulations show substantial improvements in the sensitivity at the given specificity. The data motivate the use of the ZI-NB as a powerful alternative to the NB, which results in higher detection rates for low-count data. Finally, compared with other methods, the results on small sample subsets are more reproducible when validated on their large sample complements, illustrating the importance of the type of shrinkage. © 2012 The Author.

To analyze and project age-specific mortality or morbidity rates age-period-cohort (APC) models are very popular. Bayesian approaches facilitate estimation and improve predictions by assigning smoothing priors to age, period and cohort effects. Adjustments for overdispersion are straightforward using additional random effects. When rates are further stratified, for example, by countries, multivariate APC models can be used, where differences of stratum-specific effects are interpretable as log relative risks. Here, we incorporate correlated stratum-specific smoothing priors and correlated overdispersion parameters into the multivariate APC model, and use Markov chain Monte Carlo and integrated nested Laplace approximations for inference. Compared to a model without correlation, the new approach may lead to more precise relative risk estimates, as shown in an application to chronic obstructive pulmonary disease mortality in three regions of England and Wales. Furthermore, the imputation of missing data for one particular stratum may be improved, since the new approach takes advantage of the remaining strata if the corresponding observations are available there. This is shown in an application to female mortality in Denmark, Sweden and Norway from the 20th century, where we treat for each country in turn either the first or second half of the observations as missing and then impute the omitted data. The projections are compared to those obtained from a univariate APC model and an extended Lee--Carter demographic forecasting approach using the proper Dawid--Sebastiani scoring rule., Comment: Published in at http://dx.doi.org/10.1214/11-AOAS498 the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)

Summary Continuously indexed Gaussian fields (GFs) are the most important ingredient in spatial statistical modelling and geostatistics. The specification through the covariance function gives an intuitive interpretation of the field properties. On the computational side, GFs are hampered with the big n problem, since the cost of factorizing dense matrices is cubic in the dimension. Although computational power today is at an all time high, this fact seems still to be a computational bottleneck in many applications. Along with GFs, there is the class of Gaussian Markov random fields (GMRFs) which are discretely indexed. The Markov property makes the precision matrix involved sparse, which enables the use of numerical algorithms for sparse matrices, that for fields in ℝ2 only use the square root of the time required by general algorithms. The specification of a GMRF is through its full conditional distributions but its marginal properties are not transparent in such a parameterization. We show that, using an approximate stochastic weak solution to (linear) stochastic partial differential equations, we can, for some GFs in the Matérn class, provide an explicit link, for any triangulation of ℝd, between GFs and GMRFs, formulated as a basis function representation. The consequence is that we can take the best from the two worlds and do the modelling by using GFs but do the computations by using GMRFs. Perhaps more importantly, our approach generalizes to other covariance functions generated by SPDEs, including oscillating and non-stationary GFs, as well as GFs on manifolds. We illustrate our approach by analysing global temperature data with a non-stationary model defined on a sphere.

Bayesian analysis of time-to-event data, usually called survival analysis, has received increasing attention in the last years. In Cox-type models it allows to use information from the full likelihood instead of from a partial likelihood, so that the baseline hazard function and the model parameters can be jointly estimated. In general, Bayesian methods permit a full and exact posterior inference for any parameter or predictive quantity of interest. On the other side, Bayesian inference often relies on Markov chain Monte Carlo (MCMC) techniques which, from the user point of view, may appear slow at delivering answers. In this article, we show how a new inferential tool named integrated nested Laplace approximations can be adapted and applied to many survival models making Bayesian analysis both fast and accurate without having to rely on MCMC-based inference.

Generalized linear mixed models (GLMMs) continue to grow in popularity due to their ability to directly acknowledge multiple levels of dependency and model different data types. For small sample sizes especially, likelihood-based inference can be unreliable with variance components being particularly difficult to estimate. A Bayesian approach is appealing but has been hampered by the lack of a fast implementation, and the difficulty in specifying prior distributions with variance components again being particularly problematic. Here, we briefly review previous approaches to computation in Bayesian implementations of GLMMs and illustrate in detail, the use of integrated nested Laplace approximations in this context. We consider a number of examples, carefully specifying prior distributions on meaningful quantities in each case. The examples cover a wide range of data types including those requiring smoothing over time and a relatively complicated spline model for which we examine our prior specification in terms of the implied degrees of freedom. We conclude that Bayesian inference is now practically feasible for GLMMs and provides an attractive alternative to likelihood-based approaches such as penalized quasi-likelihood. As with likelihood-based approaches, great care is required in the analysis of clustered binary data since approximation strategies may be less accurate for such data.

SummaryStructured additive regression models are perhaps the most commonly used class of models in statistical applications. It includes, among others, (generalized) linear models, (generalized) additive models, smoothing spline models, state space models, semiparametric regression, spatial and spatiotemporal models, log-Gaussian Cox processes and geostatistical and geoadditive models. We consider approximate Bayesian inference in a popular subset of structured additive regression models, latent Gaussian models, where the latent field is Gaussian, controlled by a few hyperparameters and with non-Gaussian response variables. The posterior marginals are not available in closed form owing to the non-Gaussian response variables. For such models, Markov chain Monte Carlo methods can be implemented, but they are not without problems, in terms of both convergence and computational time. In some practical applications, the extent of these problems is such that Markov chain Monte Carlo sampling is simply not an appropriate tool for routine analysis. We show that, by using an integrated nested Laplace approximation and its simplified version, we can directly compute very accurate approximations to the posterior marginals. The main benefit of these approximations is computational: where Markov chain Monte Carlo algorithms need hours or days to run, our approximations provide more precise estimates in seconds or minutes. Another advantage with our approach is its generality, which makes it possible to perform Bayesian analysis in an automatic, streamlined way, and to compute model comparison criteria and various predictive measures so that models can be compared and the model under study can be challenged.

This thesis consists of five papers, presented in chronological order. Their content is summarised in this section.Paper I introduces the approximation tool for latent GMRF models and discusses, in particular, the approximation for the posterior of the hyperparameters θ in equation (1). It is shown that this approximation is indeed very accurate, as even long MCMC runs cannot detect any error in it. A Gaussian approximation to the density of χi|θ, y is also discussed. This appears to give reasonable results and it is very fast to compute. However, slight errors are detected when comparing the approximation with long MCMC runs. These are mostly due to the fact that a possible - skewed density is approximated via a symmetric one. Paper I presents also some details about sparse matrices algorithms.The core of the thesis is presented in Paper II. Here most of the remaining issues present in Paper I are solved. Three different approximation for χi|θ, y with different degrees of accuracy and computational costs are described. Moreover, ways to assess the approximation error and considerations about the asymptotical behaviour of the approximations are also discussed. Through a series of examples covering a wide range of commonly used latent GMRF models, the approximations are shown to give extremely accurate results in a fraction of the computing time used by MCMC algorithms.Paper III applies the same ideas as Paper II to generalised linear mixed models where χ represents a latent variable at n spatial sites on a two dimensional domain. Out of these n sites k, with n >> k , are observed through data. The n sites are assumed to be on a regular grid and wrapped on a torus. For the class of models described in Paper III the computations are based on discrete Fourier transform instead of sparse matrices. Paper III illustrates also how marginal likelihood π (y) can be approximated, provides approximate strategies for Bayesian outlier detection and perform approximate evaluation of spatial experimental design.Paper IV presents yet another application of the ideas in Paper II. Here approximate techniques are used to do inference on multivariate stochastic volatility models, a class of models widely used in financial applications. Paper IV discusses also problems deriving from the increased dimension of the parameter vector θ, a condition which makes all numerical integration more computationally intensive. Different approximations for the posterior marginals of the parameters θ, π(θi)|y), are also introduced. Approximations to the marginal likelihood π(y) are used in order to perform model comparison.Finally, Paper V is a manual for a program, named inla which implements all approximations described in Paper II. A large series of worked out examples, covering many well known models, illustrate the use and the performance of the inla program. This program is a valuable instrument since it makes most of the Bayesian inference techniques described in this thesis easily available for everyone.

In an empirical Bayesian setting, we provide a new multiple testing method, useful when an additional covariate is available, that influences the probability of each null hypothesis being true. We measure the posterior significance of each test conditionally on the covariate and the data, leading to greater power. Using covariate-based prior information in an unsupervised fashion, we produce a list of significant hypotheses which differs in length and order from the list obtained by methods not taking covariate-information into account. Covariate-modulated posterior probabilities of each null hypothesis are estimated using a fast approximate algorithm. The new method is applied to expression quantitative trait loci (eQTL) data., Published in at http://dx.doi.org/10.1214/08-AOAS158 the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org)

The paper ‘Modern statistics for spatial point processes' by Jesper Møller and Rasmus P. Waagepetersen is based on a special invited lecture given by the authors at the 21st Nordic Conference on Mathematical Statistics, held at Rebild, Denmark, in June 2006. At the conference, Antti Penttinen and Eva B. Vedel Jensen were invited to discuss the paper. We here present the comments from the two invited discussants and from a number of other scholars, as well as the authors' responses to these comments. Below Figure 1, Figure 2, etc., refer to figures in the paper under discussion, while Figure A, Figure B, etc., refer to figures in the current discussion. All numbered sections and formulas refer to the paper. The paper ‘Modern statistics for spatial point processes' by Jesper Møller and Rasmus P. Waagepetersen is based on a special invited lecture given by the authors at the 21st Nordic Conference on Mathematical Statistics, held at Rebild, Denmark, in June 2006. At the conference, Antti Penttinen and Eva B. Vedel Jensen were invited to discuss the paper. We here present the comments from the two invited discussants and from a number of other scholars, as well as the authors' responses to these comments. Below Figure 1, Figure 2, etc., refer to figures in the paper under discussion, while Figure A, Figure B, etc., refer to figures in the current discussion. All numbered sections and formulas refer to the paper.

SummaryGaussian Markov random-field (GMRF) models are frequently used in a wide variety of applications. In most cases parts of the GMRF are observed through mutually independent data; hence the full conditional of the GMRF, a hidden GMRF (HGMRF), is of interest. We are concerned with the case where the likelihood is non-Gaussian, leading to non-Gaussian HGMRF models. Several researchers have constructed block sampling Markov chain Monte Carlo schemes based on approximations of the HGMRF by a GMRF, using a second-order expansion of the log-density at or near the mode. This is possible as the GMRF approximation can be sampled exactly with a known normalizing constant. The Markov property of the GMRF approximation yields computational efficiency.The main contribution in the paper is to go beyond the GMRF approximation and to construct a class of non-Gaussian approximations which adapt automatically to the particular HGMRF that is under study. The accuracy can be tuned by intuitive parameters to nearly any precision. These non-Gaussian approximations share the same computational complexity as those which are based on GMRFs and can be sampled exactly with computable normalizing constants. We apply our approximations in spatial disease mapping and model-based geostatistical models with different likelihoods, obtain procedures for block updating and construct Metropolized independence samplers.