Add unit tests for benchmarking metrics and document metric_R2 fallback

modules/benchmark/R/metric_R2.R

@@ -1,23 +1,39 @@
 ##' @name metric_R2
 ##' @title Coefficient of Determination (R2)
 ##' @export
-##' @param metric_dat dataframe
+##' @param metric_dat dataframe with columns \code{model} and \code{obvs}
 ##' @param ... ignored
-##' 
+##'
+##' @details
+##' Computes R-squared using the correlation-based formula:
+##' \eqn{R^2 = \left(\frac{\sum(obs - \bar{obs})(mod - \bar{mod})}
+##' {\sqrt{\sum(obs - \bar{obs})^2} \cdot \sqrt{\sum(mod - \bar{mod})^2}}\right)^2}
+##'
+##' If this formula returns \code{NA} (e.g. when model output is constant
+##' across all observations), the function silently falls back to an
+##' \code{lm()}-based R-squared via \code{summary(lm())$r.squared}.
+##' This fallback may produce unreliable results and triggers a warning
+##' from \code{stats::summary.lm}: "essentially perfect fit: summary may
+##' be unreliable". Consider checking for constant model output before
+##' calling this function.
+##'
 ##' @author Betsy Cowdery
-
 metric_R2 <- function(metric_dat, ...) {
   PEcAn.logger::logger.info("Metric: Coefficient of Determination (R2)")
   numer <- sum((metric_dat$obvs - mean(metric_dat$obvs)) * (metric_dat$model - mean(metric_dat$model)))
   denom <- sqrt(sum((metric_dat$obvs - mean(metric_dat$obvs)) ^ 2)) * sqrt(sum((metric_dat$model - mean(metric_dat$model)) ^ 2))
 
   out <- (numer / denom) ^ 2
 
+  # If correlation formula returns NA (e.g. constant model output),
+  # fall back to lm()-based R-squared. Note: this fallback may trigger
+  # "essentially perfect fit" warning from stats::summary.lm and
+  # produce unreliable results in edge cases.
   if(is.na(out)){
     fit <- stats::lm(metric_dat$model ~ metric_dat$obvs)
     out <- summary(fit)$r.squared
   }
 
   return(out)
 
-} # metric_R2
+} # metric_R2

modules/benchmark/tests/testthat/test-metrics.R

@@ -0,0 +1,49 @@
+# at the top of test-metrics.R, add this:
+if (!requireNamespace("PEcAn.logger", quietly = TRUE)) {
+  PEcAn.logger <- new.env()
+  PEcAn.logger$logger.info <- function(...) invisible(NULL)
+}
+
+test_that("metric_RMSE returns 0 for perfect predictions", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_RMSE(dat), 0)
+})
+
+test_that("metric_RMSE handles NA values", {
+  dat <- data.frame(model = c(1, NA, 3), obvs = c(1, 2, 3))
+  expect_true(is.numeric(metric_RMSE(dat)))
+})
+
+test_that("metric_RMSE returns numeric", {
+  dat <- data.frame(model = c(2, 4), obvs = c(1, 3))
+  expect_equal(metric_RMSE(dat), 1)
+})
+
+test_that("metric_MAE returns 0 for perfect predictions", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_MAE(dat), 0)
+})
+
+test_that("metric_MAE returns correct value", {
+  dat <- data.frame(model = c(3, 3), obvs = c(1, 1))
+  expect_equal(metric_MAE(dat), 2)
+})
+
+test_that("metric_cor returns 1 for perfect linear relationship", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_cor(dat), 1)
+})
+
+test_that("metric_R2 returns 1 for perfect predictions", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_R2(dat), 1)
+})
+
+test_that("metric_R2 silent NA fallback produces valid output", {
+  dat <- data.frame(model = c(2, 2, 2), obvs = c(1, 2, 3))
+  expect_warning(
+    result <- metric_R2(dat),
+    "essentially perfect fit"
+  )
+  expect_true(is.numeric(result))
+})