95 CI for rmse

Fix #740

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: performance
 Title: Assessment of Regression Models Performance
-Version: 0.12.0.4
+Version: 0.12.0.5
 Authors@R: 
     c(person(given = "Daniel",
              family = "Lüdecke",

NAMESPACE

@@ -148,6 +148,7 @@ S3method(display,test_performance)
 S3method(fitted,BFBayesFactor)
 S3method(format,compare_performance)
 S3method(format,performance_model)
+S3method(format,performance_rmse)
 S3method(format,test_performance)
 S3method(logLik,cpglm)
 S3method(logLik,iv_robust)
@@ -319,6 +320,7 @@ S3method(print,performance_hosmer)
 S3method(print,performance_model)
 S3method(print,performance_pcp)
 S3method(print,performance_pp_check)
+S3method(print,performance_rmse)
 S3method(print,performance_roc)
 S3method(print,performance_score)
 S3method(print,performance_simres)

NEWS.md

@@ -14,6 +14,10 @@
 * `icc()` and `r2_nakagawa()` get a `model_component` argument indicating the 
   component for zero-inflation or hurdle models.
 
+* `performance_rmse()` (resp. `rmse()`) can now compute analytical and 
+  bootstrapped confidence intervals. The function gains following new arguments:
+  `ci`, `ci_method` and `iterations`.
+
 # performance 0.12.0
 
 ## Breaking

R/icc.R

@@ -1,5 +1,7 @@
-#' Intraclass Correlation Coefficient (ICC)
+#' @title Intraclass Correlation Coefficient (ICC)
+#' @name icc
 #'
+#' @description
 #' This function calculates the intraclass-correlation coefficient (ICC) -
 #' sometimes also called *variance partition coefficient* (VPC) or
 #' *repeatability* - for mixed effects models. The ICC can be calculated for all
@@ -15,32 +17,32 @@
 #' Else, for instance for nested models, name a specific group-level effect
 #' to calculate the variance decomposition for this group-level. See 'Details'
 #' and `?brms::posterior_predict`.
-#' @param ci Confidence resp. credible interval level. For `icc()` and `r2()`,
-#' confidence intervals are based on bootstrapped samples from the ICC resp.
-#' R2 value. See `iterations`.
-#' @param by_group Logical, if `TRUE`, `icc()` returns the variance
-#' components for each random-effects level (if there are multiple levels).
-#' See 'Details'.
+#' @param ci Confidence resp. credible interval level. For `icc()`, `r2()`, and
+#' `rmse()`, confidence intervals are based on bootstrapped samples from the
+#' ICC, R2 or RMSE value. See `iterations`.
+#' @param by_group Logical, if `TRUE`, `icc()` returns the variance components
+#' for each random-effects level (if there are multiple levels). See 'Details'.
 #' @param iterations Number of bootstrap-replicates when computing confidence
-#' intervals for the ICC or R2.
+#' intervals for the ICC, R2, RMSE etc.
 #' @param ci_method Character string, indicating the bootstrap-method. Should
-#' be `NULL` (default), in which case `lme4::bootMer()` is used for
-#' bootstrapped confidence intervals. However, if bootstrapped intervals cannot
-#' be calculated this was, try `ci_method = "boot"`, which falls back to
-#' `boot::boot()`. This may successfully return bootstrapped confidence intervals,
-#' but bootstrapped samples may not be appropriate for the multilevel structure
-#' of the model. There is also an option `ci_method = "analytical"`, which tries
-#' to calculate analytical confidence assuming a chi-squared distribution.
-#' However, these intervals are rather inaccurate and often too narrow. It is
-#' recommended to calculate bootstrapped confidence intervals for mixed models.
+#' be `NULL` (default), in which case `lme4::bootMer()` is used for bootstrapped
+#' confidence intervals. However, if bootstrapped intervals cannot be calculated
+#' this way, try `ci_method = "boot"`, which falls back to `boot::boot()`. This
+#' may successfully return bootstrapped confidence intervals, but bootstrapped
+#' samples may not be appropriate for the multilevel structure of the model.
+#' There is also an option `ci_method = "analytical"`, which tries to calculate
+#' analytical confidence assuming a chi-squared distribution. However, these
+#' intervals are rather inaccurate and often too narrow. It is recommended to
+#' calculate bootstrapped confidence intervals for mixed models.
 #' @param verbose Toggle warnings and messages.
 #' @param null_model Optional, a null model to compute the random effect variances,
 #' which is passed to [`insight::get_variance()`]. Usually only required if
 #' calculation of r-squared or ICC fails when `null_model` is not specified. If
 #' calculating the null model takes longer and you already have fit the null
 #' model, you can pass it here, too, to speed up the process.
 #' @param ... Arguments passed down to `lme4::bootMer()` or `boot::boot()`
-#' for bootstrapped ICC or R2.
+#' for bootstrapped ICC, R2, RMSE etc.; for `variance_decomposition()`,
+#' arguments are passed down to `brms::posterior_predict()`.
 #'
 #' @inheritParams r2_bayes
 #' @inheritParams insight::get_variance
@@ -321,8 +323,6 @@ icc <- function(model,
 
 
 
-#' @param ... Arguments passed down to `brms::posterior_predict()`.
-#' @inheritParams icc
 #' @rdname icc
 #' @export
 variance_decomposition <- function(model,

R/performance_rmse.R

@@ -6,7 +6,7 @@
 #'
 #' @param model A model.
 #' @param normalized Logical, use `TRUE` if normalized rmse should be returned.
-#' @inheritParams model_performance.lm
+#' @inheritParams icc
 #'
 #' @details The RMSE is the square root of the variance of the residuals and indicates
 #'   the absolute fit of the model to the data (difference between observed data
@@ -30,33 +30,138 @@
 #' # normalized RMSE
 #' performance_rmse(m, normalized = TRUE)
 #' @export
-performance_rmse <- function(model, normalized = FALSE, verbose = TRUE) {
+performance_rmse <- function(model,
+                             normalized = FALSE,
+                             ci = NULL,
+                             iterations = 100,
+                             ci_method = NULL,
+                             verbose = TRUE,
+                             ...) {
   tryCatch(
     {
-      # compute rmse
-      rmse_val <- sqrt(performance_mse(model, verbose = verbose))
-
-      # if normalized, divide by range of response
-      if (normalized) {
-        # get response
-        resp <- datawizard::to_numeric(insight::get_response(model, verbose = FALSE), dummy_factors = FALSE, preserve_levels = TRUE)
-        # compute rmse, normalized
-        rmse_val <- rmse_val / (max(resp, na.rm = TRUE) - min(resp, na.rm = TRUE))
+      out <- .calculate_rmse(model, normalized, verbose)
+      # check if CIs are requested, and compute CIs
+      if (!is.null(ci) && !is.na(ci)) {
+        # analytical CI?
+        if (identical(ci_method, "analytical")) {
+          out <- .analytical_rmse_ci(out, model, ci)
+        } else {
+          # bootstrapped CI
+          result <- .bootstrap_rmse(model, iterations, normalized, ci_method, ...)
+          # CI for RMSE
+          rmse_ci <- as.vector(result$t[, 1])
+          rmse_ci <- rmse_ci[!is.na(rmse_ci)]
+          # validation check
+          if (length(rmse_ci) > 0) {
+            rmse_ci <- bayestestR::eti(rmse_ci, ci = ci)
+            out <- cbind(data.frame(RMSE = out), rmse_ci)
+            class(out) <- c("performance_rmse", "data.frame")
+          } else {
+            insight::format_warning("Could not compute confidence intervals for RMSE.")
+          }
+        }
       }
-
-      rmse_val
     },
     error = function(e) {
       if (inherits(e, c("simpleError", "error")) && verbose) {
         insight::print_color(e$message, "red")
         cat("\n")
       }
-      NA
+      out <- NA
     }
   )
+
+  out
 }
 
 
 #' @rdname performance_rmse
 #' @export
 rmse <- performance_rmse
+
+
+# methods ---------------------------------------------------------------------
+
+#' @export
+format.performance_rmse <- function(x, ...) {
+  insight::format_table(x, ...)
+}
+
+#' @export
+print.performance_rmse <- function(x, ...) {
+  cat(insight::export_table(format(x, ...), ...))
+}
+
+
+# helper function to compute RMSE ----------------------------------------------
+
+.calculate_rmse <- function(model, normalized = FALSE, verbose = FALSE, ...) {
+  # compute rmse
+  rmse_val <- sqrt(performance_mse(model, verbose = verbose))
+
+  # if normalized, divide by range of response
+  if (normalized) {
+    # get response
+    resp <- datawizard::to_numeric(
+      insight::get_response(model, verbose = FALSE),
+      dummy_factors = FALSE,
+      preserve_levels = TRUE
+    )
+    # compute rmse, normalized
+    rmse_val <- rmse_val / (max(resp, na.rm = TRUE) - min(resp, na.rm = TRUE))
+  }
+
+  rmse_val
+}
+
+
+# analytical CIs --------------------------------------------------------------
+
+.analytical_rmse_ci <- function(out, model, ci, ...) {
+  s <- insight::get_sigma(model, ci = ci, verbose = FALSE)
+  n <- insight::n_obs(model)
+  conf_ints <- c(attr(s, "CI_low"), attr(s, "CI_high")) * ((n - 1) / n)
+  out <- data.frame(
+    RMSE = out,
+    CI = ci,
+    CI_low = conf_ints[1],
+    CI_high = conf_ints[2]
+  )
+  class(out) <- c("performance_rmse", "data.frame")
+  out
+}
+
+
+# bootstrapping CIs -----------------------------------------------------------
+
+.boot_calculate_rmse <- function(data, indices, model, normalized, ...) {
+  d <- data[indices, ] # allows boot to select sample
+  fit <- suppressWarnings(suppressMessages(stats::update(model, data = d)))
+  .calculate_rmse(model = fit, normalized = normalized)
+}
+
+.bootstrap_rmse <- function(model, iterations = 100, normalized = FALSE, ci_method = NULL, ...) {
+  if (inherits(model, c("merMod", "lmerMod", "glmmTMB")) && !identical(ci_method, "boot")) {
+    # cannot pass argument "normalized" to "lme4::bootMer()"
+    if (isTRUE(normalized)) {
+      insight::format_error("Normalized RMSE cannot be used with confidence intervals. Please use `ci_method = \"boot\"`.") # nolint
+    }
+    result <- .do_lme4_bootmer(
+      model,
+      .calculate_rmse,
+      iterations,
+      dots = list(...)
+    )
+  } else {
+    insight::check_if_installed("boot")
+    result <- boot::boot(
+      data = insight::get_data(model, verbose = FALSE),
+      statistic = .boot_calculate_rmse,
+      R = iterations,
+      sim = "ordinary",
+      model = model,
+      normalized = normalized
+    )
+  }
+  result
+}