In julianfaraway/brinla: Bayesian Regression with INLA

knitr::opts_chunk$set(echo = TRUE,comment=NA, warning=FALSE, fig.path="figs/",dev = 'svglite')

See the introduction for general information. We need the following packages (which you must install first).

library(INLA)
library(brinla)

Data and Model

The data for the example comes from the Hubble space telescope. It can also be found in the gamair package. We have the relative velocity and distance for 24 galaxies. Read background info Load and plot the data:

plot(y ~ x, xlab="Distance(Mpc)", ylab="Velocity(km/s)", data=hubble)

The velocity y is related to the distance x by y=βx

where β is Hubble's constant. The age of universe is approximated by the inverse of Hubble's constant.

Least squares

We can estimate beta using the standard linear model:

lmod <- lm(y~x-1, data=hubble)
coef(lmod)

We have 60 seconds, 60 minutes, 24 hours and 365.25 days in a year and one megaparsec (Mpc) is 3.09e19 km. Here is a function to transform Hubble's constant into the age of the universe in billions of years. We apply it our estimate:

hubtoage <- function(x) 3.09e19/(x*60^2*24*365.25*1e9)
hubtoage(coef(lmod))

Our estimate is r round(hubtoage(coef(lmod)),2) billion years. We can form a 95% confidence interval for Hubble's constant.

(bci <- confint(lmod))

We can transform this into a 95% confidence interval for the age of the universe in billions of years:

hubtoage(bci)

The upper and lower bound come out in the reverse order because of the inversion of Hubble's constant.

INLA with the default prior

We can fit the same model using the default setting in INLA:

imod <- inla(y ~ x -1, family="gaussian", data=hubble)
imod$summary.fixed

The resulting posterior mean is different from the linear model estimate of beta. Since we are using a different method, we should not expect it to be the same. Even so, there is a reason for this difference. The regression parameters of an LGM have, by definition, a Gaussian prior. The default mean and precision are:

inla.set.control.fixed.default()[c('mean','prec')]

We see that the default prior on beta is normal with mean zero and precision 0.001. The precision is the inverse of the variance. We convert this to SD:

sqrt(1/0.001)

We can see that the linear model fit for beta is somewhat more than two standard deviations from prior mean. So in this case, the default prior is actually quite informative. We need to make an adjustment by setting the precision to a much smaller value.

INLA with a weakly informative prior

We use a very small precision:

imod <- inla(y ~ x -1, family="gaussian", control.fixed=list(prec=1e-9), data=hubble)
(ibci <- imod$summary.fixed)

We get a posterior mean which is the same as the least squares estimate. The 95% credibility interval is [r round(ibci[c(3,5)],1)] which is comparable but also a bit different from the 95% confidence interval. Quite apart from the numerical differences, the two intervals are conceptually different.

We can make a plot of the posterior density of beta (with the 95% credible interval added)

x <- seq(60, 100, length.out = 100)
plot(imod$marginals.fixed$x, type = "l", xlab = "beta", ylab = "density", 
    xlim = c(60, 100))
abline(v = ibci[c(3, 5)], lty = 2)

We can convert to statements about the age of the universe in terms of billions of years:

hubtoage(ibci[c(1,3,4,5,6)])

This can also be plotted:

ageden <- inla.tmarginal(hubtoage, imod$marginals.fixed$x)
plot(ageden, type = "l", xlab = "Age in billions of years", ylab = "density")
abline(v = hubtoage(ibci[c(3, 5)]), lty = 2)