brinla
Simple example with INLA
Julian Faraway 11 February 2022
See the introduction for general information. We need the following packages (which you must install first).
library(INLA)
Loading required package: Matrix
Loading required package: foreach
Loading required package: parallel
Loading required package: sp
This is INLA_22.01.25 built 2022-01-25 17:43:11 UTC.
- See www.r-inla.org/contact-us for how to get help.
- To enable PARDISO sparse library; see inla.pardiso()
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)
x
76.581
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))
x
12.786
Our estimate is 12.79 billion years. We can form a 95% confidence interval for Hubble’s constant.
(bci <- confint(lmod))
2.5 % 97.5 %
x 68.379 84.783
We can transform this into a 95% confidence interval for the age of the universe in billions of years:
hubtoage(bci)
2.5 % 97.5 %
x 14.32 11.549
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
mean sd 0.025quant 0.5quant 0.975quant mode kld
x 75.375 3.9976 67.362 75.414 83.166 75.493 5.4991e-05
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')]
$mean
[1] 0
$prec
[1] 0.001
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)
[1] 31.623
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)
mean sd 0.025quant 0.5quant 0.975quant mode kld
x 76.581 4.1854 68.284 76.581 84.868 76.581 5.4812e-05
We get a posterior mean which is the same as the least squares estimate. The 95% credibility interval is [68.3, 84.9] 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)])
mean 0.025quant 0.5quant 0.975quant mode
x 12.786 14.339 12.786 11.537 12.786
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)
The 95% credible interval is [11.54, 14.34] billion years.