diff --git a/R/get_predicted_ci.R b/R/get_predicted_ci.R
index 6bcc019e4..e1bf39332 100644
--- a/R/get_predicted_ci.R
+++ b/R/get_predicted_ci.R
@@ -253,6 +253,7 @@ get_predicted_ci.bracl <- get_predicted_ci.mlm
                                     se = NULL,
                                     ci = 0.95,
                                     ci_method = "wald",
+                                    ci_type = "confidence",
                                     data = NULL,
                                     ...) {
   # TODO: Prediction interval for binomial: https://fromthebottomoftheheap.net/2017/05/01/glm-prediction-intervals-i/
@@ -308,6 +309,16 @@ get_predicted_ci.bracl <- get_predicted_ci.mlm
       }
     }
 
+    # add sigma to standard errors, i.e. confidence or prediction intervals
+    ci_type <- validate_argument(ci_type, c("confidence", "prediction"))
+    if (ci_type == "prediction") {
+      if (is_mixed_model(x)) {
+        se <- sqrt(se^2 + get_variance_residual(x))
+      } else {
+        se <- sqrt(se^2 + get_sigma(x)^2)
+      }
+    }
+
     ci_low <- as.numeric(predictions - (se * crit_val))
     ci_high <- as.numeric(predictions + (se * crit_val))
   }
diff --git a/R/get_predicted_mixed.R b/R/get_predicted_mixed.R
index 46182c685..751e275a1 100644
--- a/R/get_predicted_mixed.R
+++ b/R/get_predicted_mixed.R
@@ -31,11 +31,12 @@ get_predicted.lmerMod <- function(x,
     random.only <- FALSE
   }
 
-  # Prediction function
-  predict_function <- function(x, ...) {
+  # Prediction function - we need to pass some arguments, because these
+  # are also required for bootstrapping
+  predict_function <- function(x, data, ...) {
     stats::predict(
       x,
-      newdata = my_args$data,
+      newdata = data,
       type = my_args$type,
       re.form = my_args$re.form,
       random.only = random.only,
@@ -46,7 +47,14 @@ get_predicted.lmerMod <- function(x,
 
   # 1. step: predictions
   if (is.null(iterations)) {
-    predictions <- predict_function(x)
+    rez <- suppressWarnings(predict_function(x, data = my_args$data, se.fit = !is.null(ci)))
+    if (is.list(rez)) {
+      predictions <- as.numeric(rez$fit)
+      se <- rez$se.fit
+    } else {
+      predictions <- as.numeric(rez)
+      se <- NULL
+    }
   } else {
     predictions <- .get_predicted_boot(
       x,
@@ -59,9 +67,14 @@ get_predicted.lmerMod <- function(x,
   }
 
   # 2. step: confidence intervals
-  ci_data <- get_predicted_ci(x, predictions,
-    data = my_args$data, ci = ci,
-    ci_method = ci_method, ci_type = my_args$ci_type,
+  ci_data <- .get_predicted_se_to_ci(
+    x,
+    predictions = predictions,
+    se = se,
+    ci = ci,
+    ci_method = ci_method,
+    ci_type = my_args$ci_type,
+    verbose = verbose,
     ...
   )
 
@@ -113,7 +126,8 @@ get_predicted.glmmTMB <- function(x,
     ...
   )
 
-  # Prediction function
+  # Prediction function - we need to pass some arguments, because these
+  # are also required for bootstrapping
   predict_function <- function(x, data, ...) {
     stats::predict(
       x,
@@ -126,9 +140,8 @@ get_predicted.glmmTMB <- function(x,
   }
 
   # 1. step: predictions
-  rez <- predict_function(x, data = my_args$data, se.fit = TRUE)
-
   if (is.null(iterations)) {
+    rez <- predict_function(x, data = my_args$data, se.fit = TRUE)
     predictions <- as.numeric(rez$fit)
   } else {
     predictions <- .get_predicted_boot(
diff --git a/tests/testthat/test-get_predicted.R b/tests/testthat/test-get_predicted.R
index 64db6d66f..ad9cf3243 100644
--- a/tests/testthat/test-get_predicted.R
+++ b/tests/testthat/test-get_predicted.R
@@ -261,7 +261,7 @@ test_that("get_predicted - glmer with matrix response", {
   grid <- get_datagrid(model, by = "period", range = "grid", preserve_range = FALSE)
   out <- as.data.frame(get_predicted(model, data = grid, ci = 0.95))
   expect_equal(out$Predicted, c(0.19808, 0.08392, 0.07402, 0.04843), tolerance = 1e-3)
-  expect_equal(out$CI_low, c(0.1357, 0.04775, 0.0404, 0.02164), tolerance = 1e-3)
+  expect_equal(out$CI_low, c(0.13647, 0.047502, 0.040128, 0.021335), tolerance = 1e-3)
 })
 
 
diff --git a/tests/testthat/test-glmgee.R b/tests/testthat/test-glmgee.R
index ced96672c..4631b211e 100644
--- a/tests/testthat/test-glmgee.R
+++ b/tests/testthat/test-glmgee.R
@@ -115,7 +115,7 @@ test_that("get_varcov", {
 
 test_that("get_predicted", {
   grid <- get_datagrid(m1, "days")
-  out <- get_predicted(m1, data = grid, ci = 0.95)
+  out <- suppressWarnings(get_predicted(m1, data = grid, ci = 0.95))
   expect_equal(
     as.numeric(out),
     c(