Skip to content
R-squared for generalized linear mixed effects models
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Type Name Latest commit message Commit time
Failed to load latest commit information.
.gitignore Casallas changes Jan 27, 2014 Added link to new function in piecewiseSEM package Aug 18, 2015
rsquaredglmm.R Calculate sigma for nested lme Jul 10, 2014

R-squared for generalized linear mixed-effects models

This function has been completely rewritten and included in the piecewiseSEM package as sem.model.fits. See updates here:


Implementation of Schielzeth and Nakagawa's R2 for generalized linear mixed effects models in R. This function improves on the r.squaredGLMM function in the MuMIn package by incorporting different link functions for GLMERs and also returning other useful information, such as the model specification, and additional fit criteria in the form of AIC values.

For more information, see:

Nakagawa, Shinichi, and Holger Schielzeth. "A general and simple method for obtaining R2 from 
generalized linear mixed‐effects models." Methods in Ecology and Evolution 4.2 (2013): 133-142.

Johnson, Paul C.D. "Extension of Nakagawa & Schielzeth's R2GLMM to random slopes models." Methods in
Ecology and Evolution.

Version: 0.2-4 (2014-07-10)

Author: Jon Lefcheck & Juan Sebastian Casallas

Original blog post:

####WARNING: Calculation of nested random effects in lme objects still in beta!!


Generate a mock dataset

data <- data.frame(y=rnorm(100, 5, 10), y.binom=rbinom(100, 1, 0.5), y.poisson=rpois(100, 5))
data <- cbind(data, 
              data.frame(fixed1=data$y+c(runif(50, 0, 5),runif(50, 10, 50)),
                         fixed2=c("Treatment1", "Treatment2"),
                         rand2=rep(LETTERS[23:26],each=25)) )


#Linear model
mod0 <- lm(y ~ fixed1, data)
#Linear mixed effects model
mod1 <- lmer(y ~ fixed1 + (1|rand2/rand1), data)
mod1.1 <- lmer(y ~ fixed1 + (fixed1|rand2/rand1), data)
mod2 <- lmer(y ~ fixed1 + fixed2 + (1|rand2/rand1), data)
rsquared.glmm(list(mod0, mod1, mod1.1, mod2))
#Generalized linear mixed effects model (binomial)
mod3 <- glmer(y.binom ~ fixed1*fixed2 + (1|rand2/rand1), family="binomial", data)
mod3.prob <- update(mod3, family = binomial(link = "probit"))
rsquared.glmm(list(mod3, mod3.prob))
#Generalized linear mixed effects model (poisson)
mod4 <- glmer(y.poisson ~ fixed1*fixed2 + (1|rand2/rand1), family="poisson", data)
mod4.sqrt <- update(mod4, family = poisson(link = "sqrt"))
rsquared.glmm(list(mod4, mod4.sqrt))
#Get values for all kinds of models
(lme4.models <- rsquared.glmm(list(mod0, mod1, mod1.1, mod2, mod3, mod3.prob, mod4, mod4.sqrt)))

Compare output to MuMIn::r.squaredGLMM

MuMIn::r.squaredGLMM is similar to rsquared.glmm but returns less information and cannot handle different kinds of link functions:

(mumin.models <-, lapply(list(mod0, mod1, mod1.1, mod2, mod3, mod3.prob, mod4, 
mod4.sqrt), r.squaredGLMM)))


blme extends lme4, but yields different coefficients for random and fixed effects, which could explain the differences between their conditional r-squared values.

#Linear mixed effects model
blme.mod1 <- blmer(y ~ fixed1 + (1|rand2/rand1), data)
blme.mod1.1 <- blmer(y ~ fixed1 + (fixed1|rand2/rand1), data)
blme.mod2 <- blmer(y ~ fixed1 + fixed2 + (1|rand2/rand1), data)
#Generalized linear mixed effects model (binomial)
blme.mod3 <- bglmer(y.binom ~ fixed1*fixed2 + (1|rand2/rand1), family="binomial", data)
blme.mod3.prob <- update(blme.mod3, family = binomial(link = "probit"))
#Generalized linear mixed effects model (poisson)
blme.mod4 <- bglmer(y.poisson ~ fixed1*fixed2 + (1|rand2/rand1), family="poisson", data)
blme.mod4.sqrt <- update(blme.mod4, family = poisson(link = "sqrt"))
#Get values for all kinds of models
(blme.models <- rsquared.glmm(list(mod0, blme.mod1, blme.mod1.1, blme.mod2, blme.mod3, blme.mod3.prob,
# blme models yield better conditional r-squared values
all.equal(lme4.models[-(1:2)], blme.models[-(1:2)])


lmerTest::lmer extends lme4::lmer, but their random and mixed effects coefficients are the same.

# Try with lmerTest package -- output should be the same as above
#Linear mixed effects model
lmerTest.mod1 <- lmer(y ~ fixed1 + (1|rand2/rand1), data)
lmerTest.mod2 <- lmer(y ~ fixed1 + fixed2 + (1|rand2/rand1), data)
rsquared.glmm(list(mod0, lmerTest.mod1, lmerTest.mod2))
(lmerTest.models <- rsquared.glmm(list(mod0, lmerTest.mod1, lmerTest.mod2)))
# Same results
all.equal(lme4.models[1:3, -(1:3)], lmerTest.models[,-(1:3)])


nlme::lme and lme4::lmer yield very similar r-squared values.

lme.mod1 <- lme(y ~ fixed1, random=~1|rand2/rand1, data)
#Change to old optimizer to solve convergence issueslme4.models2[2:4, -(1:3)]
lme.mod1.1 <- lme(y ~ fixed1, random=~fixed1|rand2/rand1, control = lmeControl(opt = "optim"), data) 
lme.mod2 <- lme(y ~ fixed1 + fixed2, random=~1|rand2/rand1, data)
(lme.models <- rsquared.glmm(list(lme.mod1, lme.mod1.1, lme.mod2)))
# Compare to lme4 models, minor differences
all.equal(lme4.models2[2:4, -(1:3)], lme.models[,-(1:3)], tol = 1e-4)
You can’t perform that action at this time.