Created
November 22, 2018 16:35
-
-
Save aghaynes/8c9decf94acdc64d0107222f1f453803 to your computer and use it in GitHub Desktop.
Compare variances from lmer and INLA for a linear mixed model (random intercept)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# | |
# 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") | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment