I'm looking for a method or function for computing R² for glmmTMB models with a beta distribution and a logit link. I am interested in a ratio (%) response in a repeated measures design.

I looked into:

  • the sem.model.fits() function from the package {piecewiseSEM}
  • the r.squaredGLMM() function from the package {MuMIn}
  • the procedure described by Shinichi Nakagawa, Paul C. D. Johnson & Holger Schielzeth in “The coefficient of determination R²2 and intra-class correlation ICC from generalized linear-mixed effects models revisited and expanded” published in September 2017 (DOI: 10.1098/rsif.2017.0213) or by Paul C. D. Johnson in "Extension of Nakagawa & Schielzeth's R2GLMM to random slopes models" published in July 2014 (https://doi.org/10.1111/2041-210X.12225).

The procedure described by Shinichi Nakagawa (for glmer models) is the following:

mod_ref=FULL MODEL (containing all random and fixed terms)
mod_0= NULL MODEL (containing only random terms)
VarF <- var(as.vector(model.matrix(mod_ref) %*% fixef(mod_ref)))
nuN <- 1/attr(VarCorr(mod_0), "sc")^2 # note that glmer report 1/nu not nu as resiudal variance
VarOdN <- 1/nuN # the delta method
VarOlN <- log(1 + 1/nuN) # log-normal approximation
VarOtN <- trigamma(nuN) # trigamma function
c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
nuF <- 1/attr(VarCorr(mod_ref), "sc")^2 # note that glmer report 1/nu not nu as resiudal variance
VarOdF <- 1/nuF # the delta method
VarOlF <- log(1 + 1/nuF) # log-normal approximation
VarOtF <- trigamma(nuF) # trigamma function
c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
R2glmmM <- VarF/(VarF + sum(as.numeric(VarCorr(mod_ref))) + VarOtF)
R2glmmC <- (VarF + sum(as.numeric(VarCorr(mod_ref))))/(VarF +sum(as.numeric(VarCorr(mod_ref)))+VarOtF)

The procedure described by Paul Johnson is the following (mod=full model):

X <- model.matrix(mod_ref)
n <- nrow(X)
Beta <- fixef(mod_ref)
Sf <- var(X %*% Beta)
Sigma.list <- VarCorr(mod_ref)
Sl <- 
      Z <-X[,rownames(Sigma)]
      sum(diag(Z %*% Sigma %*% t(Z)))/n
  Se <- attr(Sigma.list, "sc")^2
  Sd <- 0
  total.var <- Sf + Sl + Se + Sd
  (Rsq.m <- Sf / total.var) 
  (Rsq.c <- (Sf + Sl) / total.var) 

However, I was unable to find this particular combination in any of the above packages or article.

When I try to run these two scripts on my glmmTMB model : regd=glmmTMB(Ratio~block+treatment*scale(year)+(1|year)+(1|plot)+(1|year:plot),family=list(family="beta",link="logit"),data=biomass), I get the following error message: Error in model.matrix(mod_ref) %*% fixef(mod_ref) : requires numeric/complex matrix/vector arguments at line 3 of Nakagawa's method and Error in X %*% Beta : requires numeric/complex matrix/vector arguments at line 4 of Johnson's method.

When I check model.matrix(mod_ref), it yields the following matrix: enter image description here

When I check fixef(mod_ref), I get the following: enter image description here

Could this be easily adapted?

A short reproducible example:

ngrp <- 100
eps <- 0.001
nobs <- 4000
ngrp <- 200
nobs <- 800
x <- rnorm(nobs)
f <- factor(rep(1:ngrp,nobs/ngrp))
u <- rnorm(ngrp,sd=1)
eta <- 2+x+u[f]
mu <- plogis(eta)
y <- rbeta(nobs,shape1=mu/0.1,shape2=(1-mu)/0.1)
y <- pmin(1-eps,pmax(eps,y))
dd <- data.frame(x,y,f)

mod_ref <- glmmTMB(y~x+(1|f),family=list(family="beta",link="logit"),
mod_0 <- glmmTMB(y~1+(1|f),family=list(family="beta",link="logit"),
  • $\begingroup$ Is this a mixed model? Or a single level model? If you use the framework of Nakagawa, Schielzeth and Johnson, then it is possible. $\endgroup$ – Heteroskedastic Jim Jul 27 '18 at 11:22
  • $\begingroup$ @user162986 I edited the post in a more detailed way because I tried to use the framework for glmmTMB and ended up with some error messages. Maybe you have some idea as to how to adapt this? Thanks for your interest $\endgroup$ – Guillaume A2 Aug 3 '18 at 9:02
  • 1
    $\begingroup$ That error usually occurs when one of the objects you're multiplying is not a matrix or contains text. Check the result of model.matrix(mod_ref) to be sure it is what you think it is. $\endgroup$ – Heteroskedastic Jim Aug 3 '18 at 11:39
  • $\begingroup$ @user162986 Thanks again for pushing things a little further. I added the output of model.matrix(mod_ref) and fixef(mod_ref). They seem pretty typical to me but I rarely dig this deep! $\endgroup$ – Guillaume A2 Aug 3 '18 at 14:01
  • 1
    $\begingroup$ I'm not sure, doesn't glmmTMB return a list for VarCorr(), because it always returns an element for the conditional and the possible zero-inflated model. So you may need VarCorr()[[1]] here, but a reprex would make debugging-life easier. :-) $\endgroup$ – Daniel Aug 12 '18 at 10:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.