Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# aghaynes/lmer_vs_inla_lmm.R

Created Nov 22, 2018
Compare variances from lmer and INLA for a linear mixed model (random intercept)
 # # 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` hi <- imod\$summary.hyperpar\$`0.975quant` 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` hi <- imod\$summary.hyperpar\$`0.975quant` 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")
to join this conversation on GitHub. Already have an account? Sign in to comment