Implement check_autocorrelation for simulated residuals and use simulated residuals for Poisson mixed models

Co-authored-by: strengejacke <26301769+strengejacke@users.noreply.github.com>

Author: Copilot

R/check_autocorrelation.R

@@ -4,9 +4,11 @@
 #' @description Check model for independence of residuals, i.e. for autocorrelation
 #' of error terms.
 #'
-#' @param x A model object.
+#' @param x A model object, or an object returned by `simulate_residuals()`.
 #' @param nsim Number of simulations for the Durbin-Watson-Test.
-#' @param ... Currently not used.
+#' @param ... Currently not used for models. For simulated residuals, arguments are
+#' passed to `DHARMa::testTemporalAutocorrelation()`, which can include `time` (a
+#' vector with time values) to specify the temporal order of the data.
 #'
 #' @return Invisibly returns the p-value of the test statistics. A p-value < 0.05
 #' indicates autocorrelated residuals.
@@ -18,6 +20,11 @@
 #' results for the estimates, or maybe a mixed model with error term for the
 #' cluster groups should be used.
 #'
+#' For simulated residuals (from `simulate_residuals()`), the function uses
+#' `DHARMa::testTemporalAutocorrelation()` to check for temporal autocorrelation.
+#' This requires the data to be ordered by time. If the data are not ordered by
+#' time, you can provide a `time` argument to specify the temporal order.
+#'
 #' @examples
 #' m <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
 #' check_autocorrelation(m)
@@ -49,6 +56,27 @@ check_autocorrelation.default <- function(x, nsim = 1000, ...) {
 }
 
 
+#' @rdname check_autocorrelation
+#' @export
+check_autocorrelation.performance_simres <- function(x, ...) {
+  insight::check_if_installed("DHARMa")
+  
+  # Use DHARMa's temporal autocorrelation test
+  # This requires the residuals to be ordered by time
+  # DHARMa::testTemporalAutocorrelation expects a DHARMa object
+  result <- DHARMa::testTemporalAutocorrelation(x, plot = FALSE, ...)
+  
+  # Extract p-value from the result
+  p.val <- result$p.value
+  
+  class(p.val) <- c("check_autocorrelation", "see_check_autocorrelation", class(p.val))
+  p.val
+}
+
+#' @export
+check_autocorrelation.DHARMa <- check_autocorrelation.performance_simres
+
+
 # methods ------------------------------
 
 #' @export

R/check_overdispersion.R

@@ -40,10 +40,10 @@
 #' [GLMM FAQ](http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html),
 #' section *How can I deal with overdispersion in GLMMs?*. Note that this
 #' function only returns an *approximate* estimate of an overdispersion
-#' parameter. Using this approach would be inaccurate for zero-inflated or
-#' negative binomial mixed models (fitted with `glmmTMB`), thus, in such cases,
-#' the overdispersion test is based on [`simulate_residuals()`] (which is identical
-#' to `check_overdispersion(simulate_residuals(model))`).
+#' parameter. For Poisson, zero-inflated, or negative binomial mixed models
+#' (fitted with `glmmTMB` or `lme4`), the overdispersion test is based on
+#' [`simulate_residuals()`] (which is identical to
+#' `check_overdispersion(simulate_residuals(model))`).
 #'
 #' @inheritSection check_zeroinflation Tests based on simulated residuals
 #'
@@ -242,7 +242,8 @@ check_overdispersion.merMod <- function(x, ...) {
   obj_name <- insight::safe_deparse_symbol(substitute(x))
 
   # for certain distributions, simulated residuals are more accurate
-  use_simulated <- info$family == "genpois" || info$is_zero_inflated || info$is_bernoulli || info$is_binomial || (!info$is_count && !info$is_binomial) || info$is_negbin # nolint
+  # Note: now including Poisson models for mixed models (see #595, #643)
+  use_simulated <- info$family == "genpois" || info$is_zero_inflated || info$is_bernoulli || info$is_binomial || (!info$is_count && !info$is_binomial) || info$is_negbin || info$is_poisson # nolint
 
   if (use_simulated) {
     return(check_overdispersion(simulate_residuals(x, ...), object_name = obj_name, ...))

R/check_zeroinflation.R

@@ -25,9 +25,9 @@
 #' negative binomial or zero-inflated models.
 #'
 #' In case of negative binomial models, models with zero-inflation component,
-#' or hurdle models, the results from `check_zeroinflation()` are based on
-#' [`simulate_residuals()`], i.e. `check_zeroinflation(simulate_residuals(model))`
-#' is internally called if necessary.
+#' hurdle models, or Poisson mixed models, the results from `check_zeroinflation()`
+#' are based on [`simulate_residuals()`], i.e.
+#' `check_zeroinflation(simulate_residuals(model))` is internally called if necessary.
 #'
 #' @section Tests based on simulated residuals:
 #' For certain models, resp. model from certain families, tests are based on
@@ -88,9 +88,12 @@ check_zeroinflation.default <- function(x, tolerance = 0.05, ...) {
   # model classes not supported in DHARMa
   not_supported <- c("fixest", "glmx")
 
-  # for models with zero-inflation component or negative binomial families,
-  # we use simulate_residuals()
-  if (!inherits(x, not_supported) && (model_info$is_zero_inflated || model_info$is_negbin || model_info$family == "genpois")) { # nolint
+  # for models with zero-inflation component, negative binomial families,
+  # or Poisson mixed models, we use simulate_residuals()
+  # Note: now including Poisson mixed models (see #595, #643)
+  use_simulated <- model_info$is_zero_inflated || model_info$is_negbin || model_info$family == "genpois" || (model_info$is_mixed && model_info$is_poisson) # nolint
+  
+  if (!inherits(x, not_supported) && use_simulated) {
     if (missing(tolerance)) {
       tolerance <- 0.1
     }