Skip to content

Instantly share code, notes, and snippets.

@jslefche
Last active March 21, 2020 15:14
Show Gist options
  • Save jslefche/12e11380c5efa91bfedf81751afac245 to your computer and use it in GitHub Desktop.
Save jslefche/12e11380c5efa91bfedf81751afac245 to your computer and use it in GitHub Desktop.
Hierarchical variance partitioning

Hierarchical variance partitioning used linear mixed effects models

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.

EXAMPLE

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), "%")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment