Include missing values as a separate PSI bin

docs/user_guide/selection/DropHighPSIFeatures.rst

@@ -961,3 +961,32 @@ For more details about this and other feature selection methods check out these
 Both our book and course are suitable for beginners and more advanced data scientists
 alike. By purchasing them you are supporting `Sole <https://linkedin.com/in/soledad-galli>`_,
 the main developer of feature-engine.
+
+
+Missing values as a separate bin
+--------------------------------
+
+Use ``missing_values="include"`` to measure changes in missingness as well as
+changes in observed values::
+
+    selector = DropHighPSIFeatures(
+        missing_values="include", bins=5, threshold=0.25
+    )
+    selector.fit(X_train)
+    print(selector.psi_values_)
+
+Numerical bin boundaries are estimated from non-missing values in the basis
+partition only. Missing values occupy a separate bin, and all rows contribute
+to the frequency denominator. Categorical missing values also occupy their own
+bin without needing a special category string. Existing ``"raise"`` and
+``"ignore"`` modes are unchanged.
+
+If the basis partition has fewer than two distinct observed numerical values,
+there is insufficient information to estimate boundaries: all observed values
+are put in one bin and the comparison measures missingness drift only. This
+also handles an entirely missing partition without learning boundaries from
+the test partition. Infinite values are rejected.
+
+With ``threshold="auto"``, the number of populated bins in the combined
+partitions includes the missing bin when present. Empty-bin smoothing continues
+to use ``min_pct_empty_bins``.

feature_engine/selection/drop_psi_features.py

@@ -216,6 +216,12 @@ class DropHighPSIFeatures(BaseSelector):
         Takes values 'raise' or 'ignore'. If 'ignore', missing values will be dropped
         when determining the PSI for that particular feature. If 'raise' the transformer
         will raise an error and features will not be selected.
+        With 'include', missing values form a separate bin and frequencies use all
+        rows. Numerical boundaries are learned only from observed basis values.
+        If the basis has fewer than two distinct observed values, all observed
+        values are placed in one bin, so only missingness drift can be measured.
+        The automatic threshold uses the number of populated bins across both
+        partitions, including the missing bin when present.
 
     p_value: float, default = 0.001
         The p-value to test the null hypothesis that there is no feature drift. In that
@@ -366,9 +372,9 @@ def __init__(
                 f"instead."
             )
 
-        if missing_values not in ["raise", "ignore"]:
+        if missing_values not in ["raise", "ignore", "include"]:
             raise ValueError(
-                f"missing_values takes only values 'raise' or 'ignore'. Got "
+                f"missing_values takes only values 'raise', 'ignore' or 'include'. Got "
                 f"{missing_values} instead."
             )
 
@@ -430,6 +436,8 @@ def fit(self, X: pd.DataFrame, y: pd.Series = None):
             # check if dataset contains na or inf
             _check_contains_na(X, num_variables_ + cat_variables_)
             _check_contains_inf(X, num_variables_)
+        elif self.missing_values == "include":
+            _check_contains_inf(X, num_variables_)
 
         # Split the dataframe into basis and test.
         basis_df, test_df = self._split_dataframe(X)
@@ -482,14 +490,24 @@ def fit(self, X: pd.DataFrame, y: pd.Series = None):
         # Compute PSI for numerical features
         for feature in num_variables_:
             # Discretize feature
-            basis_discrete = bucketer.fit_transform(basis_df[[feature]].dropna())
-            test_discrete = bucketer.transform(test_df[[feature]].dropna())
+            if self.missing_values == "include":
+                basis_discrete, test_discrete = self._discretize_with_missing(
+                    basis_df[[feature]], test_df[[feature]], bucketer
+                )
+            else:
+                basis_discrete = bucketer.fit_transform(basis_df[[feature]].dropna())
+                test_discrete = bucketer.transform(test_df[[feature]].dropna())
 
             # Determine percentage of observations per bin
             basis_distrib, test_distrib = self._observation_frequency_per_bin(
                 basis_discrete, test_discrete
             )
 
+            if self.missing_values == "include" and self.threshold == "auto":
+                threshold_num = self._missing_bin_threshold(
+                    basis_df.shape[0], test_df.shape[0], len(basis_distrib)
+                )
+
             # Calculate the PSI value
             self.psi_values_[feature] = np.sum(
                 (test_distrib - basis_distrib) * np.log(test_distrib / basis_distrib)
@@ -516,12 +534,15 @@ def fit(self, X: pd.DataFrame, y: pd.Series = None):
 
             # Determine the appropriate threshold for the categorical feature
             if self.threshold == "auto":
-                n_bins_cat = X[feature].nunique()
-                threshold_cat = self._calculate_auto_threshold(
-                    basis_df.shape[0],
-                    test_df.shape[0],
-                    n_bins_cat,
-                )
+                if self.missing_values == "include":
+                    threshold_cat = self._missing_bin_threshold(
+                        basis_df.shape[0], test_df.shape[0], len(basis_distrib)
+                    )
+                else:
+                    n_bins_cat = X[feature].nunique()
+                    threshold_cat = self._calculate_auto_threshold(
+                        basis_df.shape[0], test_df.shape[0], n_bins_cat
+                    )
 
             # Assess if feature should be dropped
             if self.psi_values_[feature] > threshold_cat:
@@ -613,8 +634,9 @@ def _observation_frequency_per_bin(self, basis, test):
             Test Pandas Series with percentage of observations per bin.
         """
         # Compute the feature distribution for basis and test
-        basis_distrib = basis.value_counts(normalize=True)
-        test_distrib = test.value_counts(normalize=True)
+        dropna = self.missing_values != "include"
+        basis_distrib = basis.value_counts(normalize=True, dropna=dropna)
+        test_distrib = test.value_counts(normalize=True, dropna=dropna)
 
         # Align the two distributions by merging the buckets (bins). This ensures
         # the number of bins is the same for the two distributions (in case of
@@ -634,6 +656,31 @@ def _observation_frequency_per_bin(self, basis, test):
 
         return distributions.basis, distributions.test
 
+    def _discretize_with_missing(self, basis, test, bucketer):
+        """Learn boundaries from observed basis values; reserve -1 for missingness."""
+        observed = basis.dropna()
+        basis_codes = pd.DataFrame(-1.0, index=basis.index, columns=basis.columns)
+        test_codes = pd.DataFrame(-1.0, index=test.index, columns=test.columns)
+        if observed.iloc[:, 0].nunique() <= 1:
+            # No distribution boundaries can be estimated. Compare missingness
+            # against a single observed-value bin, without learning from test data.
+            return basis_codes.mask(basis.notna(), 0), test_codes.mask(test.notna(), 0)
+        basis_codes.iloc[np.flatnonzero(basis.iloc[:, 0].notna())] = (
+            bucketer.fit_transform(observed).to_numpy()
+        )
+        observed_test = test.dropna()
+        if not observed_test.empty:
+            test_codes.iloc[np.flatnonzero(test.iloc[:, 0].notna())] = (
+                bucketer.transform(observed_test).to_numpy()
+            )
+        return basis_codes, test_codes
+
+    def _missing_bin_threshold(self, n_basis, n_test, n_bins):
+        # A single populated bin cannot drift; chi-square with zero df is undefined.
+        if n_bins <= 1:
+            return 0.0
+        return self._calculate_auto_threshold(n_basis, n_test, n_bins)
+
     def _split_dataframe(self, X: pd.DataFrame):
         """
         Split dataframe according to a cut-off value and return two dataframes: the

tests/estimator_checks/sklearn_check_wrapper.py

@@ -20,6 +20,7 @@
 
 import numpy as np
 import pandas as pd
+from scipy.sparse import issparse
 from sklearn.base import BaseEstimator, TransformerMixin, clone
 from sklearn.utils.validation import check_is_fitted
 
@@ -62,6 +63,9 @@ def _to_df(X):
         if hasattr(X, "iloc"):
             return X
 
+        if issparse(X):
+            raise TypeError("Sparse input is not supported by this transformer.")
+
         arr = np.asarray(X)
         if arr.ndim == 0:
             raise ValueError(

tests/test_selection/test_psi_missing_bin.py

@@ -0,0 +1,111 @@
+import numpy as np
+import pandas as pd
+import pytest
+
+from feature_engine.selection import DropHighPSIFeatures
+
+
+@pytest.mark.parametrize("strategy", ["equal_width", "equal_frequency"])
+@pytest.mark.parametrize("categorical", [False, True])
+def test_missingness_drift_has_its_own_bin(strategy, categorical):
+    basis = [1.0, 1.0, 2.0, 2.0, np.nan, np.nan, np.nan, np.nan]
+    test = [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, np.nan, np.nan]
+    X = pd.DataFrame({"x": basis + test})
+    if categorical:
+        X["x"] = X["x"].map({1.0: "a", 2.0: "b"})
+    selector = DropHighPSIFeatures(
+        variables="all", missing_values="include", bins=2, strategy=strategy
+    ).fit(X)
+    expected = 0.25 * np.log(0.5 / 0.25) + 2 * 0.125 * np.log(0.375 / 0.25)
+    assert selector.psi_values_["x"] == pytest.approx(expected)
+    ignored = DropHighPSIFeatures(
+        variables="all", missing_values="ignore", bins=2, strategy=strategy
+    ).fit(X)
+    assert ignored.psi_values_["x"] == pytest.approx(0.0)
+
+
+@pytest.mark.parametrize(
+    "basis,test",
+    [
+        ([np.nan] * 4, [np.nan] * 4),
+        ([np.nan] * 4, [1.0, 2.0, 3.0, 4.0]),
+        ([1.0, 2.0, 3.0, 4.0], [np.nan] * 4),
+    ],
+)
+@pytest.mark.parametrize("threshold", [0.25, "auto"])
+def test_all_missing_partition_is_supported(basis, test, threshold):
+    selector = DropHighPSIFeatures(
+        missing_values="include", bins=2, threshold=threshold
+    ).fit(pd.DataFrame({"x": basis + test}))
+    assert np.isfinite(selector.psi_values_["x"])
+    if all(pd.isna(basis)) and all(pd.isna(test)):
+        assert selector.psi_values_["x"] == 0
+        assert selector.features_to_drop_ == []
+    else:
+        assert selector.features_to_drop_ == ["x"]
+
+
+def test_auto_threshold_counts_missing_category():
+    X = pd.DataFrame({"x": ["a", None, "b", None, "a", "a", "b", None]})
+    selector = DropHighPSIFeatures(
+        missing_values="include", variables="all", bins=2, threshold="auto"
+    ).fit(X)
+    expected_threshold = selector._calculate_auto_threshold(4, 4, 3)
+    assert ("x" in selector.features_to_drop_) == (
+        selector.psi_values_["x"] > expected_threshold
+    )
+
+
+def test_missing_bin_does_not_collide_with_real_category():
+    X = pd.DataFrame({"x": ["missing", None] * 4})
+    selector = DropHighPSIFeatures(missing_values="include", variables="all", bins=2)
+    left, right = selector._observation_frequency_per_bin(X.iloc[:4], X.iloc[4:])
+    assert len(left) == len(right) == 2
+    assert left.sum() == right.sum() == 1
+
+
+def test_duplicate_row_labels_do_not_duplicate_observations():
+    X = pd.DataFrame(
+        {
+            "period": [0] * 8 + [1] * 8,
+            "x": [
+                1.0,
+                1.0,
+                2.0,
+                2.0,
+                np.nan,
+                np.nan,
+                np.nan,
+                np.nan,
+                1.0,
+                1.0,
+                1.0,
+                2.0,
+                2.0,
+                2.0,
+                np.nan,
+                np.nan,
+            ],
+        },
+        index=[0] * 16,
+    )
+    selector = DropHighPSIFeatures(
+        split_col="period", cut_off=0, missing_values="include", bins=2
+    ).fit(X)
+    expected = 0.25 * np.log(2) + 0.25 * np.log(1.5)
+    assert selector.psi_values_["x"] == pytest.approx(expected)
+
+
+def test_infinite_values_are_not_treated_as_missing():
+    X = pd.DataFrame({"x": [1.0, 2.0, 3.0, np.inf] * 2})
+    with pytest.raises(ValueError, match="inf"):
+        DropHighPSIFeatures(missing_values="include", bins=2).fit(X)
+
+
+def test_all_missing_categorical_feature_has_zero_psi():
+    X = pd.DataFrame({"x": pd.Series([None] * 8, dtype="object")})
+    selector = DropHighPSIFeatures(
+        missing_values="include", variables="all", bins=2, threshold="auto"
+    ).fit(X)
+    assert selector.psi_values_["x"] == 0
+    assert selector.features_to_drop_ == []

tests/test_sklearn_compatible/test_check_input_wrapper.py

@@ -0,0 +1,17 @@
+import numpy as np
+import pytest
+from scipy import sparse
+
+from tests.estimator_checks.sklearn_check_wrapper import _SklearnCheckInputWrapper
+
+
+@pytest.mark.parametrize("format", ["bsr", "coo", "csc", "csr", "dia", "dok", "lil"])
+@pytest.mark.parametrize("kind", ["matrix", "array"])
+def test_sparse_input_raises_explicit_error(format, kind):
+    constructor = getattr(sparse, f"{format}_{kind}", None)
+    if constructor is None:
+        pytest.skip("Sparse array format unavailable in this SciPy version.")
+    X = constructor(np.eye(3))
+
+    with pytest.raises(TypeError, match="[Ss]parse.*not supported"):
+        _SklearnCheckInputWrapper._to_df(X)