add median yearly_r2

Author: michielkallenberg

cybench/evaluation/aggregated_metrics.py

@@ -40,6 +40,26 @@ def calc_nrmse(y_true: npt.ArrayLike, y_pred: npt.ArrayLike) -> float:
     return float("nan") if denom == 0 else float(rmse / denom)
 
 
+def calc_median_yearly_r2(
+    df: pd.DataFrame,
+    target_col: str,
+    model_col: str,
+    *,
+    year_col: str = KEY_YEAR,
+) -> float:
+    """Median of per-year R², where each year's R² is computed across regions."""
+    yearly_r2: list[float] = []
+    for _, year_df in df.groupby(year_col):
+        _, r2 = calc_r_r2(
+            year_df[target_col].values,
+            year_df[model_col].values,
+        )
+        yearly_r2.append(r2)
+    if not yearly_r2:
+        return float("nan")
+    return float(np.nanmedian(yearly_r2))
+
+
 def get_metrics_dict(
     df: pd.DataFrame,
     target_col: str,
@@ -79,13 +99,16 @@ def compute_report_metrics(
     Metrics used by visualize_results_aggregated / build_results_dashboard.
 
     Views:
-      - region_year: pooled region-year rows (r, R², NRMSE, anomaly r/R²)
+      - region_year: pooled region-year rows (r, R², NRMSE, med per-year R², anomaly r/R²)
       - spatial: per-location means across years
-      - temporal: per-year means across locations
+      - temporal: per-year means across locations (r, R²)
     """
     complete = df[target_col].notna() & df[model_col].notna()
     df = df.loc[complete].copy()
     region_year = get_metrics_dict(df, target_col, model_col, loc_col=loc_col)
+    region_year["median_r2"] = calc_median_yearly_r2(
+        df, target_col, model_col, year_col=year_col
+    )
 
     spatial = df.groupby(loc_col)[[target_col, model_col]].mean()
     r_spatial, r2_spatial = calc_r_r2(
@@ -115,7 +138,8 @@ def format_report_metrics(metrics: dict[str, Any]) -> str:
     sp = metrics["spatial"]
     tm = metrics["temporal"]
     return (
-        f"region-year r={ry['r']:.2f} R²={ry['r2']:.2f} NRMSE={ry['nrmse']:.2f} | "
+        f"region-year r={ry['r']:.2f} R²={ry['r2']:.2f} NRMSE={ry['nrmse']:.2f} "
+        f"med-R²/yr={ry['median_r2']:.2f} | "
         f"spatial r={sp['r']:.2f} R²={sp['r2']:.2f} | "
         f"temporal r={tm['r']:.2f} R²={tm['r2']:.2f} | "
         f"anomaly r={ry['r_res']:.2f} R²={ry['r2_res']:.2f}"

cybench/runs/analysis/benchmark_run_catalog.py

@@ -20,6 +20,7 @@
     "r",
     "r2",
     "nrmse",
+    "r2_yearly_median",
     "r_spatial",
     "r2_spatial",
     "r_temporal",
@@ -29,7 +30,19 @@
 )
 
 # For deltas: higher is better vs lower is better.
-HIGHER_IS_BETTER = frozenset({"r", "r2", "r_spatial", "r2_spatial", "r_temporal", "r2_temporal", "r_res", "r2_res"})
+HIGHER_IS_BETTER = frozenset(
+    {
+        "r",
+        "r2",
+        "r_spatial",
+        "r2_spatial",
+        "r_temporal",
+        "r2_temporal",
+        "r2_yearly_median",
+        "r_res",
+        "r2_res",
+    }
+)
 LOWER_IS_BETTER = frozenset({"nrmse"})
 
 
@@ -44,6 +57,7 @@ def flatten_report_metrics(metrics: dict[str, Any]) -> dict[str, Any]:
         "r": ry.get("r"),
         "r2": ry.get("r2"),
         "nrmse": ry.get("nrmse"),
+        "r2_yearly_median": ry.get("median_r2"),
         "r_res": ry.get("r_res"),
         "r2_res": ry.get("r2_res"),
         "r_spatial": sp.get("r"),

cybench/runs/analysis/collect_walk_forward_results.py

@@ -162,6 +162,7 @@ def summary_rows_to_dashboard_records(
         ("region_year", "r", "r"),
         ("region_year", "r2", "r2"),
         ("region_year", "nrmse", "nrmse"),
+        ("region_year", "median_r2", "r2_yearly_median"),
         ("spatial", "r", "r_spatial"),
         ("spatial", "r2", "r2_spatial"),
         ("temporal", "r", "r_temporal"),

cybench/runs/viz/build_results_dashboard.py

@@ -13,6 +13,8 @@
 import numpy as np
 import pandas as pd
 
+from cybench.evaluation.aggregated_metrics import calc_median_yearly_r2
+
 
 @dataclass
 class SourceConfig:
@@ -152,6 +154,7 @@ def load_records_from_runs_root(
                 ("region_year", "r", safe_float(m.get("r"))),
                 ("region_year", "r2", safe_float(m.get("r2"))),
                 ("region_year", "nrmse", safe_float(m.get("nrmse"))),
+                ("region_year", "median_r2", safe_float(s.get("r2_yearly_median"))),
                 ("spatial", "r", safe_float(sp.get("r"))),
                 ("spatial", "r2", safe_float(sp.get("r2"))),
                 ("temporal", "r", safe_float(s.get("r_time_model"))),
@@ -198,6 +201,7 @@ def load_records_from_runs_root(
 
         temporal = clean.groupby("year")[["yield", "pred"]].mean().sort_index()
         r_tm, r2_tm = calc_r_r2(temporal["yield"].values, temporal["pred"].values)
+        r2_yearly_median = calc_median_yearly_r2(clean, "yield", "pred", year_col="year")
 
         loc_means = clean.groupby("adm_id")["yield"].mean()
         y_res = clean["yield"] - clean["adm_id"].map(loc_means)
@@ -223,6 +227,7 @@ def load_records_from_runs_root(
             ("region_year", "r", safe_float(r)),
             ("region_year", "r2", safe_float(r2)),
             ("region_year", "nrmse", safe_float(nrmse)),
+            ("region_year", "median_r2", safe_float(r2_yearly_median)),
             ("spatial", "r", safe_float(r_sp)),
             ("spatial", "r2", safe_float(r2_sp)),
             ("temporal", "r", safe_float(r_tm)),
@@ -339,6 +344,7 @@ def load_records(sources: List[SourceConfig], output_dir: str) -> List[Dict]:
                 ("region_year", "r", safe_float(m.get("r"))),
                 ("region_year", "r2", safe_float(m.get("r2"))),
                 ("region_year", "nrmse", safe_float(m.get("nrmse"))),
+                ("region_year", "median_r2", safe_float(s.get("r2_yearly_median"))),
                 ("spatial", "r", safe_float(sp.get("r"))),
                 ("spatial", "r2", safe_float(sp.get("r2"))),
                 ("temporal", "r", safe_float(s.get("r_time_model"))),

cybench/runs/viz/dashboard_template.html

@@ -250,6 +250,7 @@ <h3>Drill-down plots</h3>
       { key: "region_year_r", label: "r", group: "Region-Year", higherBetter: true },
       { key: "region_year_r2", label: "R²", group: "Region-Year", higherBetter: true },
       { key: "region_year_nrmse", label: "NRMSE", group: "Region-Year", higherBetter: false },
+      { key: "region_year_median_r2", label: "med R²/yr", group: "Region-Year", higherBetter: true },
       { key: "spatial_r", label: "r", group: "Spatial", higherBetter: true },
       { key: "spatial_r2", label: "R²", group: "Spatial", higherBetter: true },
       { key: "temporal_r", label: "r", group: "Temporal", higherBetter: true },
@@ -263,6 +264,7 @@ <h3>Drill-down plots</h3>
       region_year_r: { min: 0, max: 1, higherBetter: true },
       region_year_r2: { min: 0, max: 1, higherBetter: true },
       region_year_nrmse: { min: 0, max: 0.5, higherBetter: false },
+      region_year_median_r2: { min: 0, max: 1, higherBetter: true },
       spatial_r: { min: 0, max: 1, higherBetter: true },
       spatial_r2: { min: 0, max: 1, higherBetter: true },
       temporal_r: { min: 0, max: 1, higherBetter: true },
@@ -425,7 +427,7 @@ <h3>Drill-down plots</h3>
       <th rowspan="2" class="sortable" data-sort="n_regions">Regions ${sortMark("n_regions")}</th>
       <th rowspan="2" class="sortable" data-sort="n_years">Years ${sortMark("n_years")}</th>
       <th rowspan="2" class="sortable" data-sort="n_samples">Samples ${sortMark("n_samples")}</th>
-      <th colspan="3" class="group-head">Region-Year</th>
+      <th colspan="4" class="group-head">Region-Year</th>
       <th colspan="2" class="group-head">Spatial</th>
       <th colspan="2" class="group-head">Temporal</th>
       <th colspan="2" class="group-head">Anomaly</th>

cybench/runs/viz/visualize_results_aggregated.py

@@ -20,7 +20,11 @@
 
 from cybench.util.geo import get_shapes_from_polygons, world_shape_path
 from cybench.config import KEY_LOC, KEY_TARGET
-from cybench.evaluation.aggregated_metrics import calc_r_r2, get_metrics_dict
+from cybench.evaluation.aggregated_metrics import (
+    calc_median_yearly_r2,
+    calc_r_r2,
+    get_metrics_dict,
+)
 
 # -----------------------------
 # Configuration
@@ -208,21 +212,21 @@ def generate_markdown_table(stats_list: List[dict]) -> str:
 
     md = (
         "| Dataset | N_regions | N_years | "
-        "Region-Year | Region-Year | Region-Year | "
+        "Region-Year | Region-Year | Region-Year | Region-Year | "
         "Spatial | Spatial | "
         "Temporal | Temporal | "
         "Anomaly | Anomaly |\n"
     )
     md += (
         "|  |  |  | "
-        "r | R² | NRMSE | "
+        "r | R² | NRMSE | med R²/yr | "
         "r | R² | "
         "r | R² | "
         "r | R² |\n"
     )
     md += (
         "| :--- | ---: | ---: | "
-        "---: | ---: | ---: | "
+        "---: | ---: | ---: | ---: | "
         "---: | ---: | "
         "---: | ---: | "
         "---: | ---: |\n"
@@ -255,6 +259,7 @@ def style(text):
             f"| {style(fmt(mod['r']))} "
             f"| {style(fmt(mod['r2']))} "
             f"| {style(fmt(mod['nrmse']))} "
+            f"| {style(fmt(s['r2_yearly_median']))} "
             f"| {style(fmt(sp['r']))} "
             f"| {style(fmt(sp['r2']))} "
             f"| {style(fmt(s['r_time_model']))} "
@@ -293,13 +298,13 @@ def fmt(val):
       <th class="left" rowspan="2">Dataset</th>
       <th rowspan="2">N_regions</th>
       <th rowspan="2">N_years</th>
-      <th colspan="3">Region-Year</th>
+      <th colspan="4">Region-Year</th>
       <th colspan="2">Spatial</th>
       <th colspan="2">Temporal</th>
       <th colspan="2">Anomaly</th>
     </tr>
     <tr>
-      <th>r</th><th>R²</th><th>NRMSE</th>
+      <th>r</th><th>R²</th><th>NRMSE</th><th>med R²/yr</th>
       <th>r</th><th>R²</th>
       <th>r</th><th>R²</th>
       <th>r</th><th>R²</th>
@@ -324,6 +329,7 @@ def fmt(val):
             f"<td>{fmt(mod['r'])}</td>"
             f"<td>{fmt(mod['r2'])}</td>"
             f"<td>{fmt(mod['nrmse'])}</td>"
+            f"<td>{fmt(s['r2_yearly_median'])}</td>"
             f"<td>{fmt(sp['r'])}</td>"
             f"<td>{fmt(sp['r2'])}</td>"
             f"<td>{fmt(s['r_time_model'])}</td>"
@@ -405,6 +411,7 @@ def generate_local_report_html(stats_list: List[dict], pdf_filename: str) -> str
         <td>{m['r']:.2f}</td>
         <td>{m['r2']:.2f}</td>
         <td>{m['nrmse']:.2f}</td>
+        <td>{s['r2_yearly_median']:.2f}</td>
         <td>{sp['r']:.2f}</td>
         <td>{sp['r2']:.2f}</td>
         <td>{s['r_time_model']:.2f}</td>
@@ -445,13 +452,13 @@ def generate_local_report_html(stats_list: List[dict], pdf_filename: str) -> str
         <th rowspan="2">Dataset</th>
         <th rowspan="2">N_regions</th>
         <th rowspan="2">N_years</th>
-        <th colspan="3">Region-Year</th>
+        <th colspan="4">Region-Year</th>
         <th colspan="2">Spatial</th>
         <th colspan="2">Temporal</th>
         <th colspan="2">Anomaly</th>
       </tr>
       <tr>
-        <th>r</th><th>R²</th><th>NRMSE</th>
+        <th>r</th><th>R²</th><th>NRMSE</th><th>med R²/yr</th>
         <th>r</th><th>R²</th>
         <th>r</th><th>R²</th>
         <th>r</th><th>R²</th>
@@ -681,6 +688,7 @@ def process_dataset(
         _as_float_array(ts[KEY_TARGET]),
         _as_float_array(ts[model]),
     )
+    r2_yearly_median = calc_median_yearly_r2(df_filtered, KEY_TARGET, model)
 
     # Spatial stats (mean over years per location)
     spatial = cast(pd.DataFrame, df_filtered.groupby(KEY_LOC)[[KEY_TARGET, model]].mean())
@@ -713,6 +721,7 @@ def process_dataset(
         "metrics_baseline": metrics_base,
         "r_time_model": r_time_model,
         "r2_time_model": r2_time_model,
+        "r2_yearly_median": r2_yearly_median,
         "r_time_base": r_time_base,
         "metrics_spatial": {"r": r_spatial_model, "r2": r2_spatial_model},
     }

tests/evaluation/test_aggregated_metrics.py

@@ -3,6 +3,7 @@
 
 from cybench.config import KEY_LOC, KEY_TARGET, KEY_YEAR
 from cybench.evaluation.aggregated_metrics import (
+    calc_median_yearly_r2,
     calc_nrmse,
     calc_r_r2,
     compute_report_metrics,
@@ -37,6 +38,21 @@ def test_get_metrics_dict_has_anomaly_columns():
     assert set(out) == {"r", "r2", "nrmse", "r_res", "r2_res"}
 
 
+def test_calc_median_yearly_r2():
+    rows = []
+    for loc in ("A", "B"):
+        rows.append({KEY_LOC: loc, KEY_YEAR: 2019, KEY_TARGET: 5.0, "model": 5.0})
+    rows.extend(
+        [
+            {KEY_LOC: "A", KEY_YEAR: 2020, KEY_TARGET: 5.0, "model": 5.0},
+            {KEY_LOC: "B", KEY_YEAR: 2020, KEY_TARGET: 7.0, "model": 6.0},
+        ]
+    )
+    df = pd.DataFrame(rows)
+    _, r2_2020 = calc_r_r2([5.0, 7.0], [5.0, 6.0])
+    assert calc_median_yearly_r2(df, KEY_TARGET, "model") == np.nanmedian([1.0, r2_2020])
+
+
 def test_compute_report_metrics_views():
     df = _make_df()
     out = compute_report_metrics(df, KEY_TARGET, "model")
@@ -46,6 +62,8 @@ def test_compute_report_metrics_views():
     assert "region_year" in out
     assert "spatial" in out
     assert "temporal" in out
+    assert "median_r2" in out["region_year"]
+    assert np.isfinite(out["region_year"]["median_r2"])
 
 
 def test_compute_report_metrics_ignores_nan_predictions():

tests/runs/test_collect_walk_forward_results.py

@@ -121,6 +121,7 @@ def test_summary_rows_to_dashboard_records(tmp_path: Path):
             "r2_spatial": -0.93,
             "r_temporal": 0.18,
             "r2_temporal": -5.10,
+            "r2_yearly_median": 0.42,
             "r_res": 0.05,
             "r2_res": -1.76,
         }
@@ -132,7 +133,7 @@ def test_summary_rows_to_dashboard_records(tmp_path: Path):
     (assets / "maize_NL_scatter.png").write_bytes(b"png")
 
     records = summary_rows_to_dashboard_records(rows, tmp_path)
-    assert len(records) == 9
+    assert len(records) == 10
     assert records[0]["view"] == "region_year"
     scatter_recs = [r for r in records if r.get("images", {}).get("scatter")]
     assert scatter_recs