diff --git a/R/check_model.R b/R/check_model.R
index 853e647a6..e763b0fe9 100644
--- a/R/check_model.R
+++ b/R/check_model.R
@@ -218,7 +218,7 @@ check_model.model_fit <- function(x,
 
 
 
-# helper ------------------------
+# compile plots for checks of linear models  ------------------------
 
 .check_assumptions_linear <- function(model, model_info) {
   dat <- list()
@@ -245,6 +245,8 @@ check_model.model_fit <- function(x,
 
 
 
+# compile plots for checks of generalized linear models  ------------------------
+
 .check_assumptions_glm <- function(model, model_info) {
   dat <- list()
 
@@ -263,6 +265,9 @@ check_model.model_fit <- function(x,
   if (isTRUE(model_info$is_binomial)) {
     dat$BINNED_RESID <- binned_residuals(model)
   }
+  if (isTRUE(model_info$is_count)) {
+    dat$OVERDISPERION <- .diag_overdispersion(model)
+  }
 
   dat <- datawizard::compact_list(dat)
   class(dat) <- c("check_model", "see_check_model")
@@ -271,6 +276,8 @@ check_model.model_fit <- function(x,
 
 
 
+# compile plots for checks of Bayesian models  ------------------------
+
 .check_assumptions_stan <- function(model) {
   if (inherits(model, "brmsfit")) {
 
diff --git a/R/check_model_diagnostics.R b/R/check_model_diagnostics.R
index a9c09faec..569bda816 100644
--- a/R/check_model_diagnostics.R
+++ b/R/check_model_diagnostics.R
@@ -1,3 +1,5 @@
+# prepare data for VIF plot ----------------------------------
+
 .diag_vif <- function(model) {
   dat <- datawizard::compact_list(check_collinearity(model))
   if (is.null(dat)) {
@@ -17,6 +19,9 @@
 }
 
 
+
+# prepare data for QQ plot ----------------------------------
+
 .diag_qq <- function(model) {
   if (inherits(model, c("lme", "lmerMod", "merMod", "glmmTMB", "gam"))) {
     res_ <- sort(stats::residuals(model), na.last = NA)
@@ -54,6 +59,8 @@
 
 
 
+# prepare data for random effects QQ plot ----------------------------------
+
 .diag_reqq <- function(model, level = .95, model_info) {
   # check if we have mixed model
   if (!model_info$is_mixed) {
@@ -117,6 +124,7 @@
 
 
 
+# prepare data for normality of residuals plot ----------------------------------
 
 .diag_norm <- function(model) {
   r <- try(stats::residuals(model), silent = TRUE)
@@ -133,6 +141,8 @@
 
 
 
+# prepare data for influential obs plot ----------------------------------
+
 .diag_influential_obs <- function(model, threshold = NULL) {
   s <- summary(model)
 
@@ -184,6 +194,7 @@
 
 
 
+# prepare data for non-constant variance plot ----------------------------------
 
 .diag_ncv <- function(model) {
   ncv <- tryCatch(
@@ -207,6 +218,9 @@
 }
 
 
+
+# prepare data for homogeneity of variance plot ----------------------------------
+
 .diag_homogeneity <- function(model) {
   faminfo <- insight::model_info(model)
   r <- tryCatch(
@@ -246,6 +260,48 @@
 
 
 
+# prepare data for homogeneity of variance plot ----------------------------------
+
+.diag_overdispersion <- function(model) {
+  faminfo <- insight::model_info(model)
+
+  # data for poisson models
+  if (faminfo$is_poisson && !faminfo$is_zero_inflated) {
+    d <- as.data.frame(insight::get_predicted(model, predict = "expectation", ci = NA))
+    d$Residuals <- insight::get_response(model) - as.vector(d$Predicted)
+    d$Res2 <- d$Residuals^2
+    d$V <- d$Predicted
+  }
+
+  # data for negative binomial models
+  if (faminfo$is_negbin && !faminfo$is_zero_inflated) {
+    d <- as.data.frame(insight::get_predicted(model, predict = "expectation", ci = NA))
+    d$Residuals <- insight::get_response(model) - as.vector(d$Predicted)
+    d$Res2 <- d$Residuals^2
+    d$V <- d$Predicted * (1 + d$Predicted / insight::get_sigma(model))
+  }
+
+  # data for zero-inflated poisson models
+  if (faminfo$is_poisson && faminfo$is_zero_inflated) {
+    d <- as.data.frame(insight::get_predicted(model, predict = "expectation", ci = NA))
+    d$Residuals <- insight::get_response(model) - as.vector(d$Predicted)
+    d$Res2 <- d$Residuals^2
+    if (inherits(model, "glmmTMB")) {
+      ptype <- "zprob"
+    } else {
+      ptype <- "zero"
+    }
+    d$Prob <- stats::predict(model, type = ptype)
+    d$V <- d$Predicted * (1 - d$Prob) * (1 + d$Predicted * d$Prob)
+  }
+
+  d
+}
+
+
+
+# helpers ----------------------------------
+
 .sigma_glmmTMB_nonmixed <- function(model, faminfo) {
   if (!is.na(match(faminfo$family, c("binomial", "poisson", "truncated_poisson")))) {
     return(1)