Make Cook's distance lines in the influential observations plot match check_outliers() (#964)

Co-authored-by: Daniel <mail@danielluedecke.de>

Author: jmgirard

NEWS.md

@@ -72,6 +72,17 @@
   Before, these functions returned inaccurate values or failed with an unclear
   error (#165).
 
+* The influential observations plot of `check_model()` and `check_outliers()`
+  for linear mixed models from package *lme4* now shows standardized residuals
+  instead of raw residuals. Before this fix, points fell outside the Cook's
+  distance lines but were not flagged as outliers (#857).
+
+* The Cook's distance lines in the influential observations plot of
+  `check_model()` and `check_outliers()` now use the threshold of
+  `check_outliers()`, including a user-supplied `threshold`. Before this fix,
+  linear models always used a threshold from the model's F statistic, so the
+  lines did not match the flagged outliers (#916).
+
 # performance 0.18.2
 
 ## Bug fixes

R/check_model_diagnostics.R

@@ -204,20 +204,35 @@
 .model_diagnostic_outlier <- function(model, threshold = NULL) {
   s <- summary(model)
 
-  if (inherits(model, "lm", which = TRUE) == 1) {
+  # use the threshold from `check_outliers()`, so that the contour lines
+  # match the flagged outliers
+  if (!is.null(threshold)) {
+    cook_levels <- threshold
+  } else if (inherits(model, "lm", which = TRUE) == 1) {
     cook_levels <- round(stats::qf(0.5, s$fstatistic[2], s$fstatistic[3]), 2)
-  } else if (is.null(threshold)) {
-    cook_levels <- c(0.5, 1)
   } else {
-    cook_levels <- threshold
+    cook_levels <- c(0.5, 1)
   }
 
   n_params <- tryCatch(model$rank, error = function(e) insight::n_parameters(model))
 
   infl <- stats::influence(model, do.coef = FALSE)
   model_resid <- as.numeric(insight::get_residuals(model))
 
-  std_resid <- tryCatch(stats::rstandard(model, infl), error = function(e) model_resid)
+  std_resid <- tryCatch(
+    stats::rstandard(model, infl),
+    error = function(e) {
+      # no rstandard() method (e.g., lmerMod): standardize Pearson residuals,
+      # so that the Cook's distance contours match `cooks.distance()`
+      .safe(
+        as.numeric(
+          stats::residuals(model, type = "pearson") /
+            (stats::sigma(model) * sqrt(1 - infl$hat))
+        ),
+        model_resid
+      )
+    }
+  )
 
   plot_data <- data.frame(
     Hat = infl$hat,

R/check_outliers.R

@@ -638,7 +638,7 @@ check_outliers.default <- function(
   attr(outlier, "text_size") <- 3
   attr(outlier, "influential_obs") <- .safe(.model_diagnostic_outlier(
     x,
-    threshold = unlist(thresholds)
+    threshold = thresholds$cook
   )) # nolint
   attr(outlier, "variables") <- "(Whole model)"
   attr(outlier, "raw_data") <- my_data

tests/testthat/test-check_outliers.R

@@ -479,6 +479,41 @@ test_that("check_outliers numeric, z-score", {
 })
 
 
+test_that("influential observations plot data uses standardized residuals for lmer models", {
+  skip_if_not_installed("lme4")
+
+  # the plot draws Cook's distance contours from leverage, standardized
+  # residuals and the number of parameters, so these must reproduce the
+  # Cook's distance that check_outliers() tests (#857)
+  m <- lme4::lmer(Reaction ~ Days + (Days | Subject), data = lme4::sleepstudy)
+  out <- check_outliers(m)
+  d <- attributes(out)$influential_obs
+  p <- attributes(d)$n_params
+  implied_cook <- d$Std_Residuals^2 * d$Hat / (p * (1 - d$Hat))
+  expect_equal(implied_cook, d$Cooks_Distance, tolerance = 1e-6, ignore_attr = TRUE)
+  expect_identical(
+    which(implied_cook > attributes(out)$threshold$cook),
+    which(as.numeric(out) == 1)
+  )
+})
+
+
+test_that("influential observations plot data uses the Cook's distance threshold for lm models", {
+  # the contour lines must follow the threshold that check_outliers() uses,
+  # including a user-supplied threshold (#916)
+  m <- lm(disp ~ mpg + hp, data = mtcars)
+  for (th in list(NULL, list(cook = 2))) {
+    out <- check_outliers(m, method = "cook", threshold = th)
+    d <- attributes(out)$influential_obs
+    expect_equal(attributes(d)$cook_levels, attributes(out)$threshold$cook)
+    expect_identical(
+      which(d$Influential == "Influential"),
+      which(as.numeric(out) == 1)
+    )
+  }
+})
+
+
 test_that("check_outliers with psych", {
   skip_if_not_installed("psych")
   data(Harman.5, package = "psych")