Compare variances from lmer and INLA for a linear mixed model (random intercept)

lmer_vs_inla_lmm.R

# 
# Compare lmer and inla for LMM
#  largely taken from Spatial and spatio-temporal bayesian models with R-INLA (Blangiardo & Cameletti, 2015), section 5.4.2
# 

m <- 10000 # N obs
set.seed(1234)
x <- rnorm(m)
group <- sample(seq(1, 100), size = m, replace = TRUE)

# simulate random intercept
tau.ri <- .25
1/sqrt(tau.ri) #SD
set.seed(4567)
v <- rnorm(length(unique(group)), 0, sqrt(1/tau.ri))
# assign random intercept to individuals
vj <- v[group]
# simulate y
tau <- 3
1/sqrt(tau) #SD
set.seed(8910)
b0 <- 5
beta1 <- 2
y <- rnorm(m, b0 + beta1*x + vj, 1/sqrt(tau))

library(lme4)
mod <- lmer(y ~ x + (1|group))
summary(mod)
vc <- VarCorr(mod)


library(INLA)
form <- y ~ x + f(group, model = "iid", param = c(1, 5e-5))
imod <- inla(form, family = "gaussian", 
             data = data.frame(y = y, x = x, group = group))
summary(imod)


cbind(truth = c(tau, tau.ri), lmer = 1/c(attr(vc, "sc")^2, unlist(vc)), inla = imod$summary.hyperpar$`0.5quant`)


plot(imod, 
     plot.fixed.effects = F,
     plot.lincomb = F,
     plot.random.effects = F,
     plot.predictor = F, 
     plot.prior = TRUE)


plot(imod$marginals.hyperpar$`Precision for the Gaussian observations`, type = "l")
# the equivalent of this for lmer is not easy to get at all. One would have to profile the deviance function (see http://lme4.r-forge.r-project.org/slides/2009-07-16-Munich/Precision-4.pdf)
lo <- imod$summary.hyperpar$`0.025quant`[1]
hi <- imod$summary.hyperpar$`0.975quant`[1]
dd <- imod$marginals.hyperpar$`Precision for the Gaussian observations`
dd <- dd[dd[,1] >= lo & dd[,1] <= hi, ]
dd <- rbind(dd, c(hi, 0))
dd <- rbind(dd, c(lo, 0))
polygon(dd, col = "blue")


plot(imod$marginals.hyperpar$`Precision for group`, type = "l")
lo <- imod$summary.hyperpar$`0.025quant`[2]
hi <- imod$summary.hyperpar$`0.975quant`[2]
dd <- imod$marginals.hyperpar$`Precision for group`
dd <- dd[dd[,1] >= lo & dd[,1] <= hi, ]
dd <- rbind(dd, c(hi, 0))
dd <- rbind(dd, c(lo, 0))
polygon(dd, col = "blue")