library(knitr) opts_chunk$set(cache=FALSE,comment=NA, fig.path="figs/", warning=FALSE, message=FALSE, optipng='-o7', pngquant='--speed=1 --quality=0-50') options(digits=5,show.signif.stars=FALSE,width=120) knit_hooks$set(optipng = hook_optipng) knit_hooks$set(pngquant = hook_pngquant) knitr::knit_hooks$set(setPch = function(before, options, envir) { if(before) par(pch = 20) }) opts_chunk$set(setPch = TRUE)
See the introduction for an overview. Load the libraries:
library(INLA) library(brinla) library(ggplot2)
We illustrate the method using the following data:
data(reeds, package="brinla") ggplot(reeds,aes(site,nitrogen))+geom_point()
We will build a model of the form:
nitrogen = intercept + site + error
where intercept
is a fixed effect modeled with a normal prior with a zero mean and large variance.
The site
is a random effect and error
is the error term. These have distributions which are
normal with mean zero and some variance. These variances are treated as hyperparameters with
priors we shall specify.
By default, INLA places a gamma prior on the precisions (inverse of the variance). We fit this model and examine the posteriors in terms of the SDs (which are easier to understand than precisions)
formula <- nitrogen ~ 1 + f(site, model="iid") imod <- inla(formula,family="gaussian", data = reeds) bri.hyperpar.summary(imod)
We can also plot the posteriors of these two SDs.
bri.hyperpar.plot(imod)
INLA has other choices for the priors. Here is an example using the penalized complexity prior. We use the SD of the response to help us set the scale of this prior (although it is better if you set this with knowledge of the problem behind the data)
sdres <- sd(reeds$nitrogen) pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdres,0.01))) formula <- nitrogen ~ f(site, model="iid", hyper = pcprior) pmod <- inla(formula, family="gaussian", data=reeds) bri.hyperpar.summary(pmod)
Plot the posterior densities.
bri.hyperpar.plot(pmod)
It is possible to set our own prior for the SD of the site
effect. The use of the
half(ie. positive part of) Cauchy is a commonly used choice which is not directly programmed
in INLA. The density for the SD ($\sigma$) with scaling $\lambda$ is:
$$
p(\sigma | \lambda) = \frac{2}{\pi\lambda(1+(\sigma/\lambda)^2)}, \quad\quad \sigma \ge 0.
$$
We calculate scaling equivalent to the PC prior scaling for future use:
(lambda = 3*sdres/tan(pi*0.99/2))
INLA works with the precision and the calculation requires the log density of the precision, $\tau$, which is $$ \log p(\tau | \lambda) = -\frac{3}{2}\log\tau - \log (\pi\lambda) - \log (1+1/(\tau\lambda^2)) $$
We then define the prior as:
halfcauchy = "expression: lambda = 0.022; precision = exp(log_precision); logdens = -1.5*log_precision-log(pi*lambda)-log(1+1/(precision*lambda^2)); log_jacobian = log_precision; return(logdens+log_jacobian);"
And use this prior with
hcprior = list(prec = list(prior = halfcauchy)) formula <- nitrogen ~ f(site, model="iid", hyper = hcprior) hmod <- inla(formula, family="gaussian", data=reeds) bri.hyperpar.summary(hmod)
Plot the posterior distributions:
bri.hyperpar.plot(hmod)
The result is quite similar to the penalized complexity prior.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.