check_autocorrelation() gets methods for DHARMa objects and objects from simulate_residuals(). (#864)

Author: strengejacke

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: performance
 Title: Assessment of Regression Models Performance
-Version: 0.15.2.1
+Version: 0.15.2.2
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",

NAMESPACE

@@ -28,7 +28,9 @@ S3method(as.data.frame,r2_nakagawa)
 S3method(as.double,check_outliers)
 S3method(as.double,item_omega)
 S3method(as.double,performance_roc)
+S3method(check_autocorrelation,DHARMa)
 S3method(check_autocorrelation,default)
+S3method(check_autocorrelation,performance_simres)
 S3method(check_collinearity,BFBayesFactor)
 S3method(check_collinearity,MixMod)
 S3method(check_collinearity,afex_aov)

NEWS.md

@@ -2,6 +2,9 @@
 
 ## Changes
 
+* `check_autocorrelation()` gets methods for `DHARMa` objects and objects from
+  `simulate_residuals()`.
+
 * Improved documentation for printing-methods.
 
 # performance 0.15.2

R/check_autocorrelation.R

@@ -4,9 +4,12 @@
 #' @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 time A vector with time values to specify the temporal order of the data.
+#' Only used if `x` is an object returned by `simulate_residuals()` or by `DHARMa`.
+#' @param ... Currently not used for models. For simulated residuals, arguments are
+#' passed to `DHARMa::testTemporalAutocorrelation()`.
 #'
 #' @return Invisibly returns the p-value of the test statistics. A p-value < 0.05
 #' indicates autocorrelated residuals.
@@ -48,6 +51,32 @@ check_autocorrelation.default <- function(x, nsim = 1000, ...) {
   p.val
 }
 
+#' @rdname check_autocorrelation
+#' @export
+check_autocorrelation.performance_simres <- function(x, time = NULL, ...) {
+  insight::check_if_installed("DHARMa")
+
+  if (is.null(time)) {
+    insight::format_warning(
+      "Data are assumed to be ordered by time. If this is not the case, please provide a `time` argument."
+    )
+    time <- seq_along(x$scaledResiduals)
+  }
+
+  # 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, time = time, 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 ------------------------------
 

man/check_autocorrelation.Rd

@@ -0,0 +1,53 @@
+test_that("check_autocorrelation works with simulated residuals", {
+  skip_if_not_installed("DHARMa")
+  skip_if_not_installed("glmmTMB")
+  skip_if_not(getRversion() >= "4.0.0")
+
+  data(Salamanders, package = "glmmTMB")
+
+  # Test with a simple Poisson GLM
+  m <- glm(count ~ spp + mined, family = poisson, data = Salamanders)
+
+  # Simulate residuals
+  set.seed(123)
+  simres <- simulate_residuals(m)
+
+  # Check autocorrelation
+  set.seed(123)
+  expect_warning(
+    {
+      out <- check_autocorrelation(simres)
+    },
+    regex = "Data are assumed"
+  )
+
+  # Should return a p-value
+  expect_type(out, "double")
+  expect_s3_class(out, "check_autocorrelation")
+
+  expect_equal(as.vector(out), 0.2211415, tolerance = 1e-3, ignore_attr = TRUE)
+})
+
+
+test_that("check_autocorrelation.DHARMa works", {
+  skip_if_not_installed("DHARMa")
+
+  # Test that the DHARMa method works
+  data(mtcars)
+  m <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars)
+
+  set.seed(123)
+  simres <- DHARMa::simulateResiduals(m, plot = FALSE)
+
+  expect_warning(check_autocorrelation(simres), regex = "Data are assumed")
+  set.seed(123)
+  expect_silent({
+    out <- check_autocorrelation(simres, time = seq_along(simres$scaledResiduals))
+  })
+
+  # Should return a p-value
+  expect_type(out, "double")
+  expect_s3_class(out, "check_autocorrelation")
+
+  expect_equal(as.vector(out), 0.4163168, tolerance = 1e-3, ignore_attr = TRUE)
+})