Scaling IGMRF-models in R-INLA
Sigrunn Holbek Sørbye
2014

Introduction
Intrinsic Gaussian Markov random fields (IGMRFs) are widely used as priors to model latent dependency
structures in R-INLA. Examples of this type of models include the rw1, rw2, besag, besag2, bym and therw2d models, see http://www.r-inla.org for further descriptions. All of these models can be described interms of having an improper Gaussian density with a scaled precision matrix that reflects the neighbourhood
structure of the model. The scaling is represented as a random precision parameter, and applying thesemodels the user has to define a hyperprior for this parameter. This is a rather delicate issue and such
hyperpriors are often selected in a rather ad-hoc and subjective manner. One problem is that a specific fixed
hyperprior for the precision parameter does not have the same interpretation for different IGMRFs.
The aim of this tutorial is to introduce and demonstrate a recently implemented option in R-INLA named
scale.model, which can be added when specifying the inla formula-call for the different IGMRFs. For
example, we might specify the formula as
formula = y ~f(idx, model = "rw2", scale.model = TRUE, hyper=..)By defining scale.model = TRUE, the rw2-model is scaled to have a generalized variance equal to 1, see
Section 2 for further explanations. The new option is also available using control.group,f(idx, model = "..", group = g, control.group = list(model = "rw2", scale.model = TRUE))
which is relevant for the rw1, rw2 and the besag-models. Note that scale.model = FALSE by default for allmodels to make the option backwards-compatible.
The purpose of scaling the models is to ensure that a fixed hyperprior for the precision parameter has asimilar interpretation for different types of IGMRFs. In general, IGMRFs can be seen as models that reflect
local deviation from an unspecified overall level. More specifically, the models considered here penalize local
deviation from a constant level (rw1, besag, bym), a line (rw2) or a plane (rw2d). A priori, the hyperpriorchosen for the precision parameter will influence how large we allow this local deviation to be. It seems
reasonable to select hyperpriors such that different random effects in a regression model are assumed tohave a similar range for this deviation. By scaling the models, this is achieved by assigning the same fixed
hyperprior to the precisions of all IGMRFs in the regression model. If we do this using unscaled IGMRFs(as often done in the literature), the influence using different hyperprior choices is less controllable as thevariances of the IGMRFs differ.
The hyperpriors for scaled models will also be invariant to the shape and size of the graph for a specificIGMRF and to different scalings of covariates. By mapping the random precision to the marginal variance ofan IGMRF, we also get an alternative interpretation of different parameter choices for a hyperprior, as thesecan be viewed in terms of the local deviation, see Section 2 for further details.
The new option is based on ideas introduced in Sørbye and Rue (2014) and these ideas have also been appliedin Papoila et al. (2014). Note that the implementation in Sørbye and Rue (2014) is slightly different as there
the hyperprior is assigned to scaled precision parameters, while the current option scales the IGMRFs directly.
The results are the same, but the new option makes the ideas in Sørbye and Rue (2014) very easy to apply.

1

We demonstrate this in Section 3, analysing two semiparametric regression models. Concluding remarks are
given in Section 4.

Scaling the IGMRFs
An IGMRF can be defined as a vector x = (x1 , . . . , xn ), having an improper Gaussian density with a sparse
precision matrix which is not of full rank. In Rue and Held (2005), the order of an IGMRF is defined as therank deficiency of its precision matrix. This implies that a zero-mean IGMRF of kth order is specified by
1π(x) = (2π)−(n−k)/2 (|Q|∗ )1/2 exp{− xT Qx},2
where the precision matrix Q is assumed to be semi-positive definite. |Q|∗ denotes the generalized determinantequal to the product of the n − k non-zero eigenvalues of Q. The marginal variances of x are given by thediagonal elements of the generalized inverse matrix Q− , that isσ 2 (xi ) = Q−ii ,i = 1, . . . , n.

Different IGMRFs can be described by expressing the precision matrix as Q = τ R. Here, τ denotes the random
precision parameter while the matrix R reflects the specific neighbourhood structure of the model which canbe represented as a labelled graph with nodes and edges. This implies that the marginal variances/standard
deviations of the IGMRF will depend on both the structure and size of R. For example, we can consider themarginal standard deviations of the rw1 and rw2-models when the precision is fixed as τ = 1. Due to thedifferent structure matrices of these two models, the levels for the marginal standard deviations are quitedifferent, see Figure 1. This illustrates that it seems unreasonable to assign the exact same hyperprior for theprecision parameters of these two models when these parameters are random.
Notice that for a fixed precision τ , the marginal variances of the components of x can be expressed as afunction of τ by
2σ{τ=1} (xi )στ2 (xi ) =, i = 1, . . . , n.
(1)τTo characterize the level of the marginal variances when τ = 1, we define the generalized variance of x as thegeometric mean!
n1X
22log(σ{τ =1} (xi )) .
σGV (x) = exp
n i=12
When scale.model = TRUE, IGMRFs are scaled such that σGV(x) = 1. This implies that the precisionparameters of different models have similar interpretation which makes it reasonable to assign the same fixed
hyperprior to precision parameters of different IGMRFs in a regression model. Note that in Sørbye and Rue2(2014), the generalized variance was referred to as the reference variance σref(x). The hyperprior was then2assigned to the scaled precision τ /σref (x) to achieve the same interpretation for different models.
# Functions to calculate the marginal standard deviations of rw1 and rw2
ss.rw1.sdref = function(values,tau) {n = length(values)
stopifnot(!missing(n))
stopifnot(n > 3)
y = NA

formula = y ~ -1 + f(values, model = "rw1", constr = TRUE,
values = values, diagonal = 1e-10,
hyper = list(prec = list(initial = log(tau), fixed = TRUE)))r = inla(formula, data = data.frame(y, values), family = "gaussian",control.family = list(fixed = TRUE),
control.compute = list(return.marginals = F), verbose = F)2

sd = r$summary.random$values$sd
return (sd)}ss.rw2.sdref = function(values,tau) {n = length(values)
stopifnot(!missing(n))
stopifnot(n > 3)
y = NA
A1 = matrix(1, n, 1)
A2 = matrix(0, n, 1)
for(i in 1:n) A2[i, ]= i
A = rbind(c(A1), c(A2))
e = rep(0, 2)
extraconstr = list(A = A, e= e)
formula = y ~ -1 + f( values, model="rw2", constr = FALSE,
values = values, diagonal = 1e-10, extraconstr = extraconstr,
hyper = list(prec = list(initial = log(tau), fixed = TRUE)))r = inla(formula, data = data.frame(y, values), family = "gaussian",control.family = list(fixed = TRUE),
control.compute = list(return.marginals = F), verbose=F)
sd = r$summary.random$values$sd
return (sd)}

8060
2040

Marginal standard deviation

5.5
5.0
4.5
4.0
3.5
3.0

Marginal standard deviation# The marginal standard deviations for n nodes
n = 100
sd.rw1 = ss.rw1.sdref(1:n, 1)
sd.rw2 = ss.rw2.sdref(1:n, 1)020406080

1000Node i

20406080

100

Node iFigure 1: The marginal standard deviations σ(xi ) for the rw1-model (left) and the rw2-model (right) calculatedusing a fixed precision τ = 1 and n = 100 nodes.
As already mentioned, IGMRFs penalize local deviation from an unspecified overall level. The marginal
standard deviation of the models, integrating out the random precision, give information on how large weallow this local deviation to be. Following Sørbye and Rue (2014), we apply the identity in (1) also when τ is3

random and define an upper limit U for the marginal standard deviation of a model by 2

σGV (x)P (σ(xi ) > U ) ≈ P τ= α,U2(2)where α is a fixed small probability. For scaled models, this implies that the parameter choices for thehyperprior of τ has the same interpretation in terms of U , for all models x.
In R-INLA, precision parameters are assigned Gamma distributions by default (but any other distribution
can be specified by the user). Assume that τ ∼ Gamma(a, b), where a denotes the shape parameter and b
the inverse-scale parameter. Then bτ ∼ Gamma(a, 1) and (2) implies that
2bσGV(x)= F −1 (α, a, 1),U2(3)where F −1 (·) denotes the inverse cumulative distribution function of the Gamma distribution. For scaled
models, we again notice that the interpretation of a and b stays the same for all x and the upper limit U iseasily calculated as
U = sqrt(b/qgamma(alpha, a, 1))
For example, if we choose the default parameter values in R-INLA, a = 1 and b = 5 · 10−5 , the upper limit
equals U ≈ 0.22 when α = 0.001. Notice that the explicit value of U has to be interpreted with care as thisdepends on the selected value of α. Also, U is not uniquely defined and different combinations of a and b
giving the same values of U , do not in general give the same estimated results.

ExamplesThe new option makes it very easy to compare results using unscaled and scaled models. Here we consider
two examples, analysing semiparametric regression models. Other examples are given in Sørbye and Rue(2014) and Papoila et al. (2014).

Semiparametric regression with different scalings of the covariate
We first consider a simple semiparametric regression model in which n observations are generated byyi = f (xi ) + i ,i = 1, . . . , n.

(4)

The noise terms, i , are assumed independent and normally distributed with mean 0 and a fixed standard
deviation σ. The aim is to estimate the underlying smooth function f (.).
Applying the rw2-model as a prior, the formula-call using the scaled model is defined asformula = y ~f(x, model = "rw2", scale.model = TRUE, hyper=...)

where the hyperprior assigned to the precision parameter is specified in hyper. Here, we just select thedefault hyperprior, that is τ ∼ Gamma(1, 5 · 10−5 ).
Let the unknown true function f (.) be defined asf (x) =x
2πx+ sin() − 0.5,ttx ∈ (0, t).

We generate n = 101 observations based on (4) in which σ = 0.5 and call inla:
result = inla(formula, family = "gaussian", data = data.frame(y, x))
We then plot the estimated function

4result$summary.random$x$mean

0.00.2

0.4

0.6

0.8

1.0

0.5
f(x)

−1.0

−0.5

0.0

0.5
0.0
−1.0

−0.5f(x)0.0−1.0−0.5

f(x)

0.5

both using the unscaled and scaled models for three different intervals of x, see Figure 2.

0204060x80

100

0200x

400

600

800

1000xFigure 2: The true function (black solid line) and the estimated functions using the scaled (blue dashed line)
and unscaled (red dashed/dotted line) rw2-models for different scalings of the covariate.
Using the scaled model, the results using a fixed hyperprior is invariant to different scalings of the covariate
and we get exactly the same estimated curve in all three cases. For the unscaled model, the marginal variancewill decrease/increase as the distance between nodes decreases/increases. By using the same fixed hyperpriorfor the three cases, we implicitly assume very different values for the upper limit U in (3), see Table 1. In theextreme case of assuming a very low value of U , the resulting estimate is just a straight line.
Interval length
U (unscaled model)

t=10.009

t = 100
9.4

t = 1000
294.8

Table 1: The upper limit U as defined in (3) for the unscaled rw2-model, calculated when a = 1, b = 5 · 10−5and α = 0.001.
As noted in Akerkar, Martino, and Rue (2010) and Sørbye and Rue (2014), a given hyperprior for the
precision of unscaled models can be adjusted to account for different scalings of the covariate, adjusting theinverse-scale parameter of the Gamma distribution. However, by scaling the models, we don’t have to think
about this at all as a given fixed prior will give the same estimated curve.

Semiparametric regression with several random effects
The linear predictor of a latent Gaussian model often consists of a combination of random effects modelled
by different IGMRFs. As an example we will take a look at the Munich rental data analysed in Section 4.2.2
in Rue and Held (2005), also given as an example in Volume I at http://www.r-inla.org. The responsevariable in this dataset is rent per square meter (in Euros) for a flat and the covariates include spatial location,floor space, year of construction and various indicator variables. Smooth effects of the covariates are modelledusing the rw2-model, while the spatially structured effect for location is modelled using the besag-model.
Previously, this data has been analysed using a fixed Gamma(1, 0.01) for all the precision parameters in themodel. The inverse-scale parameter of the Gamma distribution can be adjusted to account for different typesof IGMRFs by scaling this parameter with the generalized variance, see Sørbye and Rue (2014), but it iseasier to just use the new option of scale.model equal to TRUE or FALSE.
To simplify notation let
hyper.prec = list(prec = list(param = c(1, 0.01)))
The formula call (see http://www.r-inla.org) is then specified bydata(Munich)g = system.file("demodata/munich.graph", package="INLA")## Note that here we want to have an estimator of the effect of year5

## also the for years where we have no observation, therefore we give a## vector with all possible values assumed by the covariate year, that## is seq(1918,2001)
formula = rent ~ -1 +f(location, model = "besag", graph = g, initial = 1, hyper = hyper.prec,scale.model = option) +f(year, model = "rw2", values = seq(1918, 2001), hyper = hyper.prec,scale.model = option) +f(floor.size, model = "rw2", hyper = hyper.prec, scale.model = option) +Gute.Wohnlage + Beste.Wohnlage + Keine.Wwv + Keine.Zh +Kein.Badkach + Besond.Bad + Gehobene.Kueche +zim1 + zim2 + zim3 + zim4 + zim5 + zim6
Further, the call to inla is given as
mod = inla(formula, data = Munich, control.fixed=list(prec=1))−2−20−120412Figure 3 illustrates the resulting estimated effects of the covariates ‘Year of construction’ and ‘Floor size’
using the scaled and unscaled models. Applying the given hyperprior for the scaled models, we have assumed
that the upper limit in (3) is U ≈ 3.2. This seems to give quite reasonable results, while the estimated curves
are more wiggly using the unscaled models. The results for the spatially structured effects are very similarusing the scaled and unscaled models (results not shown).

1920

1940

1960

1980

2000

50

Year of construction

100

150

Floor size

Figure 3: The estimated effects of the covariates ’Year of construction’ and ’Floor size’ on the rental prizes inMunich, using the scaled (black solid line) and the unscaled (red dashed line) rw2-models.
Using the scaled models, we can now study sensitivity to the parameters of the Gamma-prior simultaneously
for the different random effects modelled by IGMRFs, instead of doing this separately for each of the effects.
For example, we can choose to keep the shape parameter fixed at a = 1, while the inverse-scale parameter b
is calculated to give specific values of U . The estimated functions for the effects of the two covariates ‘Yearof construction’ and ‘Floor size’ are given in Figure 4. Here, the value of U ranges from 0.001 to 30 whichcorresponds to values of b in the range (10−9 , 0.9). We conclude that the results are not very sensitive tochanges in the inverse-scale parameter when a = 1. This also applies to the spatially structured effect which
does not change significantly for these different parameter choices (results not shown). Naturally, we could
also consider other values of a.

Concluding remarks
By scaling all IGMRFs to a unit generalized variance, we avoid that the precision parameters of these models
have different interpretation. This makes it reasonable to assign the same hyperprior to precision parameters ofdifferent IGMRF-models, which greatly simplifies the issue of selecting hyperpriors and performing sensitivity

6510

3

4

U= 0.001
U= 0.1
U= 10U= 30−60−10−4012−200

U= 0.001
U= 0.1
U= 10
U= 30

1920

1940

1960

1980

2000

50

Year of construction

100150Floor size

Figure 4: The estimated effects of the covariates ’Year of construction’ and ’Floor size’ on the rental prizes inMunich, varying the upper limit for the scaled models.
analysis. These ideas were introduced in Sørbye and Rue (2014), in which the precision parameter was scaledby the generalized variance. The new option in R-INLA makes these ideas easy to apply as we don’t have tocalculate the generalized variance for each model. We simply set the scale.model-option equal to TRUE in
the inla formula-call for the relevant terms.
If you are now convinced that all IGMRF-models should in fact be scaled, you can do this in an even simpler
way than already described. By specifying the following global option
inla.setOption(scale.model.default = TRUE)
all the IGMRF-models are scaled automatically. This also includes the situation using the control.groupargument mentioned in the introduction. Hyperpriors still have to be selected with care using scaled models,but the issue of selecting hyperpriors is simplified and we believe that hyperpriors can be chosen in a moremeaningful and controlled way.

ReferencesAkerkar, R., S. Martino, and H. Rue. 2010. “Implementing Approximate Bayesian Inference for Survival
Analysis Using Integrated Nested Laplace Approximations.” Preprint Statistics, Norwegian University ofScience and Technology 1: 1–38.
Papoila, A. L., A. Riebler, A. Amaral-Turkman, R. São-João, C. Ribeiro, C. Geraldes, and A. Miranda. 2014.
“Stomach Cancer Incidence in Southern Portugal 1998-2006: A Spatio-Temporal Analysis.” Biometrical
Journal 56: 403–15.
Rue, H., and L. Held. 2005. Gaussian Markov Random Fields. Chapman & Hall/CRC, Boca Raton.
Sørbye, S. H., and H. Rue. 2014. “Scaling Intrinsic Gaussian Markov Random Field Priors in Spatial
Modelling.” Spatial Statistics 8: 39–51. https://doi.org/http://dx.doi.org/10.1016/j.spasta.2013.06.004.

7