alexa Bayesian Nonparametric Regression and Density Estimation Using Integrated Nested Laplace Approximations | Open Access Journals
ISSN: 2155-6180
Journal of Biometrics & Biostatistics
Like us on:
Make the best use of Scientific Research and information from our 700+ peer reviewed, Open Access Journals that operates with the help of 50,000+ Editorial Board Members and esteemed reviewers and 1000+ Scientific associations in Medical, Clinical, Pharmaceutical, Engineering, Technology and Management Fields.
Meet Inspiring Speakers and Experts at our 3000+ Global Conferenceseries Events with over 600+ Conferences, 1200+ Symposiums and 1200+ Workshops on
Medical, Pharma, Engineering, Science, Technology and Business

Bayesian Nonparametric Regression and Density Estimation Using Integrated Nested Laplace Approximations

Xiao-Feng Wang*

Department of Quantitative Health Sciences/Biostatistics Section, Cleveland Clinic Lerner Research Institute, Cleveland, OH 44195, USA

*Corresponding Author:
Xiao-Feng Wang
Department of Quantitative Health Sciences
Cleveland Clinic Lerner Research Institute
9500 Euclid Avenue/JJN3, Cleveland, OH 44195, USA
E-mail: [email protected]

Received Date: June 18, 2013; Accepted Date: June 20, 2013; Published Date: June 25, 2013

Citation: Wang XF (2013) Bayesian Nonparametric Regression and Density Estimation Using Integrated Nested Laplace Approximations. J Biomet Biostat 4:e125. doi: 10.4172/2155-6180.1000e125

Copyright: © 2013 Wang XF. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Visit for more related articles at Journal of Biometrics & Biostatistics


Nonparametric regression; Density estimation; Approximate Bayesian inference; Integrated nested Laplace approximations; Markov chain Monte Carlo

INLA and Approximate Bayesian Inference

Structured additive regression models have been extensively used in many fields, such as medicine, public health, and economics. In these models, the response variable yi is assumed to belong to an exponential family and its mean μi is linked to a structured additive predictor ηi through a link function g(.), where

Equation (1)

Here α0 is the intercept, αk’s are the linear effects of covariates xk’s, fl(.)’s are unknown functions of the covariates μl’s and εi’s is unstructured error term. In Bayesian statistics, the common solution to obtain estimates for the models is the use of Markov chain Monte Carlo (MCMC) techniques. Although it is always possible to implement them in theory, MCMC methods could come with a wide range of issues in practice: parameter samples could be highly correlated, computational time is very long, and estimates may content large Monte Carlo errors, etc.

Approximate Bayesian inference using integrated nested Laplace approximations (INLA) is a recently proposed method for solving the structured additive regression models where the latent field of the models is Gaussian [1]. The methodology is particularly useful, when the latent Gaussian model is a Gaussian Markov random field (GMRF) with precision matrix G controlled by a few hyperparameters. Denote that β is the vector of all unknown regression parameters in the structured additive predictor (1), assigned Gaussian priors; ξ is the vector of hyperparameters for which non-Gaussian priors are assigned in the model. In a Bayesian framework, one is to estimate the marginal posterior density

Equation (2)

given data y for each component βj of the Gaussian vector β. The key of INLA is to construct a nested approximation of (2). First, the marginal posterior Equation is approximated using


where Equation is the Gaussian approximation to the full conditional distribution of β evaluated in the mode β*(ξ) for a given ξ. Then, one computes Equation an approximation of the posterior Equation Rue et al. [1] suggested three approximation approaches: a Gaussian approximation, a full Laplace approximation, and a simplified Laplace approximation. Lastly, INLA approach combines the previous two steps with a numerical integration to reach the goal. The approximate posterior of Equation is obtained by


where the sum is over values of ξ with area weights Δq .

INLA method provides accurate approximations to the posterior marginals of the latent variables which are extremely fast to compute. No samples of the posterior marginal distributions need to be drawn using INLA, so it is a computationally convenient alternative to MCMC. All computations required by the INLA methodology have been efficiently implemented by a R package INLA (available on It integrates GMRFLib, a C-based library for fast and exact simulation of GMRF.

Rue et al. [1] demonstrated that INLA can be applied to solve a variety of popular statistical models, which includes generalized linear mixed models, stochastic volatility models, spatial disease mapping, and Log-Gaussian Cox processes. In this paper, we show that two classical nonparametric smoothing problems, nonparametric regression and probability density estimation, can be achieved by using INLA. Simulated examples and R functions are demonstrated to illustrate the use of the methods. The future research on the applications of INLA to other complex statistical models is discussed in the paper.

Nonparametric Regression Using INLA

Nonparametric regression is a classical problem in statistics, which does not require to specify a parametric functional form for the relationship between the response y and the predictor x. We assume a model of the form


where m(.) is an unknown smooth function and the errors εi’s are assumed to be independent and identically distributed as N(0, s2).

There are several popular techniques to estimate nonparametric regression models, such as, local polynomial regression methods [2], and smoothing spline approaches [3]. Here we examine the use of spline techniques. Smoothing spline is motivated by considering the penalized residual sum of squares as a criterion for the goodness of fit of m(x) to the data. It minimizes


where λ is a smoothing parameter to control the trade-off between fit and the penalty Equation. The solution to this minimization problem is a natural cubic spline with knots at the distinct observed values of Xi.

Wahba [4] showed that the smoothing spline is equivalent to Bayesian estimation with a partially improper prior. m(x) has the prior distribution which is the same as the distribution of the stochastic process


where Equation is fixed, and V(x) is the one-fold integrated Wiener process,


Thus, estimating m(x) becomes to seek the solution to the stochastic differential equation

Equation (3)

as a prior over m. Note that such a differential equation of order two is the continuous-time version of a second-order random walk. However, the solution of (3) does not have any Markov properties. The precision matrix is dense, hence it is computationally intensive. Lindgren and Rue [5] suggested a Galerkin approximation to m(x), as the solution of (3). To be specific, let x1 < x2 < ... <xn be the set of (unequal-spaced) fixed points, a finite element representation of m(x) is constructed as


for the piecewise linear basis functions ψi’s and random weights wi’s.

In order to estimate the smooth function m(x), one needs to determine the joint distribution of the weights w = (w1,....,wn)T Using the Galerkin method, w is derived as a GMRF with mean zero and precision matrix G. Let di =xi+1-xi for i =1,..., n -1 and d-1= d0 = dn = dn+1 = ∞.

The n× n symmetric matrix G is defined as


where the non-zero elements of row i are given by


with gi,i+1 Ξ gi+1,i and gi,i+2 Ξ gi+2,i due to symmetry. G is a sparse matrix with rank n - 2, making the model computationally effective.

If we assign Equationas a smoothness prior over m, the cubic smoothing spline Equation at x coincides with the posterior expectation of m(x) give the data, i.e. Equation Therefore, the nonparametric regression problem becomes to fit a latent Gaussian model. It can be accomplished using INLA since w is a GMRF. The implementation of the method needs some extra programming in R to define a userspecified GMRF. A R function called npr.INLA is listed in Appendix. The following simple code is to call the fit using INLA, where “x” and “y” are inputs of numeric vector of predictors and responses, respectively.

R> library (INLA)

R> library (splines)

R> fit1<-npr.INLA (x,y)

The outputs of npr.INLA include a vector of sorted x values at which the estimate is computed, a vector of smoothed estimates for the regression at the corresponding x, and two vectors of the corresponding 2.5% and 97.5% credible bands for the regression.

We now show two simulated examples to compare the INLA approach and the conventional cubic smoothing spline regression. Generalized cross-validation method is applied to select the smoothing parameters for smoothing splines. The Bayesian method takes the advantage that the smoothing parameter is automatically determined by the model fitting without any user input. In the first example, m(x) = 3sin(2.5x) + 2exp(-5x2), x~U(0,1.5), and ε ~ N(0,0.52) with sample size n=50. In the second example, m(x) = 3x2/5 - cos(πx), x ~U(-3.0,0), and ε ~ N(0,0.82) with sample size n=100. Figure 1 shows the estimation results. The estimates using INLA are denoted by the solid lines, and the corresponding 95% credible bands are denoted by the dashed lines. The estimates using the cubic smoothing spline regression are denoted by the dash-dotted lines. The true functions are denoted by the dotted lines. Unsurprisingly, the estimates using INLA are almost identical to those using cubic smoothing spline. The 95% credible bands from the INLA fits completely cover the true functions in these two examples. Approximate Bayesian inference through INLA allows fast Bayesian computation and makes it possible to perform analysis in an automatic way.


Figure 1: Simulated examples for nonparametric regression: (a) m(x) = 3sin(2.5x) + 2exp(−5x ), n = 50,σ = 0.5; (b) m(x) = x2/5 − cos(π x), n =100,σ = 0.8. The estimates (solid lines) with credible bands (dashed lines) using INLA are compared to the estimates (dash-dotted lines) using cubic smoothing spline regression. The true functions are denoted by the dotted lines.

Density Estimation Using INLA

Nonparametric density estimation can be also implemented using INLA. Brown et al. [6] proposed a “root-unroot” density estimation procedure, which turned density estimation into a nonparametric regression problem. The regression problem was created by binning the original observations into suitable size of bins and applying a mean-matching variance stabilizing root transform to the binned data counts. Then, a wavelet block thresholding regression was used to obtain the density estimate. Here we adopt Brown et al. [6] rootunroot procedure but use a second-order random walk model with INLA for the regression step. The second-order random walk model is particularly suitable for an equi-spaced nonparametric time series regression problem [7]. In addition, there are two advantages to use the Bayesian nonparametric approach. First, we avoid the smoothing parameter selection, where the smoothness of curve is automatically determined by the Bayesian model fitting. Second, it is straightforward to construct credible bounds of a regression curve from INLA output. As a result, constructing credible bands for the probability density function becomes a natural by-product in the density estimation. Let Equation be a random sample from a distribution with the density function fx. The estimation algorithm is summarized as follows.

1. Poissonization. Divide Equation in T equal length intervals. Let Equation be the count of observations in each of the intervals.

2. Root Transformation. Apply the mean-matching variance stabilizing root transform, Equation

3. Bayesian Smoothing with INLA. Consider the time series Equation to be the sum Equation of a smooth trend function m and a noise component ε. Fit a second-order random walk model with INLA for the equi-spaced time series to obtain an posterior mean estimate Equation of m, and α / 2 and 1-α / 2 quantiles, Equation

4. Unroot Transformation and Normalization. The density function fx is estimated by


and the 100(1-α )% credible bands of f(x) is


where Equation is a normalization constant.

The R function, density.INLA in Appendix implements the above root-unroot algorithm. The following code is to call the density fit using INLA, where “x” is a numeric vector of data values, “m” is the number of equally spaced points at which the density is to be estimated, “from” and “to” are the left and right-most points of the grid at which the density is to be estimated. If “from” and “to” are missing, “from” equals to the minima of the data values minus “cut” times the range of the data values and “to” equals to the maxima of the data values plus “cut” times the range of the data values.

R> library (INLA)

R> fit1 <-density.INLA (x, m=101, from=min (x), to=max (x), cut=0.1)

Figure 2 shows the two simulated examples to compare INLA method and conventional kernel density estimation. In the first example data were generated from the standard normal distribution, X ~ N(0,1),with sample size n=500; In the second example data were generated from a normal mixture model, X ~ 0.5N(-1.5,1) + 0.5N(2.5,0.752 ), with sample size n=1000. In the figure, the estimates using INLA are denoted by the solid lines, and the 95% credible bands are denoted by the dashed lines. The kernel density estimates with Sheather and Jones [8] plug-in bandwidth are denoted by the dash-dotted lines. The true functions are denoted by the dotted lines. We note that the INLA estimates are very close to the kernel density estimates. The INLA approach allows us to compute the credible bands of the density function without additional computational efforts.


Figure 2: Simulated examples for density estimation: (a) X ~ N(0,1), n = 500; (b) X ~ 0.5N(−1.5,1) + 0.5N(2.5,0.752 ), n =1000 . The estimates (solid lines) using INLA are compared to the estimates (dash-dotted lines) using kernel density estimation with Sheather and Jones [8]’s plug-in bandwidth. The true functions are denoted by the dotted lines.


MCMC techniques were commonly choices to fit structured additive regression models, however they may often suffer the issues of convergence and computational time. INLA approach provides a novel approximate Bayesian inference method to fit a large class of structured additive models with latent Gaussian field. It, as an alternative to MCMC, provides accurate approximations to estimate posterior marginals and avoid time-consuming sampling. INLA method has been implemented in C and a R-interfaced package is available under Linux, Windows and Macintosh. The package provides a user-friendly interface and makes latent Gaussian models applicable in a general way.

We have shown that two classical nonparametric models can be fit through the approximate Bayesian inferential procedure. Nonparametric regression is treated within a general second-order random walk model by assigning appropriate GMRF priors. Density estimation is transformed to an equally-spaced nonparametric Bayesian time series regression problem by adapting Brown et al. [6] root-unroot algorithm. The approximate Bayesian inference of the nonparametric models enjoys the advantage of fast computation, automatic selection of smoothing parameter and construction of credible bands.

Recently, Lindgren et al. [9] have addressed the explicit link between Gaussian fields and Gaussian Markov random fields through the stochastic partial differential equation approach. INLA approach has been applied to a variety of complex models, such as spatialtemporal disease mapping [10], additive mixed quantile regression models [11], and beta semiparametric mixed models [12]. INLA, as new Bayesian computation tools, have a great potential to be used in many scientific fields. There are many open problems of applying INLA to advanced statistical models. For instance, it would be of interest to use INLA for joint modeling of longitudinal and time-to-event data, functional data analysis, measurement error models, and sparse ultrahigh- dimensional modeling. Such problems need further investigation and should be detailed in the future.

Appendix: R Functions

We present here the key compute code to call INLA for nonparametric regression and density estimation. Simulation examples and code files can be obtained from the journal’s webpage or

## Nonparametric regression using INLA

npr.INLA <- function(x, y, diagonal = 1e-03, constr = T, ...){

if (any( stop(‘x’ contains missing values!”)

if (any( stop(‘y’ contains missing values!”)

if (length(x) != length(y)) stop(‘x’ and ‘y’ have different lengths!”)

y <- y[order(x)]

x <- sort(x)

B <- bs(x, degree = 1, intercept = TRUE)

idx <- 1:length(x)

Q <- Gmatrix(x) <- inla(y ~ B.1 + B.2 + f(idx, model = ‘generic0”, Cmatrix = Q,

diagonal = diagonal, constr = constr, ...) - 1,

data = = y, idx = idx, B = B)),

control.predictor = list(compute = T))

return(structure(list(x = x, y =$summary.linear. predictor[,1],

y.lower =$summary.linear.predictor[,3],

y.upper =$summary.linear.predictor[,5])))


## Density estimation using INLA

density.INLA <- function(x, m = 101, from, to, cut = 0.1,

diagonal = 1e-03, constr = T, ...){

if (any( stop(‘x’ contains missing values!”)

if (missing(from)) from <- min(x) - cut * diff(range(x))

if (missing(to)) to <- max(x) + cut * diff(range(x))

bins=seq(from, to, length.out = m)

x.bins <- hist(x, breaks = bins, plot=FALSE)

x.bins.root <- sqrt(x.bins$counts+1/4)

idx <- 1:length(x.bins.root)

Q <- Gmatrix(idx) <- inla(x.bins.root ~ f(idx, model = “generic0”, Cmatrix = Q,

diagonal = diagonal, constr = constr),

data = = x.bins.root, idx = idx)),

control.predictor = list(compute = T))

inla.est <-$summary.linear.predictor[,1]

inla.lower <-$summary.linear.predictor[,3]

inla.upper <-$summary.linear.predictor[,5]

SimpsonInt <- function (x, f, subdivisions = 256){

ap <- approx(x, f, n = 2 * subdivisions + 1)

integral <- diff(ap$x)[1] * (ap$y[2 * (1:subdivisions) - 1]

+ 4 * ap$y[2 * (1:subdivisions)] + ap$y[2 * (1:subdivisions) + 1])/3



normalized <- SimpsonInt(x.bins$mids, inla.estˆ2)

f <- inla.estˆ2/normalized

f.lower <- inla.lowerˆ2/normalized

f.upper <- inla.upperˆ2/normalized

return(structure(list(x = x.bins$mids, y = f, y.lower= f.lower, y.upper=f.upper)))


## Function to compute G matrix

Gmatrix <- function(x, sparse = TRUE){

if (any( stop(“’x’ contains missing values!”)

x <- sort(x)

n <- length(x)

d <- diff(x)

d <- c(Inf, Inf, d, Inf, Inf)

k <- 3:(n + 2)

g <- (2/((d[k - 1]ˆ2) * (d[k - 2] + d[k - 1]))

+ 2/(d[k - 1]*d[k]) * (1/d[k - 1] + 1/d[k])

+ 2/((d[k]ˆ2) * (d[k] + d[k + 1])))

k <- 4:(n + 2)

g1 <- -2/(d[k - 1]ˆ2) * (1/d[k - 2] + 1/d[k])

k <- 5:(n + 2)

g2 <- 2/(d[k - 2] * d[k - 1] * (d[k - 2]+d[k - 1]))

G <- diag(g)

G[row(G) == col(G) + 1] <- g1

G[col(G) == row(G) + 1] <- g1

G[row(G) == col(G) + 2] <- g2

G[col(G) == row(G) + 2] <- g2

if (sparse == TRUE) G <- as(G, “dgTMatrix”)




Select your language of interest to view the total content in your interested language
Post your comment

Share This Article

Relevant Topics

Recommended Conferences

Article Usage

  • Total views: 11688
  • [From(publication date):
    October-2013 - Oct 20, 2017]
  • Breakdown by view type
  • HTML page views : 7902
  • PDF downloads :3786

Post your comment

captcha   Reload  Can't read the image? click here to refresh

Peer Reviewed Journals
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals
International Conferences 2017-18
Meet Inspiring Speakers and Experts at our 3000+ Global Annual Meetings

Contact Us

Agri, Food, Aqua and Veterinary Science Journals

Dr. Krish

[email protected]

1-702-714-7001 Extn: 9040

Clinical and Biochemistry Journals

Datta A

[email protected]

1-702-714-7001Extn: 9037

Business & Management Journals


[email protected]

1-702-714-7001Extn: 9042

Chemical Engineering and Chemistry Journals

Gabriel Shaw

[email protected]

1-702-714-7001 Extn: 9040

Earth & Environmental Sciences

Katie Wilson

[email protected]

1-702-714-7001Extn: 9042

Engineering Journals

James Franklin

[email protected]

1-702-714-7001Extn: 9042

General Science and Health care Journals

Andrea Jason

[email protected]

1-702-714-7001Extn: 9043

Genetics and Molecular Biology Journals

Anna Melissa

[email protected]

1-702-714-7001 Extn: 9006

Immunology & Microbiology Journals

David Gorantl

[email protected]

1-702-714-7001Extn: 9014

Informatics Journals

Stephanie Skinner

[email protected]

1-702-714-7001Extn: 9039

Material Sciences Journals

Rachle Green

[email protected]

1-702-714-7001Extn: 9039

Mathematics and Physics Journals

Jim Willison

[email protected]

1-702-714-7001 Extn: 9042

Medical Journals

Nimmi Anna

[email protected]

1-702-714-7001 Extn: 9038

Neuroscience & Psychology Journals

Nathan T

[email protected]

1-702-714-7001Extn: 9041

Pharmaceutical Sciences Journals

John Behannon

[email protected]

1-702-714-7001Extn: 9007

Social & Political Science Journals

Steve Harry

[email protected]

1-702-714-7001 Extn: 9042

© 2008-2017 OMICS International - Open Access Publisher. Best viewed in Mozilla Firefox | Google Chrome | Above IE 7.0 version