3
$\begingroup$

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: 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:

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)

mod_ref <- glmmTMB(y~x+(1|f),family=list(family="beta",link="logit"),
          data=dd)
mod_0 <- glmmTMB(y~1+(1|f),family=list(family="beta",link="logit"),
          data=dd)
$\endgroup$
11
  • $\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$ 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$ 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$ 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$ 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
1
$\begingroup$

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.

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.