diff --git a/src/modelskill/metrics.py b/src/modelskill/metrics.py
index 5a6f008a0..0e80c5272 100644
--- a/src/modelskill/metrics.py
+++ b/src/modelskill/metrics.py
@@ -20,7 +20,6 @@
 import numpy as np
 import pandas as pd
 from numpy.typing import ArrayLike
-from scipy import stats
 
 from .settings import options
 
@@ -590,6 +589,9 @@ def peak_ratio(
     # Use total_seconds() to handle any datetime precision (ns, us, ms, s)
     dt = time[1:] - time[:-1]
     dt_seconds = dt.total_seconds().values
+    # Lazy import: scipy.stats is heavy and only used by a few metrics, not on package import
+    from scipy import stats
+
     dt_mode_seconds = float(stats.mode(dt_seconds, keepdims=False)[0])
     N_years = dt_mode_seconds / 24 / 3600 / 365.25 * len(time)
     peak_index, AAP_ = _partial_duration_series(
diff --git a/src/modelskill/plotting/_scatter.py b/src/modelskill/plotting/_scatter.py
index e8d5c0600..013aa1f84 100644
--- a/src/modelskill/plotting/_scatter.py
+++ b/src/modelskill/plotting/_scatter.py
@@ -11,7 +11,6 @@
 from matplotlib import patches
 from matplotlib.axes import Axes
 from matplotlib.ticker import MaxNLocator
-from scipy import interpolate
 import pandas as pd
 
 import modelskill.settings as settings
@@ -766,6 +765,9 @@ def __scatter_density(x, y, binsize: float = 0.1, method: str = "linear"):
     yg = yg.ravel()
 
     ## Interpolate histogram density data to scatter data
+    # Lazy import: scipy.interpolate is heavy and only needed here, not on package import
+    from scipy import interpolate
+
     Z_grid = interpolate.griddata((xg, yg), hist, (x, y), method=method)
 
     # Replace negative values (should there be some) in case of 'cubic' interpolation