User defined integration points
Haavard Rue (hrue@r-inla.org)
May 2017 (updated April 2019 and May 2023)
Introduction
This short note describe a new option that allow the user to use user-defined integration points (or “design”points), instead of the default ones. The relevant integration in INLA does
Zf (x|θ, y) π(θ|y) dθ = f (x|y)
where π(θ|y) is the approximated posterior marginal for the hyperparameters, and where f (x|θ, y) is theapproximated marginal for x for that configuration. The output of this integral is the posterior marginalf (x|y). In practice, we use a discrete set of integration points for θ, and corresponding weights w, to getX
f (x|θi , y)wi π(θi |y)f (x|y) ≈ifor which we require wi ≥ 0 andP

i wi = 1. Usually, the integration is done in a standardised scale,z = A(θ − γ)i.e. with respect to π(z|y). Here, γ is the mode of π(θ|y) and the matrix A is the negative square root of theapproximated covariance matrix for θ|y at the mode.
The relevant options areopts = control.inla(int.strategy = "user", int.design = Design)where Design is a matrix with the integration points and the integration weights. The jth row of Design
consists of the values θj = (θ1j , . . . , θmj ), and the integration weight for this configuration, wj . The values
are in the θ-scale, meaning that you have to know exactly what you are doing, including knowing the orderingof the hyperparameters.
Another version, is to define the points in the standardised scale z. To do this, use
opts = control.inla(int.strategy = "user.std", int.design = Design)instead. The meaning of Design is unchanged, except that these can be given in standardised coordinates.
This version is more relevant if you want to implement a generic new integration design instead of the ones
already provided.

First example
In this artificial example, we want to compute the change of the marginal variance of one component, x1 , ofa hidden AR(1) process, with respect to lag one correlation ρ. So we want to compute∂SD(x1 |y)∂ρ1for a fixed value of ρ = ρ0 . We have to compute a numerical approximation, using finite difference. While
doing this, it is a good idea to keep the design fixed, to avoid introducing an error for changing that part aswell.
Let us first setup the experimentn = 100
rho = 0.9
x = scale(arima.sim(n, model = list(ar = rho)))y = x + rnorm(n, sd = 0.1)
this gives the following

0−2

−1

value

1

2

plot(y, xlab = "time", ylab = "value")lines(x, lwd=2)

020

40

60
time

To compute the derivative, we dorho.0 = rho
to.theta = inla.models()$latent$ar1$hyper$theta2$to.theta
rho.0.internal = to.theta(rho.0)
r = inla(y ~ -1 + f(time, model="ar1",
hyper = list(
theta1 = list(prior = "loggamma",
param = c(1,1)),
theta2 = list(initial = rho.0.internal,fixed=TRUE))),
control.inla = list(int.strategy = "grid"),

280

100

data = data.frame(y, time = 1:n))sd.0 = r$summary.random$time[1,"sd"]print(sd.0)## [1] 0.01222806
We will now change ρ a little, while we keep the same integration points,summary(r)#### Call:##c("inla.core(formula = formula, family = family, contrasts = contrasts,##", " data = data, quantiles = quantiles, E = E, offset = offset, ", "##scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,##", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =##verbose, ", " lincomb = lincomb, selection = selection, control.compute##= control.compute, ", " control.predictor = control.predictor,##
control.family = control.family, ", " control.inla = control.inla,##control.fixed = control.fixed, ", " control.mode = control.mode,##control.expert = control.expert, ", " control.hazard = control.hazard,##control.lincomb = control.lincomb, ", " control.update =##control.update, control.lp.scale = control.lp.scale, ", "##control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,##", " inla.call = inla.call, inla.arg = inla.arg, num.threads =##num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,##working.directory = working.directory, ", " silent = silent, inla.mode##= inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =##.parent.frame)")## Time used:##
Pre = 0.606, Running = 0.338, Post = 0.0274, Total = 0.972## Random effects:##NameModel##time AR1 model#### Model hyperparameters:##mean
sd 0.025quant 0.5quant## Precision for the Gaussian observations 2.17e+04 2.35e+04
1588.695 1.44e+04## Precision for time
8.48e-01 1.19e-01
0.636 8.41e-01##0.975quant mode## Precision for the Gaussian observations
84027.47 85.88## Precision for time
1.10 0.83#### Marginal log-Likelihood: -70.33## is computed## Posterior summaries for the linear predictor and the fitted values are computed## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')nm <- nrow(r$summary.hyperpar)
Design = as.matrix(cbind(r$joint.hyper[, seq_len(nm)], 1))head(Design)##Log precision for the Gaussian observations Log precision for time 1## [1,]
5.613519-0.44610162 1## [2,]
5.613578-0.27923996 1
3

## [3,]## [4,]## [5,]## [6,]

5.613617
5.613654
5.613710
7.500788-0.16799886 1-0.06122839 1
0.09892731 1-0.66924672 1

where the last column is the (un-normalised) integration weights. The integration weights should depend on
the area associated with every point, but for simplicity we just set the weights to be a constant.
Design has dimension 35, 3. We call inla() again reusing the previous found modeh.rho = 0.01
rho.1.internal = to.theta(rho.0 + h.rho)
rr = inla(y ~ -1 + f(time, model="ar1",
hyper = list(
theta1 = list(prior = "loggamma",
param = c(1,1)),
theta2 = list(initial = rho.1.internal,fixed=TRUE))),
control.mode = list(result = r, restart=FALSE),
data = data.frame(y, time = 1:n),control.inla = list(int.strategy = "user",int.design = Design))
sd.1 = rr$summary.random$time[1,"sd"]print(sd.1)## [1] 0.01050915
and then our estimate of the derivative isderiv.1 = (sd.1 - sd.0) / h.rho
print(deriv.1)## [1] -0.1718913
PS: In the logfile of the inla()-call, the configurations are shown in the z scale even for int.strategy="user".

Second example
There is also another (experimental) option, that iscontrol.inla = list(int.stategy = "user.expert")
for which the weights includes π(θi |y). The following example show how to use it, fitting the same model inthree different ways.
n = 50x = rnorm(n)y = 1 + x + rnorm(n, sd = 0.2)
param = c(1, 0.04)
dz = 0.1
r.std = inla(y ~ 1 + x, data = data.frame(y, x),
control.inla = list(int.strategy = "grid",dz = dz,diff.logdens = 8),
control.family = list(hyper = list(prec = list(prior = "loggamma",

4param = param))))s = r.std$internal.summary.hyperpar[1,"sd"]m = r.std$internal.summary.hyperpar[1,"mean"]theta = m + s*seq(-4, 4, by = dz)
weight = dnorm(theta, mean = m, sd = s)r = rep(list(list()), length(theta))
for(k in seq_along(r)) {r[[k]] = inla(y ~ 1 + x,control.family = list(hyper = list(prec = list(initial = theta[k],fixed=TRUE))),data = data.frame(y, x))}r.merge = inla.merge(r, prob = weight)## Warning in inla.merge(r, prob = weight): This function is experimental.
## Warning in inla.merge(r, prob = weight): Merging 'cpo' and 'pit'-results## are/can be, approximate only
r.design = inla(y ~ 1 + x,
data = data.frame(y, x),control.family = list(hyper = list(prec = list(## the prior here does not really matter, as we will override## it with the user.expert in any case.
prior = "pc.prec",param = c(1, 0.01)))),
control.inla = list(int.strategy = "user.expert",
int.design = cbind(theta, weight)))for(k in 1:2) {plot(inla.smarginal(r.std$marginals.fixed[[k]]),lwd = 2, lty = 1, type = "l",xlim = inla.qmarginal(c(0.0001, 0.9999), r.std$marginals.fixed[[k]]))lines(inla.smarginal(r.design$marginals.fixed[[k]]),lwd = 2, col = "blue", lty = 1)
lines(inla.smarginal(r.merge$marginals.fixed[[k]]),lwd = 2, col = "yellow", lty = 1)
}

5

1210
8
6
4
2

inla.smarginal(r.std$marginals.fixed[[k]])$y

0

0.90

0.95

1.00

1.05

inla.smarginal(r.std$marginals.fixed[[k]])$x

6

1.10

1210
8
6
4
2

inla.smarginal(r.std$marginals.fixed[[k]])$y

0

0.90

0.95

1.00

1.05

inla.smarginal(r.std$marginals.fixed[[k]])$x

7

1.10