diff --git a/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py b/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py
index c6ca30af9..925f80823 100644
--- a/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py
+++ b/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py
@@ -19,13 +19,15 @@
 
 import numpy as np
 import numpy.typing as npt
-from sklearn.metrics import r2_score
 
 from openstef_core.types import Quantile
 
 _Q_05 = Quantile(0.05)
 _Q_95 = Quantile(0.95)
 
+# R² is undefined for fewer than this many samples (variance needs at least two points).
+_MIN_R2_SAMPLES = 2
+
 
 def completeness(
     y: npt.NDArray[np.floating],
@@ -527,7 +529,8 @@ def r2(
     Returns:
         The R² score as a float. Best possible score is 1.0, and it can be negative
         (because the model can be arbitrarily worse). A constant model that always
-        predicts the mean of y_true would get an R² score of 0.0.
+        predicts the mean of y_true would get an R² score of 0.0. Fewer than two
+        samples returns NaN, since R² is undefined there (matching scikit-learn).
 
     Example:
         Basic usage with energy load data
@@ -552,10 +555,24 @@ def r2(
         >>> isinstance(score, float)
         True
     """
-    if len(y_true) == 0 or len(y_pred) == 0:
+    # R² is undefined for fewer than two samples; match scikit-learn, which
+    # returns NaN (with a warning) in that case.
+    if len(y_true) < _MIN_R2_SAMPLES or len(y_pred) < _MIN_R2_SAMPLES:
         return float("NaN")
 
-    return float(r2_score(y_true, y_pred, sample_weight=sample_weights))
+    y_true = np.asarray(y_true, dtype=float)
+    y_pred = np.asarray(y_pred, dtype=float)
+    weights = np.ones_like(y_true) if sample_weights is None else np.asarray(sample_weights, dtype=float)
+
+    weighted_mean = np.average(y_true, weights=weights)
+    residual_sum = float(np.sum(weights * (y_true - y_pred) ** 2))
+    total_sum = float(np.sum(weights * (y_true - weighted_mean) ** 2))
+
+    # Constant y_true: match scikit-learn (1.0 for a perfect fit, else 0.0).
+    if total_sum == 0.0:
+        return 1.0 if residual_sum == 0.0 else 0.0
+
+    return float(1.0 - residual_sum / total_sum)
 
 
 def pinball_loss(
diff --git a/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py b/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py
index d73f069af..4028bc87f 100644
--- a/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py
+++ b/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py
@@ -15,6 +15,7 @@
     mape,
     pinball_loss,
     precision_recall,
+    r2,
     relative_pinball_loss,
     riqd,
     rmae,
@@ -624,3 +625,30 @@ def test_pinball_loss_various(
 
     # Assert
     assert abs(result - expected) < 1e-8, f"Expected {expected} but got {result}"
+
+
+def test_r2_perfect_and_constant_predictor():
+    """R² is 1.0 for a perfect fit and 0.0 for a constant mean predictor."""
+    y_true = np.array([1.0, 2.0, 3.0, 4.0])
+
+    assert r2(y_true, y_true) == pytest.approx(1.0)
+    assert r2(y_true, np.full_like(y_true, y_true.mean())) == pytest.approx(0.0)
+
+
+def test_r2_matches_known_values():
+    """R² matches hand-computed values, including a sample-weighted case."""
+    # 1 - 1.5 / 29.1875
+    assert r2(np.array([3.0, -0.5, 2.0, 7.0]), np.array([2.5, 0.0, 2.0, 8.0])) == pytest.approx(0.948608, abs=1e-5)
+    # Weighted: weighted mean 2.25, residual sum 2.0, total sum 2.75 -> 1 - 2 / 2.75
+    weighted = r2(
+        np.array([1.0, 2.0, 3.0]),
+        np.array([1.0, 2.0, 4.0]),
+        sample_weights=np.array([1.0, 1.0, 2.0]),
+    )
+    assert weighted == pytest.approx(0.272727, abs=1e-5)
+
+
+def test_r2_undefined_for_fewer_than_two_samples():
+    """Fewer than two samples returns NaN (R² undefined, matches scikit-learn)."""
+    assert np.isnan(r2(np.array([]), np.array([])))
+    assert np.isnan(r2(np.array([5.0]), np.array([5.0])))