In julianfaraway/brinla: Bayesian Regression with INLA

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 more about INLA. This tutorial paper is directed towards spatial problems but same methods can be profitably applied on one-dimensional smoothing problems. See also [@lind:rue:jss:15] for details. The construction is detailed in our book.

Load the packages (you may need to install the brinla package):

library(ggplot2)
library(INLA)
library(brinla)

Data

We use the fossil example from [@bral:gsa:97] and used by [@chaud:jasa:99]. We have the ratio of strontium isotopes found in fossil shells in the mid-Cretaceous period from about 90 to 125 million years ago. We rescale the response as in the SiZer paper.

data(fossil, package="brinla")
fossil$sr <- (fossil$sr-0.7)*100

Plot the data:

plot(sr ~ age, fossil, xlab="Age",ylab="Strontium Ratio")

GP fitting

The default fit uses priors based on the SD of the response and the range of the predictor to motivate sensible priors.

gpmod = bri.gpr(fossil$age, fossil$sr)

We can plot the resulting fit and 95% credible bands

plot(sr ~ age, fossil, xlab="Age",ylab="Strontium Ratio")
lines(gpmod$xout, gpmod$mean)
lines(gpmod$xout, gpmod$lcb, lty=2)
lines(gpmod$xout, gpmod$ucb, lty=2)

Basis functions

The default number of spline basis functions is 25. Let's see what happens if you increase this to 100:

gpmod = bri.gpr(fossil$age, fossil$sr, nbasis = 100)

It makes very little difference although it does take longer to compute. This demonstrates that you just need enough splines for the flexibility you want. If you use more, you won't get a rougher fit.

Spline degree

We can decrease the degree of the splines (default is 2 corresponding to cubic splines):

gpmod = bri.gpr(fossil$age, fossil$sr, degree=0)

We get a piecewise constant fit (the diagonal parts are just from interpolating the fit - could use a better grid to avoid this)

GP kernel shape

We can change the shape of the GP kernel where default is alpha = 2

gpmod = bri.gpr(fossil$age, fossil$sr, alpha = 1)

The kernel is less smooth resulting in a rougher fit.

Prior means on sigma and range

We can set the prior mean on the error to be much larger than the default (which is the SD of the response - already quite large:

gpmod = bri.gpr(fossil$age, fossil$sr, sigma0 = 10*sd(fossil$sr))

Doesn't have much affect on the outcome. Nice to have a robust choice of prior.

We can set the prior mean on the error to be much smaller than the default (which is the SD of the response)

gpmod = bri.gpr(fossil$age, fossil$sr, sigma0 = 0.1*sd(fossil$sr))

Again quite robust to this choice.

We can also experiment with the range. The default is one quarter of the range of the predictor. Let's make it equal to the range:

gpmod = bri.gpr(fossil$age, fossil$sr, rho0 = max(fossil$age) - min(fossil$age))

Again makes very little difference.

Penalized complexity priors

We can also use penalized complexity priors. We set a high value for sigma where we judge only a 5% chance it is more than this and a high value for the range where we judge only a 5% chance that it is less than this.

highsig = 10*sd(fossil$sr)
lowrho = 0.05*(max(fossil$age) - min(fossil$age))
gpmod = bri.gpr(fossil$age, fossil$sr, pcprior = c(highsig, lowrho))

Package versions

sessionInfo()

References

julianfaraway/brinla documentation built on April 6, 2023, 2:33 p.m.