Notes

A practical INLA introduction

for the INLA course in Iceland 2016

Haakon Bakka <bakka@[Link]>

Department of Mathematical Sciences, NTNU

Wide time estimate: 4h presentation, 3h exercises.

Goal
Notes

The goal of this presentation is to know how to code hierarchical

Bayesian models in the R-INLA interface.

This is not about how INLA works internally. This is not for
discussing whether the model you fit is sensible.
Good news
Notes

All the theory we will see later is wrapped up in the R-package

INLA which is easy to use.

Getting INLA
Notes

I The web page [Link] contains source-code,

worked-through examples, reports and instructions for
installing the package. An INLA tutorial is in preparation.

I The R-package INLA works on Linux, Windows and Mac and

can be installed within R by
[Link]("INLA",
repos="[Link]

Later, it can be upgraded with

[Link](testing=T)
Which INLA version do you have?
Notes

[Link]()
##
##
## INLA version ............: 0.0-1475838132
## INLA date ...............: Fri 7 Oct 13:02:12 CEST 2016
## INLA hgid ...............: hgid: 1a142f78d436 date: Fri Oct 07 12:55:12 2016 +0200
## INLA-program hgid .......: hgid: 1a142f78d436 date: Fri Oct 07 12:55:12 2016 +0200
## Maintainers .............: Havard Rue <hrue@[Link]>
## : Finn Lindgren <[Link]@[Link]>
## : Daniel Simpson <[Link]@[Link]>
## : Andrea Riebler <[Link]@[Link]>
## : Elias Teixeira Krainski <[Link]@[Link]>
## : Geir-Arne Fuglstad <fulgstad@[Link]>
## Web-page ................: [Link]
## Email support ...........: help@[Link]
## : r-inla-discussion-group@[Link]
## Source-code .............: [Link]/hrue/r-inla

How to use INLA

Notes
Example: Ski flying records
There are essentially four parts to an INLA-program:

1. Data organisation: Make an object to store response,

covariates, . . .
data = [Link](y = y, x = x)
2. Use the formula-notation to specify the model (similar to lm
and glm functions)
formula = y~x
3. Call the inla-program

res = inla(formula, data=data, family="gaussian")

4. Extract posterior information, e.g. for a first overview use
summary(res)
The INLA packagefor
INLA package forRR— How does it work?
Notes

Produces:
Input • Input files Output
• .ini files

Runs the A R object

Data frame, INLA
inla in the class
formula package
program inla

Collect
results

Gianluca Baio ( UCL) Introduction to INLA Bayes 2013, 21 May 2013 45 / 92

Main relevance for you is through top and kill, and keep=T.

What happens in the black box?

Notes

The implementation of the INLA method consists of three parts:

GMRFLib-Library: A library for GMRFs written in C
inla-program: The implementation of INLA written in C
INLA package for R: An R-interface to the inla-program

The first two are not particularly user-friendly. They are used in the
background by the INLA package.
Implementing INLA
Notes
All procedures required to perform INLA need to be carefully
implemented to achieve a good speed; easier to implement a slow
version of INLA.
I The GMRFLib-library
I Basic library written in C, user friendly for programmers
I The inla-program
I Define latent Gaussian models and interface with the
GMRFLib-library
I Avoids the need for C-programming
I Models are defined using .ini-files
I Requires to write input files in a special format
I inla-program write all the results (E/Var/marginals) to files
I The INLA package for R
I R-interface to the inla-program.
I Convert “formula”-statements into “.ini”-file definitions
The first two are not particularly user-friendly. They are used in the
background by the INLA package.

What comes out?

Notes

Call:
"inla(formula = formula, family = \"gaussian\", data = data)"

Time used:
Pre-processing Running inla Post-processing Total
0.0581 0.0161 0.0181 0.0924

Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 137.0288 1.3929 134.2798 137.0288 139.7741 137.0288 0
x 2.1259 0.0526 2.0221 2.1259 2.2295 2.1259 0
...
Data organization
Notes

The responses and covariates are collected in a list or data frame.

Assume response y, covariates x1 and x2, and time index t. Then
they can be organized with

# Option 1
data = list(y = y, x1 = x1, x2 = x2, t = t)

# Option 2
data = [Link](y = y, x1 = x1, x2 = x2, t = t)

Variable formula: specifying the linear predictor

Notes

The model is specified through a formula similar to glm:

formula = y ∼ x1 + x2 + f(t, ...)

I y is the name of the response in the data

I The fixed effects are given i.i.d. Gaussian priors
I The f function specifies random effects (e.g. temporal,
spatial, smooth effect of covariates and Besag model)
I Use -1 if you don’t want an automatic intercept
The inla function
Notes

result = inla(
## This is all that is needed for a basic call
# Description of the model
formula,
# Likelihood
family = "gaussian",
# List or data frame with response, covariates, etc.
data = data,

# check what happens

verbose = TRUE,

## ...
# 'control statements'
# [Link]
)

Likelihood functions
Notes

I "gaussian"
I "T"
I "poisson"
I "nbinomial"
I "binomial"
I "exponential"
I "weibull"
I "coxph"
I See list at [Link]
or

names([Link]()$likelihood)
Posterior inference
Notes

Main functions:
I summary(result)
I plot(result)
I result2 = [Link](result)

Example: Simple linear regression

Notes

. . . such as our ski flying example.

Stage 1: Gaussian likelihood

yi | ηi ∼ N (ηi , σo2 )

Stage 2: Covariates are connected to likelihood by

ηi = β0 + β1 xi

Stage 3: σo2 : variance of observation noise

Example: Simple linear regression
Notes

# Generate data
x = runif(10)
y = 1 + 2*x + rnorm(n = 100, sd = 0.1)

# Run inla
formula = y ~ 1 + x
result = inla(formula,
data = [Link](x = x, y = y),
family = "gaussian")

# Get summary
summary(result)

##
## Call: Notes
## c("inla(formula = formula, family = \"gaussian\", data = [Link](x = x, ", " y = y))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.5021 0.1717 0.0597 0.7335
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.9839 0.0260 0.9327 0.9839 1.0349 0.9839 0
## x 2.0394 0.0392 1.9622 2.0394 2.1165 2.0394 0
##
## The model has no random effects
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 91.72 12.42 68.69 91.24
## 0.975quant mode
## Precision for the Gaussian observations 117.24 90.10
##
## Expected number of effective parameters(std dev): 2.007(0.0011)
## Number of equivalent replicates : 49.81
##
## Marginal log-Likelihood: 67.39
Summary for the fixed effects
Notes

result$[Link]
## mean sd 0.025quant 0.5quant 0.975quant mode
## (Intercept) 0.9838558 0.02598130 0.9327357 0.983855 1.034922 0.9838557
## x 2.0393672 0.03922258 1.9621939 2.039366 2.116458 2.0393673
## kld
## (Intercept) 2.014763e-12
## x 1.732002e-12

Marginal posterior densities

Notes

The marginal posterior densities are stored as a matrices with x-

and y -values
m = result$[Link][[1]]
par(mar=c(5,5,1,1))
plot(m)
●
●
●
●●
●
15

●
● ●●
● ●
● ●
● ●
● ●
● ●
● ●
● ●
● ●
10

● ●
● ●
● ●
y

● ●

● ●
5

● ●

● ●
● ●
● ●
● ● ● ● ● ● ●● ●●● ●●● ●● ● ● ● ● ● ●
0

0.8 0.9 1.0 1.1 1.2

x
Marginal posterior densities
Notes

The rough shape can be interpolated to higher resolution:

par(mar=c(5,5,1,1))
plot([Link](m))
●
●●
●
●
●●
●

15
●●
● ●
●●
●
● ●
●
●
●● ●●
●
● ●
●
●
●
●
● ●
●
●
● ●
●
●
● ●
●
●
● ●
●
●
● ●
●
●
●
● ●
●

[Link](m)$y
●
●
● ●
●
●
● ●
●
●
●
●
● ●
●
●
● ●●
●
● ●

10
●
●
● ●
●
●
●
● ●
●
●
● ●●
●
● ●
●
●
● ●
●
●
● ●
●
●● ●
●
●
● ●
●
●
● ●
●
●
● ●
●
●● ●
●
●
● ●
●
● ●
●
●
● ●
●
●● ●●
5
●
● ●
●
●
● ●
●
●● ●
●
●
● ●
●
●
● ●
●
●● ●●
●
● ●
●
●● ●
●●
●
●
● ●
●
●● ●
●●
●
● ●
●
●●
●
● ●
●
●
●
●
●●
●
● ●
●
●
●●
●
●
●●
●
●●
● ●
●●
●
●●
●
●
●
●
●●
●●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●● ●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●●
●
●●
●
●●
●
●
●●
●
●●
●
●●
●
●●
●
●●
●
●●
●
0

0.8 0.9 1.0 1.1 1.2

[Link](m)$x

Marginal posterior densities

Notes

# Extract quantiles
[Link](0.05, m)
## [1] 0.9409562
# Distribution function
[Link](0.975, m)
## [1] 0.3686119
# Density function
[Link](1, m)
## [1] 12.68609
# Generate realizations
[Link](4, m)
## [1] 1.0411516 1.0310270 1.0137149 0.9876703
Marginal posterior densities
Notes

# Calculate expected value of x and x^2

E = [Link](function(x) c(x,x^2), m)
E
## [1] 0.9838558 0.9686471
# Calculate sd
sqrt(E[2]-E[1]^2)
## [1] 0.02597829
# Compare to estimate
round(result$[Link][,1:2], 3)
## mean sd
## (Intercept) 0.984 0.026
## x 2.039 0.039

Organisation of the returned inla-object

Notes

You find summary information (mean, sd, quantiles, . . . ) in:

## [1] "[Link]" "[Link]"
## [3] "[Link]" "[Link]"
## [5] "[Link]" "[Link]"
## [7] "[Link]" "[Link]"

For example:
result$[Link]
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 0.9838558 0.02598130 0.9327357 0.983855 1.034922 0.9838557 2.014763e-12
## x 2.0393672 0.03922258 1.9621939 2.039366 2.116458 2.0393673 1.732002e-12
Organisation of the returned inla-object
Notes

You find summary information (mean, sd, quantiles, . . . ) in:

## [1] "[Link]" "[Link]"
## [3] "[Link]" "[Link]"
## [5] "[Link]" "[Link]"
## [7] "[Link]" "[Link]"

Each object is thereby a list. Get the marginal for intercept:

head(result$[Link][[1]])
## x y
## [1,] 0.7240424 1.191303e-16
## [2,] 0.7760051 7.128905e-11
## [3,] 0.8279677 1.999550e-06
## [4,] 0.8402018 1.677450e-05
## [5,] 0.8538978 1.590287e-04
## [6,] 0.8539491 1.603295e-04

Further general information

Notes

# formula used
result$.args$formula
## y ~ 1 + x
## NULL
# data used
head(result$.args$data)
## x y
## 1 0.59185320 2.2577575
## 2 0.63805949 2.3618730
## 3 0.87653695 2.8205119
## 4 0.91609803 2.9861276
## 5 0.01756631 0.9930623
## 6 0.32969977 1.3996612
# log-file including information of INLA approximations
result$logfile
Get estimates for variance not precision
Notes

1
Goal: Posterior mean and standard deviation of σ02 = τ0 .

names(result$[Link])
## [1] "Precision for the Gaussian observations"
# get the marginal for the precision
tau0 = result$[Link][[1]]

# Calculate expected value of 1/x and 1/x^2

E = [Link](function(x) c(1/x,(1/x)^2), tau0)

# Calculate sd
mysd = sqrt(E[2] - E[1]^2)

print(c(mean=E[1], sd=mysd))
## mean sd
## 0.011092890 0.001593518

Get also the posterior marginal

Notes

# transform the postrior marginal

sigma2 = [Link](function(x){1/x}, tau0)
head(sigma2)
## x y
## [1,] 0.007201669 3.389633
## [2,] 0.007415462 6.208781
## [3,] 0.007549437 8.852326
## [4,] 0.007648673 11.388941
## [5,] 0.007728114 13.848324
## [6,] 0.007794667 16.246947
# from the marginal also 'z'ummary information can be derived
[Link](sigma2)
## Mean 0.0110905
## Stdev 0.00157556
## Quantile 0.025 0.00840267
## Quantile 0.25 0.009961
## Quantile 0.5 0.0109431
## Quantile 0.75 0.012057
## Quantile 0.975 0.0145918
Add random effects
Notes

f(name, model="...", hyper=...,

constr=FALSE, cyclic=FALSE, ...)

I name – the index of the effect (each f-function needs its own!)
I model – the type of latent model. E.g. "iid", "rw2", "ar1",
"besag", and so on
I hyper – specify the prior on the hyperparameters
I constr – sum-to-zero constraint?
I cyclic – are you cyclic?
I ...

Example: Add random effect

Notes

Add an AR(1) random effect to the linear predictor.

Stage 1:
yi |ηi ∼ N (ηi , σo2 )
Stage 2: Covariates and AR(1) component connected to
likelihood by

ηi = β0 + β1 xi + ai

Stage 3: I σo2 : variance of observation noise

I ρ : dependence in AR(1) process
I σ 2 : variance of the innovations in AR(1) process
Example: Add random effect
Notes

# Generate AR(1) sequence

[Link](580258)
t = 1:100
ar = rep(0,100)
for(i in 2:100)
ar[i] = 0.8*ar[i-1]+rnorm(n = 1, sd = 0.1)

# Generate data with AR(1) component

x = runif(100)
y = 1 + 2*x + ar + rnorm(n = 100, sd = 0.2)

# Run inla
formula = y ~ 1 + x + f(t, model="ar1")
result = inla(formula,
data = [Link](x = x, y = y, t = t),
family = "gaussian")

summary(result)

## Notes
## Call:
## c("inla(formula = formula, family = \"gaussian\", data = [Link](x = x, ", " y = y, t = t))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.6199 0.2186 0.0952 0.9336
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.0023 0.0621 0.8783 1.0026 1.1246 1.0031 0
## x 2.0071 0.0837 1.8427 2.0070 2.1719 2.0069 0
##
## Random effects:
## Name Model
## t AR1 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant
## Precision for the Gaussian observations 24.4008 5.4772 15.4221 23.8007
## Precision for t 67.6345 46.1766 16.1276 55.9859
## Rho for t 0.7214 0.1674 0.2877 0.7618
## 0.975quant mode
## Precision for the Gaussian observations 36.8175 22.6452
## Precision for t 187.6038 37.9032
## Rho for t 0.9286 0.8343
##
## Expected number of effective parameters(std dev): 22.25(12.21)
## Number of equivalent replicates : 4.494
##
## Marginal log-Likelihood: -21.02
The interpretation of NA
Notes

R-INLA uses NA differently than other packages

I NA in the response means no likelihood contribution,

i.e. response is unobserved
I NA in a fixed effect means no contribution to the linear
predictor, i.e. the covariate is set equal to zero
I NA in a random effect f(...) means no contribution to the
linear predictor

Prediction
Notes

The distribution of the linear predictor at an unobserved location

can be computed by specifying the value of the covariate x and the
desired time index t and set y to NA.
# Add new location
x = c(x, 0.3)
t = c(t, 101)
y = c(y, NA)

# Re-compute
[Link] = inla(formula,
data = [Link](x = x, t = t, y = y),
family="gaussian",
[Link] = list([Link] = "grid"),
[Link] = list(config = TRUE),
[Link] = list(compute = TRUE))
Prediction
Notes
m = [Link]$[Link][[101]]
round([Link]$[Link][101,], 3)
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## predictor.101 1.573 0.122 1.321 1.576 1.819 1.587 0
par(mar=c(5,5,0.5,0.5))
plot([Link](m), type="l")

4
3
[Link](m)$y

2
1
0

1.0 1.5 2.0 2.5

[Link](m)$x

Prediction
Notes

Caution: This is not yet the predictive distribution, as the

observation noise is missing.
One way to add is by sampling from the posterior distribution.
n = 100000
x = [Link](n, [Link])

x101 = rep(NA, n)
[Link]= rep(NA, n)
for(i in 1:n){
x101[i] = x[[i]]$latent["Predictor.101",]
[Link][i] = x[[i]]$hyperpar[1]
}
Illustration of samples
Notes

library(MASS)
par(mfrow=c(1,2), mar=c(5,5,1,1))
truehist(x101, xlab="Linear predictor 101")
lines([Link]$[Link][[101]], col=2, lwd=2)
truehist([Link], xlab=expression(sqrt(1/sigma^2)))
lines([Link]$[Link][[1]], col=2, lwd=2)

0.08
3

0.06
2

0.04
1

0.02
0.00
0

1.0 1.4 1.8 2.2 20 30 40 50

Linear predictor 101 1 σ2

Posterior predictive distribution

Notes

predDist = rnorm(n, mean=x101, sd=sqrt(1/[Link]))

par(mar=c(5,5,0.5,0.5))
truehist(predDist, ylim=c(0,5))
lines([Link]$[Link][[101]], col="orange", lwd=2)
5
4
3
2
1
0

0.5 1.0 1.5 2.0 2.5

predDist
Example: Smoothing binary time series
Notes

The data set Tokyo is available in the INLA package and consists of
the number of days in Tokyo with rainfall above 1 mm in
1983–1984.
2 | | | | | | || || || |||| ||| | | || | | || || | |

Number of days with rain

1 |||| ||||| || || | ||| || ||||||| ||||| |||||||||| ||| || ||| ||| | || |||| ||| || ||||||| |||| | || |||||||| || |||| | ||||||||| | |||||||||| | |||| |||| | | ||| | ||

0 ||||||||||||||||||||||||||||||||||||||||||||||||||||| ||| |||| || ||||||||||||||||| ||||| |||||||||||| | | |||| ||| ||||||||||||||||||| | |||||||||||| ||| | | ||||||||| | | |||||||||||||||||||||||||||||||||||||||||| ||||||||||

0 100 200 300

Observations
Notes
Each observation consists of
t: Day of year; t ∈ {1, 2, . . . , 366}
nt : Number of observations for day t in 1983–1984; nt ∈ {1, 2}
yt : Number of days with rain out of nt days for day t;
yt ∈ {0, 1, 2}

data(Tokyo)
head(Tokyo,4)
## y n time
## 1 0 2 1
## 2 0 2 2
## 3 1 2 3
## 4 1 2 4
Tokyo[60,]
## y n time
## 60 0 1 60
Hierarchical model
Notes

Stage 1: We have binomial responses with known nt , but

unknown probabilities

yt ∼ Binomial(nt , pt )

Stage 2: A cyclic second order random walk (CRW2) is

connected to the likelihood by

exp (ηt )
pt = with linear predictor ηt = CRW2t
1 + exp(ηt )

Stage 3: τ : Scale parameter in CRW2 with prior

π(τ ) ∼ Gamma(1, 5 · 10−5 )

Computations
Notes

# Read data
data(Tokyo)

# Specify linear predictor

formula = y ~ -1 + f(time, model="rw2", cyclic=TRUE)

# Run model
result = inla(formula,
family = "binomial",
Ntrials = n,
data = Tokyo)
Marginal posterior of CRW2
Notes

par(mar=c(5,5,1,0.5))
toplot = c("mean", "0.025quant", "0.5quant", "0.975quant")
matplot(result$[Link]$t[, toplot],
lty =c(1,3,2,3), type="l", col=1, ylab="")

0.5
0.0
−2.5 −2.0 −1.5 −1.0 −0.5

0 100 200 300

Marginal posterior of scale parameter

Notes

par(mar=c(5,5,1.5,0.5))
plot([Link](result$[Link][[1]]), xlim=c(0, 10^6),
xlab=expression(tau), ylab="Density", type="l")
4e−05
Density

2e−05
0e+00

0e+00 2e+05 4e+05 6e+05 8e+05 1e+06

τ
Transform to probability
Notes

result = inla(formula, family = "binomial",

Ntrials = n, data = Tokyo,
[Link] = list(compute=TRUE))
par(mar=c(5,5,1.5,0.5))
matplot(result$[Link][, toplot],
lty =c(1,3,2,3), type="l", col=1, ylim=c(0,1), ylab="")

1.0
0.8
0.6
0.4
0.2
0.0

0 100 200 300

Other choices for f-terms

Notes

names([Link]()$latent)
## [1] "linear" "iid" "mec"
## [4] "meb" "rgeneric" "rw1"
## [7] "rw2" "crw2" "seasonal"
## [10] "besag" "besag2" "bym"
## [13] "bym2" "besagproper" "besagproper2"
## [16] "fgn" "ar1" "ar"
## [19] "ou" "generic" "generic0"
## [22] "generic1" "generic2" "generic3"
## [25] "spde" "spde2" "spde3"
## [28] "iid1d" "iid2d" "iid3d"
## [31] "iid4d" "iid5d" "2diid"
## [34] "z" "rw2d" "rw2diid"
## [37] "slm" "matern2d" "copy"
## [40] "clinear" "sigm" "revsigm"
## [43] "log1exp" "logdist"
Add weight to components of a random effect
Notes

formula = y ~ ... + f(idx , weight, model = <MODEL>, ...)

extends the usual

ηi = . . . + fidxi
to
ηi = . . . + weightidxi fidxi

Changing the prior: Internal scale

Notes

I Hyperparameters are represented internally with more

well-behaved transformations, e.g. correlation ρ and precision
τ are internally represented as

θ1 = log(τ )
 
1+ρ
θ2 = log
1−ρ
I The prior must be set on the parameter in internal scale
I Initial values for the mode-search must be set in internal scale
I The functions [Link] and [Link] can be used to map
back and forth.
Changing the prior: Code
Notes

hyper = list(prec = list(prior = "loggamma",

param = c(1, 0.1),
initial = 4,
fixed = FALSE))

formula = y ~ f(idx, model = "iid", hyper = hyper) + ...

# For the iid model, default options can be seen with

[Link]("iid")
We come back to priors in the last lecture block.

EPIL example
Notes

Seizure counts in a randomised trial of anti-convulsant therapy in

epilepsy. From WinBUGS manual.
Patient y1 y2 y3 y4 Trt Base Age
1 5 3 3 3 0 11 31
2 3 5 3 3 0 11 30
3 2 4 0 5 0 6 25
....
59 1 4 3 2 1 12 37

Covariates are treatment (0,1), 8-week baseline seizure counts, and

age in years.
Repeated Poisson counts
Notes

yjk ∼ Poisson(µjk ); j = 1, . . . , 59; k = 1, . . . , 4

log(µjk ) = α0 + α1 log(Basej /4) + α2 Trtj

+α3 Trtj log(Basej /4) + α4 log(Agej )
+α5 V 4 + Indj + βjk

αi ∼ N (0, τα ) τα known (0.001)

Indj ∼ N (0, τInd ) τInd ∼ Gamma(1, 0.01)
βjk ∼ N (0, τβ ) τβ ∼ Gamma(1, 0.01)

Here, V4 is an indicator variable for the 4th visit.

Model specification in INLA

Notes
1 > data ( Epil )
2 > head ( Epil , n =3)
3 y Trt Base Age V4 rand Ind CTrt ClBase4 CV4 ClAge
4 1 5 0 11 31 0 1 1 -0.5254237 -0.75635379 -0.25 0.11420370
5 2 3 0 11 31 0 2 1 -0.5254237 -0.75635379 -0.25 0.11420370
6 3 3 0 11 31 0 3 1 -0.5254237 -0.75635379 -0.25 0.11420370
7 4 3 0 11 31 1 4 1 -0.5254237 -0.75635379 0.75 0.11420370

1 > formula = y ~ ClBase4 * CTrt + ClAge + CV4 +

2 f ( Ind , model = " iid " ,
3 hyper = list ( prec = list ( prior = " loggamma " ,
4 param = c (1 ,0.01) ) ) ) +
5 f ( rand , model = " iid " ,
6 hyper = list ( prec = list ( prior = " loggamma " ,
7 param = c (1 ,0.01) ) )

1 > result = inla ( formula , family = " poisson " , data = Epil ,
2 control . fixed = list ( prec . intercept = 0.001 ,
3 prec = 0.001) )
Bayesian hierarchical models
Notes

INLA can be used with Bayesian hierarchical models where we

model in different stages or levels:

Stage 1: What is the distribution of the responses?

Stage 2: What is the distribution of the underlying unobserved
(latent) components?
Stage 3: What are our prior beliefs about the parameters
controlling the components in the model?

Stage 1
Notes

How is our data (y ) generated from the underlying components (x)

and hyperparameters (θ) in the model:

I Gaussian response? (temperature, rainfall, fish weight ...)

I Count data? (people infected with a disease in each area)
I Point pattern? (locations of trees in a forest)
I Binary data? (yes/no response, binary image)
I Survival data? (recovery time, time to death)
(It is also important how data are collected!)

This information is placed into our likelihood π(y |x, θ)

Stage 2
Notes

The underlying unobserved components x are called latent

components and can be:

I Fixed effects for covariates

I Unstructured random effects (individual effects, group effects)
I Structured random effects (AR(1), regional effects, . . . )

These are linked to the responses in the likelihood through linear

predictors.

Stage 3
Notes

The likelihood and the latent model typically have hyperparameters

that control their behavior. The hyperparameters θ can include:

Examples likelihood:
I Variance of observation noise
I Dispersion parameter in the negative binomial model
I Probability of a zero (zero-inflated models)

Examples latent model:

I Variance of unstructured effects
I Correlation of multivariate effects
I Range and variance of spatial effects
I Autocorrelation parameter
Comparing results with MCMC
Notes

I When comparing the results of R-INLA with MCMC, it is

important to use the same model. That means, same data,
same priors, same constraints on parameters, intercept
included or not, . . . .
I Here we have compared the results with those obtained using
JAGS via the rjags package

Intercept, 0.125 minutes Age

Notes
1.5
5

1.0
4
Density

Density
3

0.5
2
1

0.0
0

1.4 1.5 1.6 1.7 1.8 1.9 −0.5 0.0 0.5 1.0 1.5

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0
Density

Density

1.0
0.5

0.5
0.0

0.0

0.5 1.0 1.5 2.0 2.5 1.5 2.0 2.5

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

 Intercept, 0.25 minutes Age
Notes

1.2
5
4

0.8
Density

Density
3

0.4
2
1

0.0
0
1.4 1.5 1.6 1.7 1.8 1.9 −0.5 0.0 0.5 1.0 1.5

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0
Density

Density

1.0
0.5

0.5
0.0

0.0
0.5 1.0 1.5 2.0 2.5 1.5 2.0 2.5 3.0

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

Intercept, 0.5 minutes Age

Notes
1.2
5
4

0.8
Density

Density
3

0.4
2
1

0.0
0

1.3 1.4 1.5 1.6 1.7 1.8 1.9 −0.5 0.0 0.5 1.0 1.5 2.0

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0
Density

Density

1.0
0.5

0.5
0.0

0.0

0.5 1.0 1.5 2.0 2.5 1.5 2.0 2.5 3.0

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

 Intercept, 1 minutes Age
Notes

1.2
5
4

0.8
Density

Density
3
2

0.4
1

0.0
0
1.3 1.4 1.5 1.6 1.7 1.8 1.9 −0.5 0.0 0.5 1.0 1.5 2.0

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0
Density

Density

1.0
0.5

0.5
0.0

0.0
0.5 1.0 1.5 2.0 2.5 1.5 2.0 2.5 3.0

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

Intercept, 2 minutes Age

Notes
1.2
5
4

0.8
Density

Density
3
2

0.4
1

0.0
0

1.3 1.4 1.5 1.6 1.7 1.8 1.9 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0
Density

Density

1.0
0.5

0.5
0.0

0.0

0.5 1.0 1.5 2.0 2.5 1.5 2.0 2.5 3.0

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

 Intercept, 4 minutes Age
Notes

1.2
5
4

0.8
Density

Density
3
2

0.4
1

0.0
0
1.3 1.4 1.5 1.6 1.7 1.8 1.9 −1.0 0.0 0.5 1.0 1.5 2.0

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0

1.0
Density

Density
0.5

0.5
0.0

0.0
0.5 1.0 1.5 2.0 2.5 3.0 1.5 2.0 2.5 3.0

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

Intercept, 8 minutes Age

Notes
0.0 0.2 0.4 0.6 0.8 1.0 1.2
5
4
Density

Density
3
2
1
0

1.3 1.4 1.5 1.6 1.7 1.8 1.9 −1.0 0.0 0.5 1.0 1.5 2.0

a0 [Link]

log([Link]) log([Link])
1.5

1.5
1.0

1.0
Density

Density
0.5

0.5
0.0

0.0

0.5 1.0 1.5 2.0 2.5 3.0 1.5 2.0 2.5 3.0

log(tau.b1) log(tau.b)

Running time of INLA < 0.5 seconds

Control statements
Notes

[Link] statements control computations

I [Link]
I prec: Default precision for all fixed effects except the
intercept. [Link]: Precision for intercept (Default:
0.0)
I [Link]
I compute: Compute posterior marginals of linear predictors
I [Link]
I dic, mlik, cpo: Compute measures of fit?
I config: Save internal GMRF approximations? (needed to use
[Link]())
I [Link]
I strategy and [Link] contain useful advanced
features
I There are various others as well; see help.

Summary: How to fit a model in R-INLA

Notes

1. Find the likelihood(s) (code:family) and link function(s)

(code:link) corresponding to the first level of you model
You may need family=list(...), link=list(...), data = list(y =
list(...))
2. Set up the formula for your linear predictor, including all the
random effects
You may need spatial models, group, copy, generic, cyclic,
[Link], ..., or rgeneric!
3. Set the priors (possibly also initial values) for your
hyper-parameters
You may need the joint specification in the inla() call
4. Summarise and plot what you want from the posterior
You may need lincombs (linear combinations)
Thank you for your attention!
Notes

If you have any doubts or questions, please write us at

[Link]

Exercise
Notes

Exercise set 1: Problem 2 and Problem 3