crop differences

Author: michielkallenberg

cybench/runs/analysis/global_insights_lib.py

@@ -600,6 +600,96 @@ def compare_horizons(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]:
     return detail, summary
 
 
+def compare_crops_pairwise(
+    df: pd.DataFrame,
+    *,
+    crop_a: str,
+    crop_b: str,
+) -> tuple[pd.DataFrame, pd.DataFrame]:
+    """Compare NRMSE for two crops in countries that have both (inner join on country × model).
+
+    Returns (per_pair_detail, per_model_summary).
+    Delta = crop_b_nrmse − crop_a_nrmse (positive ⇒ crop_a has lower NRMSE).
+    """
+    if df.empty or crop_a == crop_b:
+        return pd.DataFrame(), pd.DataFrame()
+
+    key_cols = ["country", "model"]
+    for col in key_cols + ["crop", "nrmse"]:
+        if col not in df.columns:
+            return pd.DataFrame(), pd.DataFrame()
+
+    a_df = df[df["crop"] == crop_a][key_cols + ["nrmse", "r2", "n_samples"]].rename(
+        columns={"nrmse": "crop_a_nrmse", "r2": "crop_a_r2", "n_samples": "crop_a_samples"}
+    )
+    b_df = df[df["crop"] == crop_b][key_cols + ["nrmse", "r2", "n_samples"]].rename(
+        columns={"nrmse": "crop_b_nrmse", "r2": "crop_b_r2", "n_samples": "crop_b_samples"}
+    )
+    if a_df.empty or b_df.empty:
+        return pd.DataFrame(), pd.DataFrame()
+
+    merged = a_df.merge(b_df, on=key_cols, how="inner")
+    merged = merged[merged["crop_a_nrmse"].notna() & merged["crop_b_nrmse"].notna()].copy()
+    if merged.empty:
+        return pd.DataFrame(), pd.DataFrame()
+
+    merged["crop_a"] = crop_a
+    merged["crop_b"] = crop_b
+    merged["delta_nrmse"] = merged["crop_b_nrmse"] - merged["crop_a_nrmse"]
+    merged["delta_r2"] = merged["crop_b_r2"] - merged["crop_a_r2"]
+    merged["crop_a_better"] = merged["delta_nrmse"] > 0
+    merged["dataset_a"] = crop_a + "_" + merged["country"]
+    merged["dataset_b"] = crop_b + "_" + merged["country"]
+
+    pair_weights = merged[["crop_a_samples", "crop_b_samples"]].min(axis=1).fillna(1).clip(lower=1)
+    merged["pair_weight"] = pair_weights
+
+    model_rows: list[dict[str, Any]] = []
+    for model, grp in merged.groupby("model", sort=True):
+        weights = grp["pair_weight"]
+        model_rows.append(
+            {
+                "model": model,
+                "n_pairs": int(len(grp)),
+                "n_countries": int(grp["country"].nunique()),
+                "crop_a_win_rate": float(grp["crop_a_better"].mean()),
+                "mean_delta_nrmse": float(grp["delta_nrmse"].mean()),
+                "weighted_delta_nrmse": _weighted_mean(grp["delta_nrmse"], weights),
+                "mean_crop_a_nrmse": float(grp["crop_a_nrmse"].mean()),
+                "mean_crop_b_nrmse": float(grp["crop_b_nrmse"].mean()),
+            }
+        )
+    summary = pd.DataFrame(model_rows).sort_values(
+        ["weighted_delta_nrmse", "mean_delta_nrmse"], ascending=[False, False]
+    )
+    detail = merged.sort_values(["model", "country"]).reset_index(drop=True)
+    return detail, summary
+
+
+def build_crop_comparison_payload(df: pd.DataFrame) -> dict[str, dict[str, Any]]:
+    """Per horizon, pairwise crop NRMSE on shared countries only."""
+    if df.empty:
+        return {}
+    crops = _crop_keys(df)
+    out: dict[str, dict[str, Any]] = {}
+    for hz in ("eos", "mid"):
+        hz_df = df[df["batch_horizon"] == hz]
+        pairs: dict[str, Any] = {}
+        for i, crop_a in enumerate(crops):
+            for crop_b in crops[i + 1 :]:
+                detail, summary = compare_crops_pairwise(hz_df, crop_a=crop_a, crop_b=crop_b)
+                pair_key = f"{crop_a}_vs_{crop_b}"
+                pairs[pair_key] = {
+                    "crop_a": crop_a,
+                    "crop_b": crop_b,
+                    "detail": _df_records(detail),
+                    "summary": _df_records(summary),
+                    "overall": _overall_crop_pair_stats(detail, crop_a, crop_b),
+                }
+        out[hz] = pairs
+    return out
+
+
 def build_insights_payload(output_root: Path, *, version: int = 1) -> dict[str, Any]:
     """Build JSON-serializable payload for the global insights dashboard."""
     paths = discover_summary_tables(output_root, version=version)
@@ -640,6 +730,7 @@ def build_insights_payload(output_root: Path, *, version: int = 1) -> dict[str,
         "horizon_summary": _df_records(horizon_summary),
         "horizon_detail": _df_records(horizon_detail),
         "overall_horizon": _overall_horizon_stats(horizon_detail),
+        "crop_comparison": build_crop_comparison_payload(df),
     }
 
 
@@ -659,6 +750,41 @@ def _overall_horizon_stats(detail: pd.DataFrame) -> dict[str, Any]:
     }
 
 
+def _overall_crop_pair_stats(
+    detail: pd.DataFrame, crop_a: str, crop_b: str
+) -> dict[str, Any]:
+    if detail.empty:
+        return {
+            "crop_a": crop_a,
+            "crop_b": crop_b,
+            "interpretation": (
+                f"No paired {crop_a}/{crop_b} comparisons in countries with both crops."
+            ),
+        }
+    weights = detail["pair_weight"]
+    paired_countries = sorted(detail["country"].unique())
+    crop_a_label = crop_a.replace("_", " ")
+    crop_b_label = crop_b.replace("_", " ")
+    return {
+        "crop_a": crop_a,
+        "crop_b": crop_b,
+        "n_pairs": int(len(detail)),
+        "n_countries": int(len(paired_countries)),
+        "paired_countries": paired_countries,
+        "crop_a_win_rate": float(detail["crop_a_better"].mean()),
+        "mean_delta_nrmse": float(detail["delta_nrmse"].mean()),
+        "weighted_delta_nrmse": float(_weighted_mean(detail["delta_nrmse"], weights)),
+        "mean_crop_a_nrmse": float(detail["crop_a_nrmse"].mean()),
+        "mean_crop_b_nrmse": float(detail["crop_b_nrmse"].mean()),
+        "interpretation": (
+            f"Paired comparison in countries with both crops. "
+            f"delta_nrmse = {crop_b_label} − {crop_a_label}; "
+            f"positive values mean {crop_a_label} has lower NRMSE. "
+            f"Crop A win % = share of country×model pairs where {crop_a_label} wins."
+        ),
+    }
+
+
 def _df_records(frame: pd.DataFrame) -> list[dict[str, Any]]:
     if frame.empty:
         return []

cybench/runs/viz/global_insights_template.html

@@ -246,6 +246,38 @@ <h3 style="font-size:1rem;margin:1rem 0 0.5rem;">Detail per crop × country × m
       </p>
     </div>
 
+    <div class="card">
+      <h2>Research Q3 — Which crop is easier to predict (fair comparison)?</h2>
+      <p class="muted">
+        Pooling all countries biases crop rankings (e.g. maize-heavy regions may be harder).
+        Here we keep only <strong>countries with both crops</strong> and compare NRMSE within the same
+        country and model — a paired crop×country comparison.
+      </p>
+      <div class="controls">
+        <span class="muted">Horizon</span>
+        <div class="toolbar" id="rq3-horizon-toolbar">
+          <button type="button" data-horizon="eos" class="active">End of season</button>
+          <button type="button" data-horizon="mid">Mid-season</button>
+        </div>
+        <label>Crop pair
+          <select id="rq3-pair-filter"></select>
+        </label>
+      </div>
+      <p class="muted" id="crop-pair-note"></p>
+      <div class="stat-grid" id="crop-pair-stats"></div>
+      <h3 style="font-size:1rem;margin:1rem 0 0.5rem;">By model (paired countries only)</h3>
+      <div class="table-scroll" id="crop-pair-summary-wrap"></div>
+      <h3 style="font-size:1rem;margin:1rem 0 0.5rem;">Detail per country × model</h3>
+      <div class="controls">
+        <label>Model <input type="search" id="crop-detail-model-filter" placeholder="Filter…" /></label>
+        <label>Country
+          <select id="crop-detail-country-filter"><option value="">All</option></select>
+        </label>
+      </div>
+      <div class="table-scroll" id="crop-pair-detail-wrap"></div>
+      <p class="muted note" id="crop-pair-delta-note"></p>
+    </div>
+
     <p class="muted">
       <a href="model_families.html">Model families (paper summary)</a> ·
       <a href="index.html">← Back to country dashboards</a>
@@ -774,6 +806,150 @@ <h3 style="font-size:1rem;margin:1rem 0 0.5rem;">Detail per crop × country × m
     countryFilter.addEventListener("change", renderHorizonDetail);
     modelFilter.addEventListener("input", renderHorizonDetail);
     renderHorizonDetail();
+
+    // Crop pair comparison (RQ3)
+    let rq3Horizon = "eos";
+    let rq3PairKey = "";
+    const rq3HorizonToolbar = document.getElementById("rq3-horizon-toolbar");
+    const rq3PairFilter = document.getElementById("rq3-pair-filter");
+    const cropPairNote = document.getElementById("crop-pair-note");
+    const cropPairStats = document.getElementById("crop-pair-stats");
+    const cropPairDeltaNote = document.getElementById("crop-pair-delta-note");
+    const cropDetailCountryFilter = document.getElementById("crop-detail-country-filter");
+    const cropDetailModelFilter = document.getElementById("crop-detail-model-filter");
+    const cpsWrap = document.getElementById("crop-pair-summary-wrap");
+    const cpdWrap = document.getElementById("crop-pair-detail-wrap");
+
+    function currentCropPairSlice() {
+      const hz = (DATA.crop_comparison || {})[rq3Horizon] || {};
+      const keys = Object.keys(hz).sort();
+      if (!rq3PairKey && keys.length) rq3PairKey = keys[0];
+      if (rq3PairKey && !hz[rq3PairKey] && keys.length) rq3PairKey = keys[0];
+      return hz[rq3PairKey] || null;
+    }
+
+    function populateCropPairOptions() {
+      const hz = (DATA.crop_comparison || {})[rq3Horizon] || {};
+      const keys = Object.keys(hz).sort();
+      rq3PairFilter.innerHTML = keys.length
+        ? keys.map(k => {
+            const p = hz[k];
+            const label = `${p.crop_a} vs ${p.crop_b}`;
+            return `<option value="${k}">${label}</option>`;
+          }).join("")
+        : `<option value="">No crop pairs</option>`;
+      if (keys.length && (!rq3PairKey || !hz[rq3PairKey])) rq3PairKey = keys[0];
+      rq3PairFilter.value = rq3PairKey || "";
+    }
+
+    let cpsSort = "weighted_delta_nrmse";
+    let cpsDir = "desc";
+    let cpdSort = "delta_nrmse";
+    let cpdDir = "desc";
+
+    function renderCropPairComparison() {
+      const slice = currentCropPairSlice();
+      const overall = slice ? (slice.overall || {}) : {};
+      const cropA = overall.crop_a || "crop A";
+      const cropB = overall.crop_b || "crop B";
+      const cropALabel = cropA.charAt(0).toUpperCase() + cropA.slice(1);
+      const cropBLabel = cropB.charAt(0).toUpperCase() + cropB.slice(1);
+
+      cropPairNote.textContent = overall.interpretation || "No paired crop comparison for this horizon.";
+      cropPairDeltaNote.textContent = slice
+        ? `Δ NRMSE = ${cropBLabel} − ${cropALabel}. Positive ⇒ ${cropALabel} has lower NRMSE. `
+          + `${cropALabel} win % = share of country×model pairs where ${cropALabel} wins.`
+        : "";
+
+      if (!overall.n_pairs) {
+        cropPairStats.innerHTML = `<p class="muted">No paired comparisons available.</p>`;
+        cpsWrap.innerHTML = "";
+        cpdWrap.innerHTML = "";
+        return;
+      }
+
+      const countries = overall.paired_countries || [];
+      cropPairStats.innerHTML = `
+        <div class="stat"><div class="label">Paired countries</div><div class="value">${overall.n_countries || countries.length}</div></div>
+        <div class="stat"><div class="label">Country×model pairs</div><div class="value">${overall.n_pairs}</div></div>
+        <div class="stat"><div class="label">${cropALabel} wins</div><div class="value good">${pct(overall.crop_a_win_rate)}</div></div>
+        <div class="stat"><div class="label">Mean ${cropALabel} NRMSE</div><div class="value">${fmt(overall.mean_crop_a_nrmse)}</div></div>
+        <div class="stat"><div class="label">Mean ${cropBLabel} NRMSE</div><div class="value">${fmt(overall.mean_crop_b_nrmse)}</div></div>
+        <div class="stat"><div class="label">Weighted Δ NRMSE</div><div class="value">${fmt(overall.weighted_delta_nrmse)}</div></div>`;
+
+      const cpsCols = [
+        { key: "model", label: "Model" },
+        { key: "n_countries", label: "Countries" },
+        { key: "n_pairs", label: "Pairs" },
+        { key: "crop_a_win_rate", label: `${cropALabel} win %`, format: pct },
+        { key: "weighted_delta_nrmse", label: "Weighted Δ NRMSE", bg: fixedDeltaBg },
+        { key: "mean_delta_nrmse", label: "Mean Δ NRMSE", bg: fixedDeltaBg },
+        { key: "mean_crop_a_nrmse", label: `Mean ${cropALabel} NRMSE`, bg: fixedNrmseBg },
+        { key: "mean_crop_b_nrmse", label: `Mean ${cropBLabel} NRMSE`, bg: fixedNrmseBg },
+      ];
+      const summaryRows = [...(slice.summary || [])].sort((a, b) => {
+        const dir = cpsDir === "asc" ? 1 : -1;
+        const cmp = compareSortValues(a, b, cpsSort);
+        if (cmp !== 0) return dir * cmp;
+        return a.model.localeCompare(b.model);
+      });
+      renderSortableTable(cpsWrap, cpsCols, summaryRows, cpsSort, cpsDir, (key) => {
+        if (cpsSort === key) cpsDir = cpsDir === "asc" ? "desc" : "asc";
+        else { cpsSort = key; cpsDir = key === "crop_a_win_rate" || key.includes("delta") ? "desc" : "asc"; }
+        renderCropPairComparison();
+      });
+
+      const detailCountries = [...new Set((slice.detail || []).map(r => r.country))].sort();
+      cropDetailCountryFilter.innerHTML = `<option value="">All</option>` +
+        detailCountries.map(c => `<option value="${c}">${c}</option>`).join("");
+
+      const cpdCols = [
+        { key: "country", label: "Country" },
+        { key: "model", label: "Model" },
+        { key: "crop_a_nrmse", label: `${cropALabel} NRMSE`, bg: fixedNrmseBg },
+        { key: "crop_b_nrmse", label: `${cropBLabel} NRMSE`, bg: fixedNrmseBg },
+        { key: "delta_nrmse", label: "Δ NRMSE", bg: fixedDeltaBg },
+        { key: "crop_a_r2", label: `${cropALabel} R²` },
+        { key: "crop_b_r2", label: `${cropBLabel} R²` },
+        { key: "crop_a_samples", label: `${cropALabel} samples`, format: (v) => v == null ? "—" : String(v) },
+        { key: "crop_b_samples", label: `${cropBLabel} samples`, format: (v) => v == null ? "—" : String(v) },
+      ];
+      const cq = cropDetailCountryFilter.value;
+      const mq = cropDetailModelFilter.value.trim().toLowerCase();
+      const detailRows = (slice.detail || []).filter(r => {
+        if (cq && r.country !== cq) return false;
+        if (mq && !r.model.toLowerCase().includes(mq)) return false;
+        return true;
+      }).sort((a, b) => {
+        const dir = cpdDir === "asc" ? 1 : -1;
+        const cmp = compareSortValues(a, b, cpdSort);
+        if (cmp !== 0) return dir * cmp;
+        return a.model.localeCompare(b.model);
+      });
+      renderSortableTable(cpdWrap, cpdCols, detailRows, cpdSort, cpdDir, (key) => {
+        if (cpdSort === key) cpdDir = cpdDir === "asc" ? "desc" : "asc";
+        else { cpdSort = key; cpdDir = key === "country" || key === "model" ? "asc" : "desc"; }
+        renderCropPairComparison();
+      });
+    }
+
+    rq3HorizonToolbar.querySelectorAll("button").forEach(btn => {
+      btn.addEventListener("click", () => {
+        rq3HorizonToolbar.querySelectorAll("button").forEach(b => b.classList.remove("active"));
+        btn.classList.add("active");
+        rq3Horizon = btn.dataset.horizon;
+        populateCropPairOptions();
+        renderCropPairComparison();
+      });
+    });
+    rq3PairFilter.addEventListener("change", () => {
+      rq3PairKey = rq3PairFilter.value;
+      renderCropPairComparison();
+    });
+    cropDetailCountryFilter.addEventListener("change", renderCropPairComparison);
+    cropDetailModelFilter.addEventListener("input", renderCropPairComparison);
+    populateCropPairOptions();
+    renderCropPairComparison();
   </script>
 </body>
 </html>