diff --git a/ecosci/data/spatial.py b/ecosci/data/spatial.py
index 6716e30..c5ccbca 100644
--- a/ecosci/data/spatial.py
+++ b/ecosci/data/spatial.py
@@ -45,4 +45,22 @@ def assign_spatial_blocks(
     km = KMeans(n_clusters=n_blocks, random_state=random_state, n_init=10)
     labels = km.fit_predict(coords)
 
+    counts = pd.Series(labels).value_counts()
+    min_count = int(counts.min())
+    if min_count < 2:
+        n_singleton = int((counts < 2).sum())
+        print(
+            f"Warning: {n_singleton} of {n_blocks} spatial block(s) have fewer than "
+            f"2 samples (smallest: {min_count}). A CV fold built from one of these "
+            f"has an undefined test-set R^2/variance and will show up as NaN in "
+            f"per-fold metrics and permutation-based feature importance. Consider "
+            f"lowering n_blocks (currently {n_blocks})."
+        )
+    elif min_count < 5:
+        print(
+            f"Warning: smallest spatial block has only {min_count} samples out of "
+            f"{n_blocks} blocks. Per-fold metrics for that block may be noisy/unstable. "
+            f"Consider lowering n_blocks if this is unintentional."
+        )
+
     return pd.Series(labels, index=df.index, name="spatial_block")
diff --git a/ecosci/evaluation/feature_importance.py b/ecosci/evaluation/feature_importance.py
index 75da92e..d12a43b 100644
--- a/ecosci/evaluation/feature_importance.py
+++ b/ecosci/evaluation/feature_importance.py
@@ -119,9 +119,12 @@ def extract_feature_importance(
             
             # Average importances across all models
             if all_importances:
-                # Average across all seeds
-                avg_feature_importances = np.mean(all_importances, axis=0)
-                std_feature_importances = np.std(all_importances, axis=0)
+                # Average across all seeds (nan-aware: for regression, a
+                # single-sample test set makes permutation importance's R^2
+                # scorer undefined for that seed, and a plain mean would let
+                # that one NaN wipe out every feature's average)
+                avg_feature_importances = np.nanmean(all_importances, axis=0)
+                std_feature_importances = np.nanstd(all_importances, axis=0)
                 
                 print(f"\n{model_name.upper()}:")
                 print(f"{'-'*80}")
@@ -272,11 +275,15 @@ def extract_cv_feature_importance(
                     fold_avg_importance = np.mean(fold_importances_list, axis=0)
                     fold_importances_dict[fold_id] = fold_avg_importance
             
-            # Average importances across all folds
+            # Average importances across all folds (nan-aware: for regression,
+            # a fold whose test set has fewer than 2 samples makes permutation
+            # importance's R^2 scorer undefined for that fold, and a plain
+            # mean would let that one NaN-only fold wipe out every feature's
+            # average)
             if fold_importances_dict:
                 all_fold_importances = list(fold_importances_dict.values())
-                avg_feature_importances = np.mean(all_fold_importances, axis=0)
-                std_feature_importances = np.std(all_fold_importances, axis=0)
+                avg_feature_importances = np.nanmean(all_fold_importances, axis=0)
+                std_feature_importances = np.nanstd(all_fold_importances, axis=0)
                 
                 print(f"\n{model_name.upper()}:")
                 print(f"{'-'*80}")
diff --git a/ecosci/evaluation/metrics.py b/ecosci/evaluation/metrics.py
index 41c1c54..d0ae381 100755
--- a/ecosci/evaluation/metrics.py
+++ b/ecosci/evaluation/metrics.py
@@ -47,7 +47,13 @@ def compute_classification_metrics(y_true, y_pred, y_proba=None) -> Dict[str, An
 
     out = {}
     observed_labels = np.unique(np.concatenate([np.asarray(y_true), np.asarray(y_pred)]))
-    average_mode = "binary" if len(observed_labels) == 2 else "macro"
+    # "binary" average is only valid when the task truly has exactly 2 classes.
+    # Use y_proba's column count when available to detect multiclass tasks even if
+    # a CV fold happens to contain only 2 labels (e.g. a spatial block missing a class).
+    if y_proba is not None and getattr(y_proba, "ndim", 0) == 2:
+        average_mode = "binary" if y_proba.shape[1] == 2 else "macro"
+    else:
+        average_mode = "binary" if set(observed_labels.tolist()) == {0, 1} else "macro"
     out["accuracy"] = accuracy_score(y_true, y_pred)
     out["balanced_accuracy"] = balanced_accuracy_score(y_true, y_pred)
     out["precision"] = precision_score(