* compute metrics on residuals

* compute median, rather than mean
* add pearson-R

Author: michielkallenberg

cybench/evaluation/eval.py

@@ -1,5 +1,6 @@
 import numpy as np
 from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error, r2_score
+from scipy.stats import pearsonr
 
 from cybench.models.model import BaseModel
 from cybench.datasets.dataset import Dataset
@@ -13,9 +14,14 @@ def metric(func):
     implemented_metrics[func.__name__] = func
     return func
 
+def get_default_metrics(residual: bool = False):
+    if residual:
+        # Only metrics that make sense on residuals
+        return ("r", "r2")
+    else:
+        # Full set
+        return ("mape", "normalized_rmse", "r", "r2")
 
-def get_default_metrics():
-    return ("normalized_rmse", "mape", "r2")
 
 
 def evaluate_model(model: BaseModel, dataset: Dataset, metrics=get_default_metrics()):
@@ -58,7 +64,7 @@ def evaluate_predictions(
         metric_function = implemented_metrics.get(metric_name)
         if metric_function:
             result = metric_function(y_true, y_pred)
-            results[metric_name] = result
+            results[metric_name] = np.round(result, 4)
         else:
             raise ValueError(f"Metric function '{metric_name}' not implemented.")
 
@@ -114,3 +120,21 @@ def r2(y_true: np.ndarray, y_pred: np.ndarray):
     """
 
     return r2_score(y_true, y_pred)
+
+
+@metric
+def r(y_true: np.ndarray, y_pred: np.ndarray):
+    """
+    Calculate Pearson correlation coefficient.
+
+    Args:
+    - y_true (numpy.ndarray): True values.
+    - y_pred (numpy.ndarray): Predicted values.
+
+    Returns:
+    - float: Pearson's R.
+    """
+    try:
+        return pearsonr(y_true, y_pred)[0]
+    except Exception:
+        return float("nan")

cybench/runs/process_results.py

@@ -103,9 +103,9 @@ def write_results_to_table(output_file: str):
     df_metrics = results_to_metrics()
     default_metrics = get_default_metrics()
     metrics = [m for m in default_metrics if m in df_metrics.columns]
-    df_metrics = df_metrics.groupby(["crop", KEY_COUNTRY, "model"], observed=True).agg(
-        {m: "mean" for m in metrics}
-    )
+    df_metrics = df_metrics.groupby(
+        ["crop", KEY_COUNTRY, "model"], observed=True
+    ).agg({m: "median" for m in metrics})
     crops = df_metrics.index.get_level_values("crop").unique()
     metrics = df_metrics.columns.unique()
     tables = {}
@@ -116,8 +116,7 @@ def write_results_to_table(output_file: str):
             tables[crop][metric] = crop_df.reset_index().pivot_table(
                 index=["crop", KEY_COUNTRY], columns="model", values=metric
             )
-
-    # Open a file to write Markdown content
+    print(f"write to {os.path.join(PATH_OUTPUT_DIR, output_file)}")
     with open(os.path.join(PATH_OUTPUT_DIR, output_file), "w") as file:
         for crop, metrics in tables.items():
             for metric, values in metrics.items():
@@ -132,12 +131,12 @@ def write_results_to_table(output_file: str):
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(
-        prog="process_results.py", description="Output markdown tables with summary of metrics"
+        prog="process_results.py",
+        description="Output markdown tables with summary of metrics",
     )
     parser.add_argument("-o", "--output_file")
     args = parser.parse_args()
     output_file = "output_tables.md"
-    if (args.output_file is not None):
+    if args.output_file is not None:
         output_file = args.output_file
-
     write_results_to_table(output_file=output_file)

cybench/runs/run_benchmark.py

@@ -16,7 +16,7 @@
 )
 
 from cybench.datasets.dataset import Dataset
-from cybench.evaluation.eval import evaluate_predictions
+from cybench.evaluation.eval import evaluate_predictions, get_default_metrics
 from cybench.models.naive_models import AverageYieldModel
 from cybench.models.trend_models import TrendModel
 from cybench.models.sklearn_models import SklearnRidge, SklearnRandomForest
@@ -252,55 +252,76 @@ def get_prediction_residuals(run_name: str, model_names: dict) -> pd.DataFrame:
     return df_all
 
 
+def _prepare_targets_preds(df_yr, model_name, y_loc_mean=None, residual=False):
+    """Prepare y_true and y_pred, optionally using residuals."""
+    y_true = df_yr[KEY_TARGET].values
+    y_pred = df_yr[model_name].values
+
+    if residual and y_loc_mean is not None:
+        y_true = y_true - df_yr[KEY_LOC].map(y_loc_mean)
+        y_pred = y_pred - df_yr[KEY_LOC].map(y_loc_mean)
+
+    return y_true, y_pred
+
+
 def compute_metrics(
     run_name: str,
     model_names: list = None,
+    residual: bool = False,
 ) -> pd.DataFrame:
     """
     Compute evaluation metrics on saved predictions.
+
     Args:
-        run_name (str): The name of the run. Will be used to store log files and model results
-        model_names (list) : names of models
+        run_name (str): The name of the run. Will be used to store log files and model results.
+        model_names (list): Names of models to evaluate. If None, all model columns are used.
+        residual (bool): If True, compute metrics on residuals (values adjusted per location).
 
     Returns:
-        a pd.DataFrame containing evaluation metrics
+        pd.DataFrame containing evaluation metrics
     """
     df_all = load_results(run_name)
     if df_all.empty:
         return pd.DataFrame(columns=[KEY_COUNTRY, "model", KEY_YEAR])
 
     rows = []
     country_codes = df_all[KEY_COUNTRY].unique()
+
     for cn in country_codes:
         df_cn = df_all[df_all[KEY_COUNTRY] == cn]
         all_years = sorted(df_cn[KEY_YEAR].unique())
+
+        # Precompute location means for residuals
+        y_loc_mean = df_cn.groupby(KEY_LOC)[KEY_TARGET].mean() if residual else None
+
         for yr in all_years:
             df_yr = df_cn[df_cn[KEY_YEAR] == yr]
-            y_true = df_yr[KEY_TARGET].values
+
             if model_names is None:
                 model_names = [
-                    c
-                    for c in df_yr.columns
+                    c for c in df_yr.columns
                     if c not in [KEY_COUNTRY, KEY_LOC, KEY_YEAR, KEY_TARGET]
                 ]
 
             for model_name in model_names:
-                metrics = evaluate_predictions(y_true, df_yr[model_name].values)
+                y_true, y_pred = _prepare_targets_preds(df_yr, model_name, y_loc_mean, residual)
+
+                # Select metrics based on residual mode
+                metrics_to_use = get_default_metrics(residual=residual)
+                metrics = evaluate_predictions(y_true, y_pred, metrics=metrics_to_use)
+
                 metrics_row = {
                     KEY_COUNTRY: cn,
                     "model": model_name,
                     KEY_YEAR: yr,
+                    **metrics
                 }
-
-                for metric_name, value in metrics.items():
-                    metrics_row[metric_name] = value
-
                 rows.append(metrics_row)
 
-    df_all = pd.DataFrame(rows)
-    df_all.set_index([KEY_COUNTRY, "model", KEY_YEAR], inplace=True)
+    df_out = pd.DataFrame(rows)
+    df_out.set_index([KEY_COUNTRY, "model", KEY_YEAR], inplace=True)
 
-    return df_all
+    return df_out
 
 
 def run_benchmark_on_all_data():