fix

easystats · Jun 13, 2024 · 6e687b3 · 6e687b3
1 parent b2dcc63
commit 6e687b3
Show file tree

Hide file tree

Showing 6 changed files with 129 additions and 18 deletions.
diff --git a/R/compute_variances.R b/R/compute_variances.R
@@ -4,7 +4,8 @@
                                name_full = NULL,
                                verbose = TRUE,
                                tolerance = 1e-8,
-                               model_component = "conditional") {
+                               model_component = "conditional",
+                               model_null = NULL) {
   ## Original code taken from GitGub-Repo of package glmmTMB
   ## Author: Ben Bolker, who used an cleaned-up/adapted
   ## version of Jon Lefcheck's code from SEMfit
@@ -44,7 +45,9 @@
 
   # we also need necessary model information, like fixed and random effects,
   # variance-covariance matrix etc. for the null model
-  model_null <- .safe(null_model(x, verbose = FALSE))
+  if (is.null(model_null)) {
+    model_null <- .safe(null_model(x, verbose = FALSE))
+  }
   vals_null <- .get_variance_information(
     model_null,
     faminfo = faminfo,
@@ -55,7 +58,9 @@
 
   # Test for non-zero random effects ((near) singularity)
   no_random_variance <- FALSE
-  if (performance::check_singularity(x, tolerance = tolerance) && !(component %in% c("slope", "intercept"))) {
+  singular_fit <- isTRUE(.safe(performance::check_singularity(x, tolerance = tolerance)))
+
+  if (singular_fit && !(component %in% c("slope", "intercept"))) {
     if (verbose) {
       format_warning(
         sprintf("Can't compute %s. Some variance components equal zero. Your model may suffer from singularity (see `?lme4::isSingular` and `?performance::check_singularity`).", name_full), # nolint
@@ -101,7 +106,7 @@
   }
 
   # Variance of random effects for NULL model
-  if (!performance::check_singularity(model_null, tolerance = tolerance) && !is.null(vals_null)) {
+  if (!singular_fit && !is.null(vals_null)) {
     # Separate observation variance from variance of random effects
     nr <- vapply(vals_null$re, nrow, numeric(1))
     not.obs.terms_null <- names(nr[nr != n_obs(model_null)])
@@ -294,7 +299,7 @@
     )
     names(vals$re) <- re_names[seq_along(vals$re)]
 
-    # nlme
+    # nlme / glmmPQL
     # ---------------------------
   } else if (inherits(x, "lme")) {
     re_names <- find_random(x, split_nested = TRUE, flatten = TRUE)
@@ -697,6 +702,7 @@
       beta = ,
       ordbeta = stats::plogis(mu),
       poisson = ,
+      quasipoisson = ,
       nbinom = ,
       nbinom1 = ,
       nbinom2 = ,
@@ -730,6 +736,7 @@
         nbinom = ,
         nbinom1 = ,
         nbinom2 = ,
+        quasipoisson = ,
         `negative binomial` = sig,
         `zero-inflated negative binomial` = ,
         genpois = .variance_family_nbinom(x, mu, sig, faminfo),
@@ -760,6 +767,7 @@
         `negative binomial` = (1 / mu) + (1 / sig),
         nbinom = ,
         poisson = ,
+        quasipoisson = ,
         nbinom1 = vv / mu,
         vv / mu^2
       )

diff --git a/R/get_nested_lme_varcorr.R b/R/get_nested_lme_varcorr.R
@@ -49,5 +49,9 @@
 
 
 .is_nested_lme <- function(x) {
-  sapply(find_random(x), function(i) any(grepl(":", i, fixed = TRUE)))
+  if (inherits(x, "glmmPQL")) {
+    length(find_random(x, flatten = TRUE)) > 1
+  } else {
+    sapply(find_random(x), function(i) any(grepl(":", i, fixed = TRUE)))
+  }
 }
diff --git a/R/get_sigma.R b/R/get_sigma.R
@@ -198,7 +198,7 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
 
 
 .get_sigma.glmmPQL <- function(x, ...) {
-  switch(stats::family(x)$family,
+  switch(x$family$family,
     gaussian = ,
     Gamma = x$sigma,
     x$sigma^2

diff --git a/R/get_variances.R b/R/get_variances.R
@@ -10,15 +10,17 @@
 #'
 #' @param x A mixed effects model.
 #' @param component Character value, indicating the variance component that should
-#'   be returned. By default, all variance components are returned. The
-#'   distribution-specific (`"distribution"`) and residual (`"residual"`)
-#'   variance are the most computational intensive components, and hence may
-#'   take a few seconds to calculate.
+#' be returned. By default, all variance components are returned. The
+#' distribution-specific (`"distribution"`) and residual (`"residual"`)
+#' variance are the most computational intensive components, and hence may
+#' take a few seconds to calculate.
 #' @param verbose Toggle off warnings.
 #' @param tolerance Tolerance for singularity check of random effects, to decide
-#'   whether to compute random effect variances or not. Indicates up to which
-#'   value the convergence result is accepted. The larger tolerance is, the
-#'   stricter the test will be. See [performance::check_singularity()].
+#' whether to compute random effect variances or not. Indicates up to which
+#' value the convergence result is accepted. The larger tolerance is, the
+#' stricter the test will be. See [performance::check_singularity()].
+#' @param null_model Optional, a null-model to be used for the calculation of
+#' random effect variances. If `NULL`, the null-model is computed internally.
 #' @param ... Currently not used.
 #'
 #' @return A list with following elements:
@@ -135,7 +137,6 @@ get_variance <- function(x,
   UseMethod("get_variance")
 }
 
-
 #' @export
 get_variance.default <- function(x,
                                  component = c(
@@ -152,6 +153,7 @@ get_variance.default <- function(x,
 }
 
 
+#' @rdname get_variance
 #' @export
 get_variance.merMod <- function(x,
                                 component = c(
@@ -161,6 +163,7 @@ get_variance.merMod <- function(x,
                                 ),
                                 verbose = TRUE,
                                 tolerance = 1e-5,
+                                null_model = NULL,
                                 ...) {
   component <- match.arg(component)
   .safe(.compute_variances(
@@ -169,7 +172,8 @@ get_variance.merMod <- function(x,
     name_fun = "get_variance",
     name_full = "random effect variances",
     verbose = verbose,
-    tolerance = tolerance
+    tolerance = tolerance,
+    model_null = null_model
   ))
 }
 
@@ -205,6 +209,7 @@ get_variance.lme <- get_variance.merMod
 get_variance.brmsfit <- get_variance.merMod
 
 
+#' @rdname get_variance
 #' @export
 get_variance.glmmTMB <- function(x,
                                  component = c(
@@ -214,6 +219,7 @@ get_variance.glmmTMB <- function(x,
                                  ),
                                  verbose = TRUE,
                                  tolerance = 1e-5,
+                                 null_model = NULL,
                                  model_component = NULL,
                                  ...) {
   component <- match.arg(component)
@@ -224,7 +230,8 @@ get_variance.glmmTMB <- function(x,
     name_full = "random effect variances",
     verbose = verbose,
     tolerance = tolerance,
-    model_component = model_component
+    model_component = model_component,
+    model_null = null_model
   ))
 }
 
@@ -240,6 +247,7 @@ get_variance.mixed <- function(x,
                                  "slope", "rho01", "rho00"
                                ),
                                verbose = TRUE,
+                               null_model = NULL,
                                tolerance = 1e-5,
                                ...) {
   component <- match.arg(component)
@@ -249,7 +257,8 @@ get_variance.mixed <- function(x,
     name_fun = "get_variance",
     name_full = "random effect variances",
     verbose = verbose,
-    tolerance = tolerance
+    tolerance = tolerance,
+    model_null = null_model
   ))
 }
 

diff --git a/WIP/nakagawa_suppl.R b/WIP/nakagawa_suppl.R
@@ -1,3 +1,5 @@
+vsCodeSnippets::load_debug_pkg()
+
 # installing glmmADMB
 library(R2admb)
 library(glmmADMB)
@@ -93,7 +95,11 @@ ColourL <- with(DataMale, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 *
 DataMale$Colour <- rbinom(length(ColourL), 1, plogis(ColourL))
 
 
+
+# ==============================================================
 # Quasi-Poisson with log link (page 5) -------------------------
+# glmmadmb -----------------------------------------------------
+# ==============================================================
 
 # Fit null model without fixed effects (but including all random effects)
 fecmodADMBr <- glmmadmb(Egg ~ 1 + (1 | Population) + (1 | Container), family = "nbinom1",
@@ -139,3 +145,61 @@ ICCadjCont <- VarCorr(fecmodADMBf)$Container[1]/(sum(as.numeric(VarCorr(fecmodAD
 c(R2glmmM = R2glmmM, R2glmmC = R2glmmC, ICCrawPop = ICCrawPop, ICCadjPop = ICCadjPop, ICCrawCont = ICCrawCont, ICCadjCont = ICCadjCont)
 
 performance::r2_nakagawa(fecmodADMBf)
+
+
+# ==============================================================
+# Quasi-Poisson with log link (page 5) -------------------------
+# glmmPQL- -----------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+fecmodPQLr <- glmmPQL(Egg ~ 1,
+  random = list(~ 1 | Population, ~ 1 | Container),
+  family = "quasipoisson", data = DataFemale
+)
+# Fit alternative model including fixed and all random effects
+fecmodPQLf <- glmmPQL(Egg ~ Treatment + Habitat, random = list(
+  ~ 1 | Population,
+  ~ 1 | Container
+), family = "quasipoisson", data = DataFemale)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(~ Treatment + Habitat, data = DataFemale) %*% fixef(fecmodPQLf)))
+# getting the observation-level variance Null model
+omegaN <- as.numeric(VarCorr(fecmodPQLr)[5, 1]) # overdispersion omega is residual variance in glmmPQL
+lambda <- as.numeric(exp(fixef(fecmodPQLr) + 0.5 * (as.numeric(VarCorr(fecmodPQLr)[2, 1]) + as.numeric(VarCorr(fecmodPQLr)[4, 1]))))
+# lambda2 <- mean(DataFemale$Egg)
+VarOdN <- omegaN / lambda
+
+VarOlN <- log(1 + omegaN / lambda)
+VarOtN <- trigamma(lambda / omegaN)
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+omegaF <- as.numeric(VarCorr(fecmodPQLf)[5, 1]) # overdispersion omega is residual variance in glmmPQL
+VarOdF <- omegaF / lambda
+VarOlF <- log(1 + omegaF / lambda)
+VarOtF <- trigamma(lambda / omegaF)
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1]))) / (VarF + sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# Raw unadjusted ICC[Population]
+ICCrawPop <- as.numeric(VarCorr(fecmodPQLr)[2, 1]) / (sum(as.numeric(VarCorr(fecmodPQLr)[c(2, 4), 1])) + VarOlN)
+# adjusted ICC[Population]
+ICCadjPop <- as.numeric(VarCorr(fecmodPQLf)[2, 1]) / (sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# Raw unadjusted ICC[Container]
+ICCrawCont <- as.numeric(VarCorr(fecmodPQLr)[4, 1]) / (sum(as.numeric(VarCorr(fecmodPQLr)[c(2, 4), 1])) + VarOlN)
+# adjusted ICC[Container]
+ICCadjCont <- as.numeric(VarCorr(fecmodPQLf)[4, 1]) / (sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# comparing the results
+c(
+  R2glmmM = R2glmmM, R2glmmC = R2glmmC, ICCrawPop = ICCrawPop, ICCadjPop = ICCadjPop,
+  ICCrawCont = ICCrawCont, ICCadjCont = ICCadjCont
+)
+
+performance::r2_nakagawa(fecmodPQLf, null_model = fecmodPQLr)
+get_sigma(fecmodPQLf)
diff --git a/man/get_variance.Rd b/man/get_variance.Rd