diff --git a/R/binned_residuals.R b/R/binned_residuals.R
index 286323968..2614b67a6 100644
--- a/R/binned_residuals.R
+++ b/R/binned_residuals.R
@@ -62,7 +62,15 @@
 #' }
 #'
 #' @export
-binned_residuals <- function(model, term = NULL, n_bins = NULL, show_dots = NULL, ...) {
+binned_residuals <- function(model,
+                             term = NULL,
+                             n_bins = NULL,
+                             show_dots = NULL,
+                             residuals = c("gaussian", "exact", "boot"),
+                             ci = 0.95,
+                             ...) {
+
+  residuals <- match.arg(residuals)
   fv <- stats::fitted(model)
   mf <- insight::get_data(model, verbose = FALSE)
 
@@ -78,7 +86,8 @@ binned_residuals <- function(model, term = NULL, n_bins = NULL, show_dots = NULL
     show_dots <- is.null(n) || n <= 1e5
   }
 
-  y <- .recode_to_zero(insight::get_response(model, verbose = FALSE)) - fv
+  y0 <- .recode_to_zero(insight::get_response(model, verbose = FALSE))
+  y <- y0 - fv
 
   if (is.null(n_bins)) n_bins <- round(sqrt(length(pred)))
 
@@ -95,24 +104,28 @@ binned_residuals <- function(model, term = NULL, n_bins = NULL, show_dots = NULL
     n <- length(items)
     sdev <- stats::sd(y[items], na.rm = TRUE)
 
-    data.frame(
+    r <- switch(residuals,
+      gaussian = stats::qnorm(c((1 - ci) / 2, (1 + ci) / 2), mean = ybar, sd = sdev / sqrt(n)),
+      exact = stats:::binom.test(sum(y0[items]), n)$conf.int - fv,
+      boot = Hmisc::smean.cl.boot(y[items], conf.int = ci)[c("Lower", "Upper")]
+    )
+    names(r) <- c("CI_low", "CI_high")
+
+    d0 <- data.frame(
       xbar = xbar,
       ybar = ybar,
       n = n,
       x.lo = model.range[1],
       x.hi = model.range[2],
-      se = stats::qnorm(0.975) * sdev / sqrt(n),
+      se = stats::qnorm((1 + ci) / 2) * sdev / sqrt(n),
       ci_range = sdev / sqrt(n)
     )
+    cbind(d0, rbind(r))
   }))
 
   d <- do.call(rbind, d)
   d <- d[stats::complete.cases(d), ]
 
-  # CIs
-  d$CI_low <- d$ybar - stats::qnorm(0.975) * d$ci_range
-  d$CI_high <- d$ybar + stats::qnorm(0.975) * d$ci_range
-
   gr <- abs(d$ybar) > abs(d$se)
   d$group <- "yes"
   d$group[gr] <- "no"