update

R/compute_variances.R

@@ -767,7 +767,7 @@
 
   # transform expected mean of the null model
   if (is.null(faminfo$family)) {
-    mu <- exp(mu)
+    mu <- link_inverse(x)(mu)
   } else {
     # transform mu
     mu <- switch(faminfo$family,
@@ -868,7 +868,7 @@
 
   switch(approx_method,
     delta = cvsquared,
-    trimamma = trigamma(cvsquared),
+    trigamma = trigamma(cvsquared),
     log1p(cvsquared)
   )
 }

WIP/nakagawa_suppl.R

@@ -384,3 +384,48 @@ c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
 
 performance::r2_nakagawa(parmodGLMERf)
 MuMIn::r.squaredGLMM(parmodGLMERf)
+
+
+
+
+
+
+# Fit null model without fixed effects (but including all random effects)
+glmmTMBr <- glmmTMB::glmmTMB(
+  count ~ (1 | site),
+  family = glmmTMB::nbinom1(),
+  data = Salamanders, REML = TRUE
+)
+glmmTMBf <- glmmTMB::glmmTMB(
+  count ~ mined + spp + (1 | site),
+  family = glmmTMB::nbinom1(),
+  data = Salamanders, REML = TRUE
+)
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% fixef(glmmTMBf)$cond))
+# getting the observation-level variance Null model
+thetaN <- sigma(glmmTMBr)
+lambda <- as.numeric(exp(fixef(glmmTMBr)$cond + 0.5 * (as.numeric(VarCorr(glmmTMBr)$cond[1]))))
+# lambda2 <- mean(DataAll$Parasite)
+VarOdN <- 1 / lambda + 1 / thetaN # the delta method
+VarOlN <- log(1 + (1 / lambda) + (1 / thetaN)) # log-normal approximation
+VarOtN <- trigamma((1 / lambda + 1 / thetaN)^(-1)) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
+
+performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+MuMIn::r.squaredGLMM(glmmTMBf)
+get_variance(glmmTMBf, null_model = glmmTMBr)