DHARMa implementation for new `check_residuals()` function

[strengejacke]

Fixes #595

I'm not sure, but I think this PR could also be related to following issues: #688, #654, #501, #500, #464, #632, #613, #408, #405, #376, #367, #274, #263

[strengejacke]

@mccarthy-m-g I invited you to get write access, which should allow you to push to this PR, so we can collaboratively work on this.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 58.77193% with 94 lines in your changes are missing coverage. Please review.

Project coverage is 59.31%. Comparing base (afddb29) to head (4127d9d).

Files	Patch %	Lines
R/check_model_diagnostics.R	0.00%	37 Missing :warning:
R/check_model.R	40.00%	18 Missing :warning:
R/simulate_residuals.R	56.00%	11 Missing :warning:
R/check_zeroinflation.R	69.23%	8 Missing :warning:
R/check_residuals.R	75.00%	7 Missing :warning:
R/check_distribution.R	50.00%	5 Missing :warning:
R/check_normality.R	44.44%	5 Missing :warning:
R/check_overdispersion.R	94.87%	2 Missing :warning:
R/check_outliers.R	95.65%	1 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #643      +/-   ##
==========================================
+ Coverage   57.95%   59.31%   +1.35%     
==========================================
  Files          84       86       +2     
  Lines        6058     6233     +175     
==========================================
+ Hits         3511     3697     +186     
+ Misses       2547     2536      -11     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[mccarthy-m-g]

Thanks! Accepted.

It will probably be a month or two before I can dedicate time to this, but I'll look at the issues you linked above today.

[strengejacke]

I don't think you need to look into all those issues I linked to, that was more a reminder for myself ;-) We should focus on #595, and then I'll check to which extent other issues might be resolved.

[mccarthy-m-g]

Sounds good, although I already started looking... 😅 (DHARMa doesn't support objects of class "lrm", so you can drop #471 from consideration)

I think the performance::check_model.performance_simres method that will be added could probably address some of the issues indirectly (e.g., #501), but that could be dealt with after this PR is merged.

[strengejacke]

but that could be dealt with after this PR is merged

Agreed.

[strengejacke]

So the basic workflow would be:

1. out <- simulate_residuals(model)
2. check_model(out) or check_<whatever>(out)?

Then we wouldn't need particular plot() or print() methods, except for some additional functions (like the suggested check_residuals()).

[mccarthy-m-g]

So the basic workflow would be:

1. out <- simulate_residuals(model)
2. check_model(out) or check_<whatever>(out)?

Then we wouldn't need particular plot() or print() methods, except for some additional functions (like the suggested check_residuals()).

Yes and no. That's the basic workflow, but we'll want a print method for performance_simres and the check_() functions (at least under my vision for how this should get implemented), and we don't want to return the residuals directly in simulate_residuals()

I'll add some rough commits with comments for context.

[mccarthy-m-g]

If you look at the changes from those commits, particularly the vignette I added, does my rationale for why we don't want to prematurely return residuals in simulate_residuals() make sense?

[strengejacke]

If you look at the changes from those commits, particularly the vignette I added, does my rationale for why we don't want to prematurely return residuals in simulate_residuals() make sense?

I see, yes, makes sense!

I'm note sure how to align the plot() method. For testing purposes, I opened a PR (https://github.com/easystats/see/pull/312), but the plots look rather different.

library(performance)
library(glmmTMB)

model <- glmmTMB(
  count ~ mined + spp + (1 | site),
  family = poisson,
  data = Salamanders
)

res <- simulate_residuals(model)
x <- check_normality(res)
plot(x)

plot(DHARMa::simulateResiduals(model))
#> DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

^{Created on 2023-10-27 with reprex v2.0.2}

[mccarthy-m-g]

I'm note sure how to align the plot() method. For testing purposes, I opened a PR, but the plots look rather different.

DHARMa uses a uniform distribution, not a normal distribution, for its left plot. CRAN is down right now and I don't have qqplotr installed so I can't share a reprex, but this code should be what we want as the plot method for check_residuals():

library(glmmTMB)

model <- glmmTMB(
  count ~ mined + spp + (1 | site),
  family = poisson,
  data = Salamanders
)

simulated_residuals <- simulate_residuals(model)

plot_simulated_residuals <- function(x) {

  dp <- list(min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)
  ggplot2::ggplot(
    tibble::tibble(scaled_residuals = residuals(x)),
    ggplot2::aes(sample = scaled_residuals)
  ) +
    qqplotr::stat_qq_band(distribution = "unif", dparams = list(min = 0, max = 1), alpha = .2) +
    qqplotr::stat_qq_line(distribution = "unif", dparams = dp, size = .8, colour = "#3aaf85") +
    qqplotr::stat_qq_point(distribution = "unif", dparams = dp, size = .5, alpha = .8, colour = "#1b6ca8") +
    ggplot2::labs(
      title = "Uniformity of Residuals",
      subtitle = "Dots should fall along the line",
      x = "Standard Uniform Distribution Quantiles",
      y = "Sample Quantiles"
    ) +
    see::theme_lucid()
}

plot_simulated_residuals(simulated_residuals)

I'm a bit skeptical about having a check_normality() method here too, unless we restrict it to the models DHARMa supports where that would make sense to do (e.g., lm). Otherwise it might be better to just return a message and NULL value pointing the user to the check_residuals() function.

[bwiernik]

I would suggest we do check_residuals(model) as the workflow so that we can have correct titles for the plot instead of "Normality of residuals"

[strengejacke]

ok, so check_residuals() would be an additional check we add to the existing functions, right? What about check_model()? Which residual plots do we want there? Or should we add the check_residuals() plot as additional panel? But then it might be confusing when we have different plots checking for "normality" of residuals (i.e. we have two different qq-plots).

[mccarthy-m-g]

ok, so check_residuals() would be an additional check we add to the existing functions, right?

Yeah. We could discuss using check_normality() on simulated residuals later (since I think these residuals can be used for that for models where normality checks are a thing), but for now the main focus should be the check_residuals() implementation since that’s the main feature add for this PR.

What about check_model()? Which residual plots do we want there?

My idea for this was that you would run check_model() on the object returned by simulate_residuals() rather than the model object itself, and the residuals plot would be the uniform distribution one DHARMa uses, with an appropriate title that doesn’t mention “normality” anywhere. So using simulated residuals would be more of an opt-in workflow than an opinionated default. I think this is a good starting point for the PR.

After that basic workflow is implemented I think it would be easier to figure out whether or not this should be the only way to get simulated residual plots, or if they should be also be available through something like check_model(model, residual_type = “simulated”).

[bwiernik]

ok, so check_residuals() would be an additional check we add to the existing functions, right? What about check_model()? Which residual plots do we want there? Or should we add the check_residuals() plot as additional panel? But then it might be confusing when we have different plots checking for "normality" of residuals (i.e. we have two different qq-plots).

For non-normal Gaussian models, we would replace the normality and plot with the residuals plot

[strengejacke]

library(performance)
library(glmmTMB)
#> Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
#> glmmTMB was built with TMB version 1.9.9
#> Current TMB version is 1.9.10
#> Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)

model <- glmmTMB(
  count ~ mined + spp + (1 | site),
  family = poisson,
  data = Salamanders
)

check_model(model)
#> `check_outliers()` does not yet support models of class `glmmTMB`.

check_model(model, detrend = FALSE)
#> `check_outliers()` does not yet support models of class `glmmTMB`.

^{Created on 2024-03-15 with reprex v2.1.0}

[strengejacke]

library(performance)
library(glmmTMB)
#> Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
#> glmmTMB was built with TMB version 1.9.9
#> Current TMB version is 1.9.10
#> Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)

model <- glmmTMB(
  count ~ mined + spp + (1 | site),
  family = poisson,
  data = Salamanders
)

plot(simulate_residuals(model))

plot(simulate_residuals(model), detrend = TRUE)

^{Created on 2024-03-15 with reprex v2.1.0}

[strengejacke]

Not sure if we should set a different range (y limits) for the y axis when detrend = TRUE?

[strengejacke]

@bwiernik and @mccarthy-m-g - I think the implementation works quite well now for the first methods. The plot() method works fine for the Q-Q plots (https://github.com/easystats/see/pull/329). Things we should consider:

• We would have to update the plot for overdispersion/zero-inflation checks (see https://github.com/easystats/performance/pull/643#issuecomment-1999544659). This is still based on the classical residuals, not the simulated ones. I have played around with revising the current code, copying the function .diag_overdispersion into a new .new_diag_overdispersion (see https://github.com/easystats/performance/blob/ad60db15550c406e533ad037107c81d81a66590d/R/check_model_diagnostics.R#L296). This did not really work - do you have any ideas how we can have new plots for overdispersion/zero-inflation? I found it quite informative.

• check_zeroinflation() and check_overdispersion() now rely on simulate_residuals() for zero-inflated or negative binomial models etc., so only the really "simple" models that returned identical results to DHARMa-tests still use the old code. Only Poisson mixed models also use the "old" code, which still could be inaccurate (see #464) - do we want to use simulate_residuals() in general for mixed models when calling check_zeroinflation() and check_overdispersion()?

• How do we want to deal with y axis limits when detrend = TRUE in Q-Q plots? (see https://github.com/easystats/performance/pull/643#issuecomment-1999548863)

[strengejacke]

Here's the implementation for check_outliers():

library(performance)

# For statistical models ---------------------------------------------
# select only mpg and disp (continuous)
mt1 <- mtcars[, c(1, 3, 4)]
# create some fake outliers and attach outliers to main df
mt2 <- rbind(mt1, data.frame(
  mpg = c(37, 40), disp = c(300, 400),
  hp = c(110, 120)
))
# fit model with outliers
model <- lm(disp ~ mpg + hp, data = mt2)

res <- simulate_residuals(model)
check_outliers(res)
#> # Outliers detection
#> 
#>   Proportion of observed outliers: 2.94%
#>   Proportion of expected outliers: 0.80%, 95% CI [0.07, 15.33]
#> No outliers were detected (p = 0.238).
DHARMa::testOutliers(res, plot = FALSE)
#> 
#>  DHARMa outlier test based on exact binomial test with approximate
#>  expectations
#> 
#> data:  res
#> outliers at both margin(s) = 1, observations = 34, p-value = 0.2381
#> alternative hypothesis: true probability of success is not equal to 0.007968127
#> 95 percent confidence interval:
#>  0.0007443642 0.1532676696
#> sample estimates:
#> frequency of outliers (expected: 0.00796812749003984 ) 
#>                                             0.02941176

DHARMa::testOutliers(res, type = "bootstrap", plot = FALSE)
#> 
#>  DHARMa bootstrapped outlier test
#> 
#> data:  res
#> outliers at both margin(s) = 1, observations = 34, p-value = 0.78
#> alternative hypothesis: two.sided
#>  percent confidence interval:
#>  0.00000000 0.05882353
#> sample estimates:
#> outlier frequency (expected: 0.015 ) 
#>                           0.02941176
check_outliers(res, type = "bootstrap", iterations = 200)
#> # Outliers detection
#> 
#>   Proportion of observed outliers: 2.94%
#>   Proportion of expected outliers: 1.32%, 95% CI [0.00, 5.88]
#> No outliers were detected (p = 0.690).

^{Created on 2024-03-18 with reprex v2.1.0}

[strengejacke]

The more I dive into the DHARMa topic, the more I realize that this approach is quite similar to what we aimed at with check_predictions().