The function VarCorrCI takes a merMod object and returns variance components and 95% confidence intervals.
Modified from: http://rpubs.com/bbolker/varwald and various other places.
library(lme4)
# Fit model
fit = lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
# Using ML
VarCorrCI(fit)
# Using REML
VarCorrCI(fit, REML = T)
| require(lme4) | |
| require(numDeriv) | |
| VarCorrCI = function(fit, CI = TRUE, REML = FALSE, ...) { | |
| object = VarCorr(fit) | |
| if (CI == FALSE & REML == FALSE & isREML(fit)) { | |
| warning("Refitting model with ML", call. = FALSE) | |
| fit = refitML(fit) | |
| } | |
| if (!require("numDeriv")) stop("numDeriv package required") | |
| useSc = attr(object, "useSc") | |
| dd = lme4:::devfun2(fit, useSc = useSc, signames = FALSE) | |
| vdd = as.data.frame(object, order = "lower.tri") | |
| nms = apply(vdd[, 1:3], 1, function(x) paste(na.omit(x), collapse = ".")) | |
| pars = vdd[, "sdcor"] | |
| if(CI == FALSE) { | |
| wci = data.frame(level = nms, vcov = pars) | |
| rownames(wci) = NULL | |
| return(wci) | |
| } else { | |
| npar0 = length(pars) | |
| if (isGLMM(fit)) { | |
| pars = c(pars, fixef(fit)) | |
| } | |
| hh1 = hessian(dd, pars) | |
| vv2 = 2 * solve(hh1) | |
| if (isGLMM(fit)) { | |
| vv2 = vv2[1:npar0, 1:npar0, drop = FALSE] | |
| } | |
| names(pars) = nms | |
| dimnames(vv2) = list(nms, nms) | |
| vhack = list(coefficients = pars, vcov = vv2) | |
| wci = confint.default(vhack) | |
| wci = data.frame(level = nms, vcov = pars, wci) | |
| rownames(wci) = NULL | |
| return(wci) | |
| } | |
| } | |
| confint.default = function (object, parm, level = 0.95, ...) { | |
| cf = coef(object) | |
| pnames = names(cf) | |
| if (missing(parm)) | |
| parm = pnames else | |
| if (is.numeric(parm)) | |
| parm = pnames[parm] | |
| a = (1 - level)/2 | |
| a = c(a, 1 - a) | |
| pct = format.perc(a, 3) | |
| fac = qnorm(a) | |
| ci = array(NA, dim = c(length(parm), 2L), dimnames = list(parm, c("lower", "upper"))) | |
| ses = sqrt(diag(object$vcov))[parm] | |
| ci[] = cf[parm] + ses %o% fac | |
| ci | |
| } | |
| format.perc = function(probs, digits) { | |
| paste(format(100 * probs, trim = TRUE, scientific = FALSE, digits = digits), "%") | |
| } |