{'fraction_non_alphanumeric': 0.045066798265712575, 'fraction_numerical': 0.014223271881278862, 'mean_word_length': 4.537673238336912, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Computing with INLA: A Review

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.


INTRODUCTION
A key obstacle in Bayesian statistics is to actually do the Bayesian inference. From a mathematical point of view, the inference step is easy, transparent and defined by first principles: We simply update prior beliefs about the unknown parameters with available information in observed data, and obtain the posterior distribution for the parameters. Based on the posterior, we can compute relevant statistics for the parameters of interest, including marginal distributions, means, variances, quantiles, credibility intervals, etc. In practice, this is much easier said than done.
The introduction of simulation based inference, through the idea of Markov chain Monte Carlo (Robert and Casella, 1999), hit the statistical community in the early 1990's and represented a major break-through in Bayesian inference. MCMC provided a general recipe to generate samples from posteriors by constructing a Markov chain with the target posterior as the stationary distribution. This made it possible (in theory) to extract and compute whatever one could wish for. Additional major developments have paved the way for popular user-friendly MCMC-tools, like WinBUGS (Spiegelhalter et al., 1995), JAGS (Plummer, 2016), and the new initiative Stan (Stan Development Team, 2015), which uses Hamiltonian Monte Carlo. Armed with these and similar tools, Bayesian statistics has quickly grown in popularity and Bayesian statistics is now well-represented in all the major research journals in all branches of statistics.
In our opinion, however, from the point of view of applied users, the impact of the Bayesian revolution has been less apparent. This is not a statement about how Bayesian statistics itself is viewed by that community, but about its rather "cumbersome" inference, which still requires a lot of CPU -and hence human time-as well as tweaking of simulation and model parameters to get it right. Re-running a lot of alternative models gets even more cumbersome, making the iterative process of model building in statistical analysis impossible (Box and Tiao, 1973, Sec. 1.1.4). For this reason, simulation based inference (and hence in most cases also Bayesian statistics) has too often been avoided as being practically infeasible.
In this paper, we review a different take on doing Bayesian inference that recently has facilitated the uptake of Bayesian modelling within the community of applied users. The given approach is restricted to the specific class of latent Gaussian models (LGMs) which, as will be clear soon, includes a wide variety of commonly applied statistical models making this restriction less limiting than it might appear at first sight. The crucial point here is that we can derive integrated nested Laplace approximation (INLA methodology) for LGMs, a deterministic approach to approximate Bayesian inference. Performing inference within a reasonable time-frame, in most cases INLA is both faster and more accurate than MCMC alternatives. Being used to trading speed for accuracy this might seem like a contradiction to most readers. The corresponding R-package (R-INLA, see www.r-inla.org), has turned out to be very popular in applied sciences and applied statistics, and has become a versatile tool for quick and reliable Bayesian inference.
Recent examples of applications using the R-INLA package for statistical analysis, include disease mapping (Schrödle and Held, 2011b,a;Ugarte et al., 2014Ugarte et al., , 2016Papoila et al., 2014;Goicoa et al., 2016;, age-period-cohort models (Riebler and Held, 2016), evolution of the Ebola virus (Santermans et al., 2016), studies of relationship between access to housing, health and well-being in cities (Kandt et al., 2016), study of the prevalence and correlates of intimate partner violence against men in Africa (Tsiko, 2015), search for evidence of gene expression heterosis (Niemi et al., 2015), analysis of traffic pollution and hospital admissions in London (Halonen et al., 2016), early transcriptome changes in maize primary root tissues in response to moderate water deficit conditions by RNA-Sequencing (Opitz et al., 2016), performance of inbred and hybrid genotypes in plant breeding and genetics (Lithio and Nettleton, 2015), a study of Norwegian emergency wards (Goth et al., 2014), effects of measurement errors (Kröger et al., 2016;Muff and Keller, 2015), network meta-analysis (Sauter and Held, 2015), time-series analysis of genotyped human campylobacteriosis cases from the Manawatu region of New Zealand (Friedrich et al., 2016), modeling of parrotfish habitats (Roos et al., 2015b), Bayesian outbreak detection (Salmon et al., 2015), studies of long-term trends in the number of Monarch butterflies (Crewe and Mccracken, 2015), long-term effects on hospital admission and mortality of road traffic noise (Halonen et al., 2015), spatio-temporal dynamics of brain tumours (Iulian et al., 2015), ovarian cancer mortality (García-Pérez et al., 2015), the effect of preferential sampling on phylodynamic inference (Karcher et al., 2016), analysis of the impact of climate change on abundance trends in central Europe (Bowler et al., 2015), investigation of drinking patterns in US Counties from 2002 to 2012 (Dwyer-Lindgren et al., 2015), resistance and resilience of terrestrial birds in drying climates (Selwood et al., 2015), cluster analysis of population amyotrophic lateral sclerosis risk (Rooney et al., 2015), malaria infection in Africa (Noor et al., 2014), effects of fragmentation on infectious disease dynamics (Jousimo et al., 2014), soil-transmitted helminth infection in sub-Saharan Africa (Karagiannis-Voules et al., 2015), analysis of the effect of malaria control on Plasmodium falciparum in Africa between 2000 and 2015 (Bhatt et al., 2015), adaptive prior weighting in generalized regression , analysis of hand, foot, and mouth disease surveillance data in China (Bauer et al., 2016), estimate the biomass of anchovies in the coast of Perú (Quiroz et al., 2015), and many others.
We review the key components that make up INLA in Section 2 and in Section 3 we combine these to outline why -and in which situations -INLA works. In Section 4 we show some examples of the use of R-INLA, and discuss some special features that expand the class of models that R-INLA can be applied to. In Section 5, we discuss a specific challenge in Bayesian methodology, and, in particular, reason why it is important to provide better suggestions for default priors. We conclude with a general discussion and outlook in Section 6.

BACKGROUND ON THE KEY COMPONENTS
In this section, we review the key components of the INLA-approach to approximate Bayesian inference. We introduce these concepts using a top-down approach, starting with latent Gaussian models (LGMs), and what type of statistical models may be viewed as LGMs. We also discuss the types of Gaussians/Gaussian-processes that are computationally efficient within this formulation, and illustrate Laplace approximation to perform integration -a method that has been around for a very long time yet proves to be a key ingredient in the methodology we review here.
Due to the top-down structure of this text we occasionally have to mention specific concepts before properly introducing and/or defining them -we ask the reader to bear with us in these cases.

Latent Gaussian Models (LGMs)
The concept of latent Gaussian models represents a very useful abstraction subsuming a large class of statistical models, in the sense that the task of statistical inference can be unified for the entire class (Rue et al., 2009). This is obtained using a three-stage hierarchical model formulation, in which observations y can be assumed to be conditionally independent, given a latent Gaussian random field x and hyperparameters θ 1 , The versatility of the model class relates to the specification of the latent Gaussian field: which includes all random terms in a statistical model, describing the underlying dependence structure of the data. The hyperparameters θ = (θ 1 , θ 2 ), control the Gaussian latent field and/or the likelihood for the data, and the posterior reads We make the following critical assumptions : 1. The number of hyperparameters |θ| is small, typically 2 to 5, but not exceeding 20.
2. The distribution of the latent field, x|θ is Gaussian and required to be a Gaussian Markov random field (GMRF) (or do be close to one) when the dimension n is high (10 3 to 10 5 ).
3. The data y are mutually conditionally independent of x and θ, implying that each observation y i only depends on one component of the latent field, e.g. x i . Most components of x will not be observed.
These assumptions are required both for computational reasons and to ensure, with a high degree of certainty, that the approximations we describe below are accurate.

Additive Models
Now, how do LGMs relate to other better-known statistical models? Broadly speaking, they are an umbrella class generalising the large number of related variants of "additive" and/or "generalized" (linear) models. For instance, interpreting the likelihood π(y i |x i , θ), so that "y i only depends on its linear predictor x i ", yields the generalized linear model setup. We can interpret {x i , i ∈ I} as η i (the linear predictor), which itself is additive with respect to other effects, Here, µ is the overall intercept and z are fixed covariates with linear effects {β j }. The difference between this formulation and an ordinary generalized linear model are the terms {f k }, which are used to represent specific Gaussian processes. We label each f k as a model component, in which element j contributes to the ith linear predictor. Examples of model components f k include auto-regressive time-series models, stochastic spline models and models for smoothing, measurement error models, random effects models with different types of correlations, spatial models etc. We assume that the model components are a-priori independent, the fixed effects (µ, β) have a joint Gaussian prior and that the fixed effects are a-priori independent of the model components.
The key is now that the model formulation in (2) and LGMs relate to the same class of models when we assume Gaussian priors for the intercept and the parameters of the fixed effects. The joint distribution of is then Gaussian, and also non-singular if we add a tiny noise term in (2). This yields the latent field x in the hierarchical LGM formulation. Clearly, dim(x) = n can easily get large, as it equals the number of observations, plus the intercept(s) and fixed effects, plus the sum of the dimension of all the model components.
The hyperparameters θ comprise the parameters of the likelihood and the model components. A likelihood family and each model component, typically has between zero and two hyperparameters. These parameters often include some kind of variance, scale or correlation parameters. Nicely, the number of hyperparameters is typically small and further, does not depend on the dimension of the latent field n nor the number of observations. This is crucial for computational efficiency, as even with a big dataset, the number of hyperparameters remains constant and assumption 1. still holds.

Gaussian Markov Random Fields (GMRFs)
In practice, the latent field should not only be Gaussian, but should also be a (sparse) Gaussian Markov random field (GMRF); see Held (2005, 2010);  for an introduction to GMRFs. A GMRF x is simply a Gaussian with additional conditional independence properties, meaning that x i and x j are conditionally independent given the remaining elements x −ij , for quite a few {i, j}'s. The simplest non-trivial example is the first-order auto-regressive model, x t = φx t−1 + t , t = 1, 2, . . . , m, having Gaussian innovations . For this model, the correlation between x t and x s is φ |s−t| and the resulting m × m covariance matrix is dense. However, x s and x t are conditionally independent given x −st , for all |s − t| > 1. In the Gaussian case, a very useful consequence of conditional independence is that this results in zeros for pairs of conditionally independent values in the precision matrix (the inverse of the covariance matrix). Considering GMRFs provides a huge computational benefit, as calculations involving a dense m × m matrix are much more costly than when a sparse matrix is used. In the auto-regressive example, the precision matrix is tridiagonal and can be factorized in O(m) time, whereas we need O(m 3 ) in the general dense case. Memory requirement is also reduced, O(m) compared to O(m 2 ), which makes it much easier to run larger models. For models with a spatial structure, the cost is O(m 3/2 ) paired with a O(m log(m)) memory requirement. In general, the computational cost depends on the actual sparsity pattern in the precision matrix, hence it is hard to provide precise estimates.

Additive Models and GMRFs
In the construction of additive models including GMRFs the following fact provides some of the "magic" that is exploited in INLA: The joint distribution for x in (3) is also a GMRF and its precision matrix consists of sums of the precision matrices of the fixed effects and the other model components.
We will see below that we need to form the joint distribution of the latent field many times, as it depends on the hyperparameters θ. Hence, it is essential that this can be done efficiently avoiding computationally costly matrix operations. Being able to simply treat the joint distribution as a GMRF with a precision matrix that is easy to compute, is one of the key reasons why the INLAapproach is so efficient. Also, the sparse structure of the precision matrix boosts computationally efficiency, compared with operations on dense matrices.
To illustrate more clearly what happens, let us consider the following simple example, where we have added a small amount of noise i . The two model components f 1j 1 (i) and f 2j 2 (i) have sparse precision matrices Q 1 (θ) and Q 2 (θ), of dimension m 1 × m 1 and m 2 × m 2 , respectively. Let τ µ and τ β be the (fixed) prior precisions for µ and β. We can express (4) using matrices, Here, A 1 , and similarly for A 2 , is a n × m 1 sparse matrix, which is zero except for exactly one 1 in each row. The joint precision matrix of (η, f 1 , f 2 , β, µ) is straight forward to obtain by rewriting The dimension is n + m 1 + m 2 + 2. Concretely, the above-mentioned "magic" implies that the only matrices that need to be multiplied are the A-matrices, which are extremely sparse and contain only one non-zero element in each row. These matrix products do not depend on θ and hence they only need to be computed once. The joint precision matrix only depends on θ through Q 1 (θ) and Q 2 (θ) and as θ change, the computational cost of re-computing Q joint (θ) is negligible.
The sparsity of Q joint (θ) illustrates how the additive structure of the model facilitates computational efficiency. For simplicity, assume n = m 1 = m 2 , and denote by e 1 and e 2 the average number of non-zero elements in a row of Q 1 (θ) and Q 2 (θ), respectively. An upper bound for the number of non-zero terms in Q joint (θ) is n(19 + e 1 + e 2 ) + 4. Approximately, this gives on average only (19 + e 1 + e 2 )/3 non-zero elements for a row in Q joint (θ), which is very sparse.

Laplace Approximations
The Laplace approximation or method, is an old technique for the approximation of integrals; see Barndorff-Nielsen and Cox (1989, Ch. 3.3) for a general introduction. The setting is as follows. The aim is to approximate the integral, as n → ∞. Let x 0 be the point in which f (x) has its maximum, then The idea is simple but powerful: Approximate the target with a Gaussian, matching the mode and the curvature at the mode. By interpreting nf (x) as the sum of log-likelihoods and x as the unknown parameter, the Gaussian approximation will be exact as n → ∞, if the central limit theorem holds. The extension to higher dimensional integrals, is immediate and the error turns out to be This is a good result for two reasons. The error is relative and with rate n −1 , as opposed to an additive error and a rate n −1/2 , which are common in simulation-based inference. The Laplace approximation used to be a key tool for doing high-dimensional integration in pre-MCMC times, but quickly went out of fashion when MCMC entered the stage. But how does it relate to what we endeavour to do here? Lets assume that we would like to compute a marginal distribution π(γ 1 ) from a joint distribution π(γ) where we have exploited the fact that we approximate π(γ −1 |γ 1 ) with a Gaussian. In the context of the LGMs we have γ = (x, θ). Tierney and Kadane (1986) show that if π(γ) ∝ exp(nf n (γ)), i.e. if f n (γ) is the average log likelihood, the relative error of the normalized approximation (7), within a O(n −1/2 ) neighbourhood of the mode, is O(n −3/2 ). In other words, if we have n replicated data from the same parameters, γ, we can compute posterior marginals with a relative error of O(n −3/2 ), assuming the numerical error to be negligible. This is an extremely positive result, but unfortunately the underlying assumptions usually do not hold.
1. Instead of replicated data from the same model, we may have one replicate from one model (as is common in spatial statistics), or several observations from similar models.
2. The implicit assumption in the above result is also that |γ| is fixed as n → ∞. However, there is only one realisation for each observation/location in the random effect(s) in the model, implying that |γ| grows with n.
Is it still possible to gain insight into when the Laplace approximation would give good results, even if these assumptions do not hold? First, let's replace replicated observations from the same model, with several observations from similar models -where we deliberately use the term "similar" in a loose sense. We can borrow strength across variables that we a-priori assume to be similar, for example in smoothing over time or over space. In this case, the resulting linear predictors for two observations could differ in only one realisation of the random effect. In addition, borrowing strength and smoothing can reduce the effect of the model dimension growing with n, since the effective dimension can then grow much more slowly with n. Another way to interpret the accuracy in computing posterior marginals using Laplace approximations, is to not look at the error-rate but at the implicit constant upfront. If the posterior is close to a Gaussian density, the results will be more accurate compared to a density that is very different from a Gaussian. This is similar to the convergence for the central limit theorem where convergence is faster if relevant properties such as uni-modality, symmetry and tail behaviour are satisfied; see for example Baghishani and Mohammadzadeh (2012). Similarly, in the context here uni-modality is necessary since we approximate the integrand with a Gaussian. Symmetry helps since the Gaussian distribution is symmetric, while heavier tails will be missed by the Gaussian. For example, assume with centred covariates z. We then expect better accuracy for π(γ 1 ), having high counts compared with low counts. With high counts, the Poisson distribution is approximately Gaussian and almost symmetric. Low counts are more challenging, since the likelihood for y i = 0 and z i = 0, is proportional to exp(− exp(γ 1 )), which has a maximum value at γ 1 = −∞. The situation is similar for binomial data of size m, where low values of m are more challenging than high values of m. Theoretical results for the current rather "vague" context are difficult to obtain and constitute a largely unsolved problem; see for example Shun and McCullagh (1995); Kauermann et al. (2009);Ogden (2016).
Let us now discuss a simplistic, but realistic, model in two dimensions for a constant c > 0 and ρ ≥ 0. This is the same functional form as we get from two Bernoulli successes, using a logit-link. Using the constant c is an alternative to scaling the Gaussian part, and the case where ρ < 0 is similar. The task now is to approximate π(x 1 ) = π(x 1 , x 2 )/π(x 2 |x 1 ), using (7). Here, the Gaussian approximation is indexed by x 1 and we use one Laplace approximation for each value of x 1 . The likelihood term has a mode at (∞, ∞), hence the posterior is a compromise between this and the Gaussian prior centred at (0, 0). We first demonstrate that even if the Gaussian approximation matching the mode of π(x) is not so good, the Laplace approximation which uses a sequence of Gaussian approximations, can do much better. Let ρ = 1/2 and c = 10 (which is an extreme value). The resulting marginal for x 1 (solid), the Laplace approximation of it (dashed) and Gaussian approximation (dot-dashed), are shown in Figure 1. The Gaussian approximation fails both to locate the marginal correctly and, of course, it also fails to capture the skewness that is present. In spite of this, the sequence of Gaussian approximations used in the Laplace approximation performs much better and only seems to run into slight trouble where the curvature of the likelihood changes abruptly.
An important feature of (7) are its properties in the limiting cases ρ → 0 and ρ → 1. When ρ = 0, x 1 and x 2 become independent and π(x 2 |x 1 ) does not depend on x 1 . Hence, (7) is exact up to a numerical approximation of the normalising constant. In the other limiting case, ρ → 1, π(x 2 |x 1 ) is the point-mass at x 2 = x 1 , and (7) is again exact up numerical error. This illustrates the good property of (7), being exact in the two limiting cases of weak and strong dependence, respectively. This indicates that the approximation should not fail too badly for intermediate dependence.

Putting It All Together: INLA
With all the key components at hand, we now can put all these together to illustrate how they are combined to from INLA. The main aim of Bayesian inference is to approximate the posterior marginals π(θ j |y), j = 1, . . . , |θ|, π(x i |y), i = 1, . . . , n.
Our approach is tailored to the structure of LGMs, where |θ| is low-dimensional, x|θ is a GMRF and the likelihood is conditional independent in the sense that y i only depends on one x i and θ. From the discussion in Section 2.5, we know that we should aim to apply Laplace approximation only to near-Gaussian densities. For LGMs, it turns out that we can reformulate our problem as series of subproblems that allows us to use Laplace approximations on these. To illustrate the general principal, consider an artificial model is some well-behaved monotone function, and u ∼ N (0, Q −1 ). The index mapping j(i) is made such that the dimension of u is fixed and does not depend on n, and all u j s are observed roughly the same number of times. Computation of the posterior marginals for β and all u j is problematic, since we have a product of a Gaussian and a non-Gaussian (which is rather far from a Gaussian). Our strategy is to break down the approximation into smaller subproblems and only apply the Laplace approximation where the densities are almost Gaussian. They key idea is to use conditioning, here on β. Then The integral we need to approximate should be close to Gaussian, since the integrand is a Poissoncount correction of a Gaussian prior. The marginals for each u j , can be expressed as π(u j |y) = π(u j |β, y) × π(β|y) dβ. Note that we can compute the integral directly, since β is one-dimensional. Similar to (10), we have that which should be close to a Gaussian. Approximating π(u j |β, y) involves approximation of the integral of this density in one dimension less, since u j is fixed. Again, this is close to Gaussian.
The key lesson learnt, is that we can break down the problem into three sub-problems.
3. Compute π(u j |y) for all j using the results from the two first steps, combined with numerical integration (11).
The price we have to pay for taking this approach is increased complexity; for example step 2 needs to be computed for all values of β's that are required. We also need to integrate out the β's in (11), numerically. If we remain undeterred by the increased complexity, the benefit of this procedure is clear; we only apply Laplace approximations to densities that are near-Gaussians, replacing complex dependencies with conditioning and numerical integration. The big question is whether we can pursue the same principle for LGMs, and whether we can make it computationally efficient by accepting appropriate trade-offs that allow us to still be sufficiently exact. The answer is Yes in both cases. The strategy outlined above can be applied to LGMs by replacing β with θ, and u with x, and then deriving approximations to the Laplace approximations and the numerical integration. The resulting approximation is fast to compute, with little loss of accuracy. We will now discuss the main ideas for each step -skipping some practical and computational details that are somewhat involved but still relatively straight forward using "every trick in the book" for GMRFs.

Approximating the Posterior Marginals for the Hyperparameters
Since the aim is to compute a posterior for each θ j , it is tempting to use the Laplace approximation directly, which involves approximating the distribution of (θ −j , x)|(y, θ j ) with a Gaussian. Such an approach will not be very successful, since the target is and will not be very close to Gaussian; it will typically involve triplets like τ x i x j . Instead we can construct an approximation to π(θ|y) ∝ π(θ)π(x|θ)π(y|x, θ) π(x|θ, y) , in which the Laplace approximation requires a Gaussian approximation of the denominator Here, P (θ) = Q(θ) + diag(c(θ)), while µ(θ) is the location of the mode. The vector c(θ) contains the negative second derivatives of the log-likelihood at the mode, with respect to x i . There are two important aspects of (15).
1. It is a GMRF with respect to the same graph as from a model without observations y, so computationally it does not cost anything to account for the observations since their impact is a shift in the mean and the diagonal of the precision matrix.
2. The approximation is likely to be quite accurate since the impact of conditioning on the observations, is only on the "diagonal"; it shifts the mean, reduces the variance and might introduce some skewness into the marginals etc. Importantly, the observations do not change the Gaussian dependency structure through the terms x i x j Q ij (θ), as these are untouched.
Since |θ| is of low dimension, we can derive marginals for θ j |y directly from the approximation to θ|y. Thinking traditionally, this might be costly since every new θ would require an evaluation of (15) and the cost of numerical integration would still be exponential in the dimension. Luckily, the problem is somewhat more well-behaved, since the latent field x introduces quite some uncertainty and more "smooth" behaviour on the θ marginals. In situations where the central limit theorem starts to kick in, π(θ|y) will be close to a Gaussian. We can improve this approximation using variance-stabilising transformations of θ, like using log(precisions) instead of precisions, the Fisher transform of correlations etc. Additionally, we can use the Hessian at the mode to construct almost independent linear combinations (or transformations) of θ. These transformations really simplify the problem, as they tend to diminish long tails and reduce skewness, which gives much simpler and better-behaved posterior densities.
The task of finding a quick and reliable approach to deriving all the marginal distributions from an approximation to the posterior density (13), while keeping the number of evaluation points low, was a serious challenge. We did not succeed on this until several years after Rue et al. (2009), and after several failed attempts. It was hard to beat the simplicity and stability of using the (Gaussian) marginals derived from a Gaussian approximation at the mode. However, we needed to do better as these Gaussian marginals were not sufficiently accurate. The default approach used now is outlined in Martins et al. (2013, Sec. 3.2), and involves correction of local skewness (in terms of difference in scale) and an integration-free method to approximate marginals from a skewness-corrected Gaussian. How this is technically achieved is somewhat involved and we refer to Martins et al. (2013) for details. In our experience we now balance accuracy and computational speed well, with an improvement over Gaussian marginals while still being exact in the Gaussian limit.
In some situations, our approximation to (13) can be a bit off. This typically happens in cases with little smoothing and/or no replications, for example when η i = µ + β z z i + u i , for a random-effect u, and a binary likelihood . With vague priors model like this verge on being improper. Ferkingstad and Rue (2015) discuss these cases and derive a correction term which clearly improves the approximation to π(θ|y).

Approximating the Posterior Marginals for the Latent Field
We will now discuss how to approximate the posterior marginals for the latent field. For linear predictors with no attached observations, the posterior marginals are also the basis to derive the predictive densities, as the linear predictor itself is a component of the latent field. Similar to (11), we can express the posterior marginals as π(x i |y) = π(x i |θ, y) π(θ|y) dθ, hence we are faced with two more challenges.
1. We need to integrate over π(θ|y), but the computational cost of standard numerical integration is exponential in the dimension of θ. We have already ruled out such an approach in Section 3.1, since it was too costly computationally, except when the dimension is low. The key to success is to come up with efficient approximate solutions for each of these problems. Classical numerical integration is only feasible in lower dimensions. If we want to use 5 integration points in each dimension, the cost would be 5 k to cover all combinations in k dimensions, which is 125 (k = 3) and 625 (k = 4). Using only 3 integration points in each dimension, we get 81 (k = 4) and 729 (k = 6). This is close to the practical limits. Beyond these limits we cannot aim to do accurate integration, but should rather aim for something that is better than avoiding the integration step, like an empirical Bayes approach which just uses the mode. In dimensions > 2, we borrow ideas from central composite design (Box and Wilson, 1951) and use integration points on a sphere around the centre; see Figure 3 which illustrates the procedure in dimension 2 (even though we do not suggest using this approach in dimension 1 and 2). The integrand is approximately spherical (after rotation and scaling), and the integration points will approximately be located on an appropriate level set for the joint posterior of θ. We can weight the spherical integration points equally, and determine the relative weight with the central point requiring the correct expectation of θ T θ, if the posterior is standard Gaussian (Rue et al., 2009, Sec. 6.5). It is our experience that this approach balances computational costs and accuracy well, and it is applied as the default integration scheme. More complex integration schemes could be used with increased computational costs.
For the second challenge, we need to balance the need for improved approximations beyond the Gaussian for π(x i |θ, y), with the fact that we (potentially) need to do this n times. Since n can be large, we cannot afford doing too heavy computations for each i to improve on the Gaussian approximations. The default approach is to compute a Taylor expansion around the mode of the Laplace approximation, which provides a linear and a cubic correction term to the (standarized) Gaussian approximation, log π(x i |θ, y) ≈ − 1 2 We match a skew-Normal distribution (Azzalini and Capitanio, 1999) to (17), such that the linear term provides a correction term for the mean, while the cubic term provides a correction for skewness. This means that we approximate (16) with a mixture of skew-Normal distributions. This approach, termed simplified Laplace approximation, gives a very good trade-off between accuracy and computational speed. Additional to posterior marginals, we can also provide estimates of the deviance information criterion (DIC) (Spiegelhalter et al., 2002), Watanabe-Akaike information criterion (WAIC) (Wantanabe, 2010;Gelman et al., 2014), marginal likelihood and conditional predictive ordinates (CPO) . Other predictive criteria such as the ranked probability score (RPS) or the Dawid-Sebastiani-Score (DSS) (Gneiting and Raftery, 2007) can also be derived in certain settings (Riebler et al., 2012;Schrödle et al., 2012). Martins and Rue (2014) discuss how the INLA-framework can be extended to a class of near-Gaussian latent models.

THE R-INLA PACKAGE: EXAMPLES
The R-INLA package (see www.r-inla.org) provides an implementation of the INLA-approach, including standard and non-standard tools to define models based on the formula concept in R. In this section, we present some examples of basic usage and some special features of R-INLA.

A Simple Example
We first show the usage of the package through a simple simulated example, y|η ∼ Poisson(exp(η)) where η i = µ + βw i + u j(i) , i = 1, . . . , n, w are covariates, u ∼ N m (0, τ −1 I), and j(i) is a known mapping from 1 : n to 1 : m. We generate data as follows set.seed(123456L) n = 50; m = 10 w = rnorm(n, sd = 1/3) u = rnorm(m, sd = 1/4) intercept = 0; beta = 1 idx = sample(1:m, n, replace = TRUE) y = rpois(n, lambda = exp(intercept + beta * w + u[idx])) giving > table(y, dnn=NULL) 0 1 2 3 5 17 18 9 5 1 We use R-INLA to do the inference for this model, by library(INLA) my.data = data.frame(y, w, idx) formula = y~1 + w + f(idx, model="iid"), r = inla(formula, data = my.data, family = "poisson") The formula defines how the response depends on covariates, as usual, but the term f(idx, model="iid") is new. It corresponds to the function f that we have met above in (2), one of many implemented GMRF model components. The iid term refers to the N (0, τ −1 I) model, and idx is an index that specifies which elements of the model component go into the linear predictor. Figure 4a shows three estimates of the posterior marginal of u 1 . The solid line is the default estimate, the simplified Laplace approximation, as outlined in Section 3 (and with the R-commands given above). The dashed line is the simpler Gaussian approximation which avoids integration over θ, r.ga = inla(formula, data = my.data, family = "poisson", control.inla = list(strategy = "gaussian", int.strategy = "eb")) The dotted line represents the (almost) true Laplace approximations and accurate integration over θ, and is the best approximation we can provide with the current software, r.la = inla(formula, data = my.data, family = "poisson", control.inla = list(strategy = "laplace", int.strategy = "grid", dz=0.1, diff.logdens=20)) It is hard to see as it almost entirely covered by the solid line, meaning that our mixture of skew-Normals is very close to being exact in this example. We also note that by integrating out θ, the uncertainty increases, as it should. To compare the approximations with a simulation based approach, Figure 4b shows the corresponding histogram for 10 5 samples using JAGS, together with the default estimate from Figure 4a. The fit is quite accurate. The CPU time used by R-INLA with default options, was about 0.16 seconds on a standard laptop where 2/3 of this time was used for administration.

A Less Simple Example Including Measurement Error
We continue with a measurement error extension of the previous example, assuming that the covariate w is only observed indirectly through z, where , i = 1, . . . , n, with intercept γ. In this case, the model needs to be specified using two likelihoods and also a special feature called copy. Each observation can have its own type of likelihood (i.e. family), which is coded using a matrix (or list) of observations, where each "column" represents one family. A linear predictor can only be associated with one observation. The copy feature allows us to have additional identical copies of the same model component in the formula, and we have the option to scale it as well. An index NA is used to indicate if there is no contribution to the linear predictor and this is used to zero-out contributions from model components. This is done in the code below: ## generate observations that we observe for 'w' m = 2 z = rbinom(n, size = m, prob = 1/(1+exp(-(0 + w)))) ## create the response. since we have two families, poisson and ## binomial, we use a matrix, one column for each family Y = matrix(NA, 2*n, 2) Y[1:n , 1] = y Y[n + 1:n, 2] = z ## we need one intercept for each family. this is an easy way to achive that Intercept = as.factor(rep(1:2, each=n)) ## say that we have 'beta*w' only for 'y' and 'w' only for 'z'. the formula ## defines the joint model for both the observations, 'y' and 'z' NAs = rep(NA, n) idx = c(NAs, 1:n) idxx = c(1:n, NAs) formula2 = Y~-1 + Intercept + f(idx, model="iid") + f(idxx, copy="idx", hyper = list(beta = list(fixed = FALSE))) ## need to use a 'list' since 'Y' is a matrix my.data2 = list(Y=Y, Intercept = Intercept, idx = idx, idxx = idxx) ## we need to define two families and give the 'size' for the binomial r2 = inla(formula2, data = my.data2, family = c("poisson", "binomial"), Ntrials = c(NAs, rep(m, n))) We refer to  for more details on measurement error models using INLA, and to the specific latent Gaussian models termed mec and meb that are available in R-INLA to facilitate the implementation of classical error models and Berkson error models, respectively.

A Spatial Example
The R-INLA package has extensive support for spatial Gaussian models, including intrinsic GMRF models on regions (often called "CAR" models, (Hodges, 2013, Ch. 5.2)), and a subclass of continuously indexed Gaussian field models. Of particular interest are Gaussian fields derived from stochastic partial differential equations (SPDEs). The simplest cases are Matérn fields in dimension d, which can be described as the solution to where ∆ is the Laplacian, κ > 0 is the spatial scale parameter, α controls the smoothness, τ controls the variance, and W(s) is a Gaussian spatial white noise process. Whittle (1954Whittle ( , 1963 shows that its solution is a Gaussian field with a Matérn covariance function having smoothness ν = α − d/2. The smoothness is usually kept fixed based on prior knowledge of the underlying process. A formulation of Matérn fields as solutions to (18) might seem unnecessarily complicated, since we already know the solution. However, Lindgren et al. (2011) showed that by using a finite basis-function representation of the continuously indexed solution, one can derive (in analogy to the well known Finite Element Method ) a local representation with Markov properties. This means that the joint distribution for the weights in the basis-function expansion is a GMRF, and the distribution follows directly from the basis functions and the triangulation of space. The main implication of this result is that it allows us to continue to think about and interpret the model using marginal properties like covariances, but at the same time we can do fast computations since the Markov properties make the precision matrix very sparse. It also allows us to add this component in the R-INLA framework, like any other GMRF model-component. The dual interpretation of Matérn fields, both using covariances and also using its Markov properties, is very convenient both from a computational but also from a statistical modeling point of view (Simpson et al., 2011(Simpson et al., , 2012Lindgren and Rue, 2015). The same ideas also apply to non-stationary Gaussian fields using non-homogeneous versions of an appropriate SPDE Fuglstad et al., 2015a,b;Yue et al., 2014), Gaussian fields that treats land as a barrier to spatial correlation (Bakka et al., 2016), multivariate random fields (Hu and Steinsland, 2016), log-Gaussian Cox processes (Simpson et al., 2016a), and in the near future also to non-separable space-time models.
We end this section with a simple example of spatial survival analysis taken from Henderson et al. (2002), studying spatial variation in leukaemia survival data in north-west England in the period 1982-1998. The focus of the example is to see how and how easily, the spatial model integrates into the model definition (Martino et al., 2010). We therefore omit further details about the dataset and refer to the original article.
First, we need to load the data and create the mesh, i.e. a triangulation of the area of interest to represent the finite dimensional approximation to (18). library(INLA) data(Leuk) loc <-cbind(Leuk$xcoord, Leuk$ycoord) bnd1 <-inla.nonconvex.hull(loc, convex=0.05) bnd2 <-inla.nonconvex.hull(loc, convex=0.25) mesh <-inla.mesh.2d(loc, boundary=list(bnd1, bnd2), max.edge=c(0.05, 0.2), cutoff=0.005) Figure 5a displays the study area and the locations of the events, while Figure 5b shows the associated mesh with respect to which we define the SPDE model. We use an additional rougher mesh to reduce boundary effects. The next step is to create a mapping matrix from the mesh onto the locations where the data are observed. Then we define the SPDE model, to define the statistical model including covariates like sex, age, white blood-cell counts (wbc) and the Townsend deprivation index (tpi), and to call a book-keeping function which keeps the indices in correct order. Finally, we call inla() to do the analysis, assuming a Weibull likelihood. Note that application of a Cox proportional hazard model will give similar results.

Special Features
In addition to standard analyses, the R-INLA package also contains non-standard features that really boost the complexity of models that can be specified and analysed. Here, we give a short summary of these, for more details see Martins et al. (2013).
replicate Each model component given as a f()-term can be replicated, creating nrep iid replications with shared hyperparameters. For example, f(time, model="ar1", replicate=person) defines one AR(1) model for each person sharing the same hyperparameters.
group Each model component given as a f()-term, can be grouped, creating ngroup dependent replications with a separable correlation structure. To create a separable space-time model, with an AR(1) dependency in time, we can specify f(space, model=spde, group=time, control.group = list(model = "ar1")) Riebler et al. (2012) used grouped smoothing priors in R-INLA to impute missing mortality rates for a specific country by taking advantage from similar countries where these data are available. The authors provide the corresponding R-code in the supplementary material. We can both group and replicate model components. A-matrix We can create a second layer of linear predictors where η is defined by the formula, but where η * = Aη is connected to the observations. Here, A is a constant (sparse) matrix; see the above spatial example.
Linear combinations We can also compute posterior marginals of v = Bx where x is the latent field and B is a fixed matrix. This could for example be β 1 − β 2 for two fixed effects, or any other linear combinations. Here is an example computing the posterior for the difference between two linear effects, β u − β v lc = inla.make.lincomb(u=1, v=-1) r = inla(y~u + v, data = d, lincomb = lc) Remote server It is easy to set up a remote MacOSX/Linux server to host the computations while doing the R-work at your local laptop. The job can be submitted and the results can be retrieved later, or we can use it interactively. This is a very useful feature for larger models. It also ensures that computational servers will in fact be used, since we can work in a local R-session but use a remote server for the computations. Here is an example running the computations on a remote server r = inla(formula, family, data = data, inla.call = "remote") To submit a job we specify r = inla(formula, family, data = data, inla.call = "submit") and we can check the status and retrieve the results when the computations are done, by inla.qstat(r) r = inla.qget(r) R-support Although the core inla-program is written in C, it is possible to pass a user-defined latent model component written in R, and use that as any other latent model component.
The R-code will be evaluated within the C-program. This is very useful for more specialised model components or re-parameterisations of existing ones, even though it will run slower than a proper implementation in C. As a simple example, the code below implements the model component iid, which is just independent Gaussian random effects N n (0, (τ I) −1 ). The skeleton of the function is predefined, and must return the graph, the Q-matrix, initial values, the mean, the log normalising constant and the log prior for the hyperparameter.

A CHALLENGE FOR THE FUTURE: PRIORS
Although the R-INLA project has been highly successful, it has also revealed some "weak points" in general Bayesian methodology from a practical point of view. In particular, our main concern is how we think about and specify priors in LGMs. We will now discuss this issue and our current plan to provide good sensible "default" priors.
Bayesian statistical models require prior distributions for all the random elements of the model. Working within the class of LGMs, this involves choosing priors for all the hyperparameters θ in the model, since the latent field is by definition Gaussian. We deliberately wrote priors since it is common practice to define independent priors for each θ j , while what we really should aim for is a joint prior for all θ, when appropriate.
The ability to incorporate prior knowledge in Bayesian statistics is a great tool and potentially very useful. However, except for cases where we do have "real/experimental" prior knowledge, for example through results from previous experiments, it is often conceptually difficult to encode prior knowledge through probability distributions for all model parameters. Examples include priors for precision and overdispersion parameters, or the amount of t-ness in the Student-t distribution. Simpson et al. (2016b) discuss these aspects in great detail.
In R-INLA we have chosen to provide default prior distributions for all parameters. We admit that currently these have been chosen partly based on the priors that are commonly used in the literature and partly out of the blue. It might be argued that this is not a good strategy, and that we should force the user to provide the complete model including the joint prior. This is a valid point, but all priors in R-INLA can easily be changed, allowing the user to define any arbitrary prior distribution. So the whole argument boils down to a question of convenience.
Do we have a "Houston, we have a problem"-situation with priors? Looking at the current practice within the Bayesian society, we came to the conclusion; we do. We will argue for this through a simple example, showing what can go wrong, how we can think about the problem and how we can fix it. We only discuss proper priors.
Consider the problem of replacing a linear effect of the Townsend deprivation index tpi with a smooth effect of tpi in the Leukaemia example in Section 4.3. This is easily implemented by replacing tpi with f(tpi, model="rw2"). Here, rw2 is a stochastic spline, simply saying that the second derivative is independent Gaussian noise (Rue and Held, 2005;Lindgren and Rue, 2008). By default, we constrain the smooth effect to also sum to zero, so that these two model formulations are the same in the limit as the precision parameter τ tends to infinity, and a vague Gaussian prior is used for the linear effect. The question is which prior should be used for τ . An overwhelming majority of cases in the literature uses some kind of a Gamma(a, b) prior for τ , implying that π(τ ) ∝ τ a−1 exp(−bτ ), for some a, b > 0. This prior is flexible, conjugate with the Gaussian, and seems like a convenient choice. Since almost everyone else is using it, how wrong can it be?
If we rewind to the point where we replaced the linear effect with a smooth effect, we realise that we do this because we want a more flexible model than the linear effect, i.e. we also want to capture deviations from the linear effect. Implicitly, if there is a linear effect, we do want to retrieve that with enough data. Measuring the distance between the straight line and the stochastic spline using the Kullback-Leibler divergence, we find that KLD ∝ 1/τ meaning that the (unidirectional) distance is d ∝ 1/τ . For simplicity, choose a = b = 1 in the Gamma-prior, then the derived prior for the distance d is Figure 7a displays this prior on the distance scale, revealing two surprising features. First, the mode is around d ≈ 0.82, and second, the prior appears to be zero for a range of positive distances. The second feature is serious as it simply prevents the spline from getting too close to the linear effect. It is clear from (19) that the effect is severe, and in practice, π(d) ≈ 0 even for positive d. This is an example of what Simpson et al. (2016b) call prior overfitting; the prior prevents the simpler model to be located, even when it is the true model. Choosing different parameters in the Gamma-prior does not change the overfitting issue. For all a, b > 0, the corresponding prior for the distance tends to 0 as d → 0. For a (well-behaved) prior to have π(d = 0) > 0, we need E(τ ) = ∞.
If we are concerned about the behaviour of the distance between the more flexible and the simpler model component, we should define the prior directly on the distance, as proposed in Simpson et al. (2016b). A prior for the distance should be decaying with the mode at distance zero. This makes the simpler model central and the point of attraction. The exponential prior is recommended as a generic choice since it has a constant rate penalisation, π(d) = λ exp(−λd). The value of λ could be chosen by calibrating some property of the model component under consideration. Note that this way of defining the prior is invariant to reparameterisations, as it is defined on the distance and not for a particular parametersation. Let us return to the stochastic spline example, assigning the exponential prior to the distance. The parameter λ can be calibrated by imposing the knowledge that the effect of tpi is not likely to be above 1 on the linear predictor scale, ..+ f(tpi, model="rw2", scale.model = TRUE, hyper = list(prec = list(prior="pc.prec", param=c(1, 0.01)))) Here, scale.model is required to ensure that the parameter τ represents the precision, not just a precision parameter (Sørbye and Rue, 2014). The estimated results are given in Figure 7b, illustrating the point-wise posterior mean, median and the 2.5% and 97.5% credibility intervals, for the effect of tpi on the mean survival time.
Here, we have only briefly addressed the important topic of constructing well-working priors, and currently we are focusing a lot of activity on this issue to take the development further. Besides others we plan to integrate automatic tests for prior sensitivity, following the work of Roos and Held (2011);Roos et al. (2015a). The final goal is to use the above ideas to construct a joint default prior for LGMs, which can be easily understood and interpreted. A main issue is how to decompose and control the variance of the linear predictor, an issue we have not discussed here. For further information about this issue, please see Simpson et al. (2016b) for the original report which introduces the class of penalised complexity (PC) priors. Some examples on application of these priors include disease mapping , bivariate meta-analysis , age-period-cohort models (Riebler and Held, 2016), Bayesian P-splines (Ventrucci and Rue, 2016), structured additive distributional regression (Klein and Kneib, 2016), Gaussian fields in spatial statistics (Fuglstad et al., 2016), modeling monthly maxima of instantaneous flow (Ferkingstad et al., 2016) and autoregressive processes .
Interestingly, the framework and ideas of PC priors, are also useful for sensitivity analysis of model assumptions and developing robust models, but it is too early to report this here. Stay tuned!

DISCUSSION
We hope we have convinced the reader that the INLA approach to approximate Bayesian inference for LGMs is a useful addition to the applied statistician's toolbox; the key components just play so nicely together, providing a very exact approximation while reducing computation costs substantially.
The key benefit of the INLA approach is that it is central to our long-term goal of making LGMs a class of models that we (as a community) can use and understand.
Developing, writing and maintaining the code-base for a such large open-source project, is a huge job. Nearly all the R/C/C++ code is written and maintained by F. Lindgren (20%) and H. Rue (80%), and is a result of a substantial amount of work over many years. Many more have contributed indirectly by challenging the current practice and implementation. The current version of this project is a result of the cumulative effort of the many users, and their willingness to share, challenge and question essentially everything. Documentation is something we could and should improve upon, but the recent book by Blangiardo and Cameletti (2015) does a really good job.
The current status of the package is good, but we have to account for the fact that the software has been developed over many years, and is basically the version we used while developing the methods. Hence, while the software works well it less streamlined and less easy to maintain than it ought to be. We are now at a stage where we know what we want the package to do and software to be, hence a proper rewrite by skilled people would really be a useful project for the society. If this would happen, we would be more than happy to share all our knowledge into a such "version 2.0" project! Another use of R-INLA is to use it purely as computational back-end. The generality of R-INLA comes with a prize of complexity for the user, hence a simplified interface for a restricted set of models can be useful to improve accessibility for a specific target audience or provide additional tools that are mainly relevant for these models. Examples of such projects, are AnimalINLA (Holand et al., 2013), ShrinkBayes (Van De Wiel et al., 2013a,b, 2014Riebler et al., 2014), meta4diag , BAPC (Riebler and Held, 2016), diseasemapping and geostatp (Brown, 2015), and Bivand et al. (2015). Similarly, the excursions package for calculating joint exceedance probabilities in GMRFs Lindgren, 2015, 2016) includes an interface to analyse LGMs estimated by R-INLA. Recent work on methodology for filtered spatial point patterns in the context of distance sampling (Yuan et al., 2016) has initiated the construction of wrapper software for fitting other complex spatial models such as those resulting from plot sampling data or for point process models within R-INLA. There is also an interesting line of research using R-INLA to do approximate inference on a sub-model within a larger model, see Guihenneuc-Jouyaux and Rousseau (2005) for a theoretical justification and Li et al. (2012) for an early application of this idea. One particular application here, is how to handle missing data in cases where the joint model is not an LGM.
Please visit us at www.r-inla.org!