Add autocorrelation testing for simulated residuals and use simulated residuals for Poisson mixed models

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: performance
 Title: Assessment of Regression Models Performance
-Version: 0.15.2
+Version: 0.15.2.1
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -159,7 +159,7 @@ Suggests:
     withr (>= 3.0.0)
 Encoding: UTF-8
 Language: en-US
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Roxygen: list(markdown = TRUE)
 Config/testthat/edition: 3
 Config/testthat/parallel: true

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

@@ -1,3 +1,8 @@
+# performance (devel)
+
+* `check_autocorrelation()` gets methods for `DHARMa` objects and objects from
+  `simulate_residuals()`.
+
 # performance 0.15.2
 
 ## Bug fixes

R/binned_residuals.R

@@ -75,16 +75,18 @@
 #' }
 #'
 #' @export
-binned_residuals <- function(model,
-                             term = NULL,
-                             n_bins = NULL,
-                             show_dots = NULL,
-                             ci = 0.95,
-                             ci_type = "exact",
-                             residuals = "deviance",
-                             iterations = 1000,
-                             verbose = TRUE,
-                             ...) {
+binned_residuals <- function(
+  model,
+  term = NULL,
+  n_bins = NULL,
+  show_dots = NULL,
+  ci = 0.95,
+  ci_type = "exact",
+  residuals = "deviance",
+  iterations = 1000,
+  verbose = TRUE,
+  ...
+) {
   ci_type <- insight::validate_argument(
     ci_type,
     c("exact", "gaussian", "boot")
@@ -98,7 +100,9 @@ binned_residuals <- function(model,
   if (isFALSE(insight::model_info(model)$is_bernoulli)) {
     ci_type <- "gaussian"
     if (verbose) {
-      insight::format_alert("Using `ci_type = \"gaussian\"` because model is not bernoulli.")
+      insight::format_alert(
+        "Using `ci_type = \"gaussian\"` because model is not bernoulli."
+      )
     }
   }
 
@@ -121,18 +125,23 @@ binned_residuals <- function(model,
   y0 <- .recode_to_zero(insight::get_response(model, verbose = FALSE))
 
   # calculate residuals
-  y <- switch(residuals,
+  y <- switch(
+    residuals,
     response = y0 - fitted_values,
     pearson = .safe((y0 - fitted_values) / sqrt(fitted_values * (1 - fitted_values))),
     deviance = .safe(stats::residuals(model, type = "deviance"))
   )
 
   # make sure we really have residuals
   if (is.null(y)) {
-    insight::format_error("Could not calculate residuals. Try using `residuals = \"response\"`.")
+    insight::format_error(
+      "Could not calculate residuals. Try using `residuals = \"response\"`."
+    )
   }
 
-  if (is.null(n_bins)) n_bins <- round(sqrt(length(pred)))
+  if (is.null(n_bins)) {
+    n_bins <- round(sqrt(length(pred)))
+  }
 
   breaks.index <- floor(length(pred) * (1:(n_bins - 1)) / n_bins)
   breaks <- unique(c(-Inf, sort(pred)[breaks.index], Inf))
@@ -151,8 +160,13 @@ binned_residuals <- function(model,
     if (n == 0) {
       conf_int <- stats::setNames(c(NA, NA), c("CI_low", "CI_high"))
     } else {
-      conf_int <- switch(ci_type,
-        gaussian = stats::qnorm(c((1 - ci) / 2, (1 + ci) / 2), mean = ybar, sd = sdev / sqrt(n)),
+      conf_int <- switch(
+        ci_type,
+        gaussian = stats::qnorm(
+          c((1 - ci) / 2, (1 + ci) / 2),
+          mean = ybar,
+          sd = sdev / sqrt(n)
+        ),
         exact = {
           out <- stats::binom.test(sum(y0[items]), n)$conf.int
           # center CIs around point estimate

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.
@@ -18,6 +21,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,11 +57,41 @@ check_autocorrelation.default <- function(x, nsim = 1000, ...) {
 }
 
 
+#' @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 ------------------------------
 
 #' @export
 plot.check_autocorrelation <- function(x, ...) {
-  insight::format_warning("There is currently no `plot()` method for `check_autocorrelation()`.")
+  insight::format_warning(
+    "There is currently no `plot()` method for `check_autocorrelation()`."
+  )
 }
 
 

R/check_clusterstructure.R

@@ -37,10 +37,7 @@
 #'   assessment of (cluster) tendency. In Proceedings of the 2002 International
 #'   Joint Conference on Neural Networks. IJCNN02 (3), 2225-2230. IEEE.
 #' @export
-check_clusterstructure <- function(x,
-                                   standardize = TRUE,
-                                   distance = "euclidean",
-                                   ...) {
+check_clusterstructure <- function(x, standardize = TRUE, distance = "euclidean", ...) {
   if (standardize) {
     x <- as.data.frame(scale(x))
   }
@@ -64,13 +61,22 @@ check_clusterstructure <- function(x,
 
   out <- list(
     H = H,
-    dissimilarity_matrix = .clusterstructure_dm(x, distance = distance, method = "ward.D2")
+    dissimilarity_matrix = .clusterstructure_dm(
+      x,
+      distance = distance,
+      method = "ward.D2"
+    )
   )
 
   attr(out, "text") <- res_text
   attr(out, "color") <- color
   attr(out, "title") <- "Clustering tendency"
-  class(out) <- c("see_check_clusterstructure", "check_clusterstructure", "easystats_check", class(out))
+  class(out) <- c(
+    "see_check_clusterstructure",
+    "check_clusterstructure",
+    "easystats_check",
+    class(out)
+  )
   out
 }
 
@@ -80,8 +86,10 @@ plot.check_clusterstructure <- function(x, ...) {
   # Can be reimplemented with ggplot in see
   stats::heatmap(
     x$dissimilarity_matrix,
-    Rowv = NA, Colv = NA,
-    labRow = FALSE, labCol = FALSE,
+    Rowv = NA,
+    Colv = NA,
+    labRow = FALSE,
+    labCol = FALSE,
     col = grDevices::colorRampPalette(c("#2196F3", "#FAFAFA", "#E91E63"))(100)
   )
 }