Skip to content

Instantly share code, notes, and snippets.

@jslefche
Last active January 19, 2023 05:21
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jslefche/00be84936c0d7c3299d2d205a6d1dedc to your computer and use it in GitHub Desktop.
Save jslefche/00be84936c0d7c3299d2d205a6d1dedc to your computer and use it in GitHub Desktop.
Hierarchically partition variance in GLMM

Compute semi-variance (partial R[2]) for GLMMs fit used lme4.

library(lme4)

example <- data.frame(
  y = rnorm(100),
  x = rnorm(100),
  nested1 = letters[1:20],
  nested2 = rep(letters[1:5], each = 4),
  crossed = letters[4:1]
)

model <- lmer(y ~ x + (1 | nested1/nested2) + (1 | crossed), example)

part_GLMM(model)
#' Function to hierarchically partition variance in GLMM
#'
#' @param model a GLMM fitted using lmer in the lme4 package
#' @param fixed whether fixed effects variance should be partitioned (using partR2)
#'
#' @required piecewiseSEM lme4 partR2
#'
part_GLMM <- function(model, fixed = FALSE) {
# Get variance components of random effects
randvar <- as.data.frame(lme4::VarCorr(model))[, c("grp", "vcov")]
# grp = the different levels of the random effects
# vcov = variance components (aka, the random variance apportioned to each level of the random effects)
# Remove residual variance from table & store for later
residvar <- randvar[nrow(randvar), , drop = F]
names(residvar)[2] <- "variance"
randvar <- randvar[-nrow(randvar), , drop = F]
# Scale by total variance
randvar$vcov <- randvar$vcov / (
# random effects variance
sum(randvar$vcov) +
# residual variance
residvar$variance +
# fixed effects variance
var(as.vector(lme4::fixef(model) %*% t(model.matrix(model)))) )
names(randvar)[2] <- "Partial_R2"
# Get fixed effects variance
if(fixed == TRUE) fixedvar <- partR2::partR2(model, partvars = names(fixef(model))[-1])
# Return output
l <- list(
R2 = piecewiseSEM::rsquared(model),
RandomVariance = randvar,
ResidualVariance = residvar)
if(fixed == TRUE) l <- append(l, list(FixedVariance = fixedvar))
return(l)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment