Take over PR #3888: Add unit tests for benchmarking metrics and document metric_R2

modules/benchmark/R/metric_R2.R

@@ -1,23 +1,20 @@
 ##' @name metric_R2
-##' @title Coefficient of Determination (R2)
+##' @title Squared Pearson Correlation (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}
 ##' 
+##' Note: Because this is a correlation-based R2, it is invariant to bias and slope.
+##' A model that is perfectly correlated but biased (e.g., model = obs + 100) will still
+##' score 1. This is distinct from variance explained or Nash-Sutcliffe Efficiency (NSE).
+##'
 ##' @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(is.na(out)){
-    fit <- stats::lm(metric_dat$model ~ metric_dat$obvs)
-    out <- summary(fit)$r.squared
-  }
-
-  return(out)
-
+  PEcAn.logger::logger.info("Metric: Squared Pearson Correlation (R2)")
+  stats::cor(metric_dat$model, metric_dat$obvs, use = "pairwise.complete.obs")^2
 } # metric_R2

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

@@ -0,0 +1,43 @@
+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_equal(metric_RMSE(dat), 0)
+})
+
+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 returns NA for constant model output", {
+  dat <- data.frame(model = c(2, 2, 2), obvs = c(1, 2, 3))
+  expect_warning(
+    result <- metric_R2(dat),
+    "the standard deviation is zero"
+  )
+  expect_equal(result, NA_real_)
+})