Use standard errors from predict-method for lmerMod

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))
   }

R/get_predicted_mixed.R

@@ -31,11 +31,12 @@
     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 @@
 
   # 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 @@
   }
 
   # 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,
     ...
   )
 
@@ -93,14 +106,14 @@
     predict <- "expectation"
     if (verbose) {
       format_warning(
         "\"prediction\" and \"classification\" are currently not supported by the `predict` argument for `glmmTMB` models.",
         "Changing to `predict=\"expectation\"`."
       )
     }
   }
 
   # TODO: prediction intervals
   # https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#predictions-andor-confidence-or-prediction-intervals-on-predictions
 
   # Sanitize input
   my_args <- .get_predicted_args(
@@ -113,7 +126,8 @@
     ...
   )
 
-  # 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 @@
   }
 
   # 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(
@@ -206,7 +219,7 @@
     predict <- "expectation"
     if (verbose) {
       format_warning(
         "\"prediction\" and \"classification\" are currently not supported by the `predict` argument for `GLMMadaptive` models.",
         "Changing to `predict=\"expectation\"`."
       )
     }
@@ -255,7 +268,7 @@
 
   if (my_args$scale == "response" && my_args$info$is_zero_inflated) {
     # 2. and 3. step: confidence intervals and back-transform
     ci_data <- .simulate_zi_predictions(model = x, newdata = data, predictions = predictions, nsim = iterations, ci = ci)
     out <- list(predictions = predictions, ci_data = ci_data)
   } else {
     # 2. step: confidence intervals

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)
 })
 
 

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(