# R-squared for glmmTMB with beta distribution and logit link

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 <-
sum(
sapply(Sigma.list,
function(Sigma)
{
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: When I check fixef(mod_ref), I get the following: A short reproducible example:

library(glmmTMB)
set.seed(101)
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)


• 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. Aug 3 '18 at 11:39
• 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()[] here, but a reprex would make debugging-life easier. :-) Aug 12 '18 at 10:29
Both the results from fixef and VarCorr are lists with $cond for conditioning variables, $zi zero inflation, and $disp for the Beta dispersion (parametrized as mean + dispersion). Using fixef(mod_ref)$cond and Sigma.list\$cond gives you no error.