Merge pull request #125 from maximtrp/fix/pandas-nullable-numeric-boundaries

Fix nullable numeric inputs on older pandas

Author: maximtrp

.github/workflows/package-pull.yml

@@ -30,7 +30,7 @@ jobs:
       run: |
         pytest .
 
-  pandas-categorical-compatibility:
+  pandas-compatibility:
     runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v4
@@ -43,10 +43,6 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         python -m pip install ".[test]" "numpy<2" "pandas==1.5.3"
-    - name: Testing categorical and duplicate-index inputs with pandas 1.5
+    - name: Testing full suite with pandas 1.5
       run: |
-        pytest \
-          tests/test_known_input_regressions.py::test_tukey_categorical_group_order_matches_scipy \
-          tests/test_known_input_regressions.py::test_tukey_duplicate_index_matches_scipy \
-          tests/test_posthocs.py::TestPosthocs::test_posthoc_tukey_hsd_repeated_index \
-          tests/test_independent_input_paths.py::test_unused_categories_are_not_treatments_and_ordered_categories_align_values
+        pytest .

scikit_posthocs/_posthocs.py

@@ -8,7 +8,17 @@
 import scipy.stats as ss
 from scipy import integrate, special
 from statsmodels.stats.multitest import multipletests
-from pandas import DataFrame, Series, MultiIndex
+from pandas import DataFrame, Series, MultiIndex, concat
+
+
+def __numeric_array(values: Series) -> np.ndarray:
+    """Expose nullable numeric data without object arrays or integer rounding."""
+    dtype = getattr(values.dtype, "numpy_dtype", None)
+    if dtype is None:
+        return values.to_numpy()
+    if dtype.kind in "biu" and values.hasnans:
+        dtype = float
+    return values.to_numpy(dtype=dtype, na_value=np.nan)
 
 
 def __validate_n_permutations(value: int, name: str) -> int:
@@ -153,6 +163,14 @@ def __convert_to_block_df(
         x = x.reset_index().melt(
             id_vars=new_block_id_col, var_name=new_group_col, value_name=new_y_col
         )
+        # pandas 1.5 melt loses nullable measurement dtypes. Concatenating the
+        # original columns preserves their common dtype and missing values.
+        if x[new_y_col].dtype == object and any(
+            getattr(dtype, "numpy_dtype", None) is not None for dtype in a.dtypes
+        ):
+            x[new_y_col] = concat(
+                [a.iloc[:, column] for column in range(a.shape[1])], ignore_index=True
+            )
         x[new_block_col] = a.index.take(x[new_block_id_col]).to_numpy()
         x[new_group_col] = a.columns.take(x[new_group_col]).to_numpy()
 
@@ -1015,7 +1033,7 @@ def posthoc_npm_test(
     k = groups.size
     n = x.shape[0]
     sigma = np.sqrt(n * (n + 1) / 12.0)
-    values = x[_val_col].to_numpy()
+    values = __numeric_array(x[_val_col])
     group_values = x[_group_col].to_numpy()
     group_codes = np.empty(n, dtype=int)
     for code, group in enumerate(groups):
@@ -1558,7 +1576,7 @@ def posthoc_anderson(
     groups = x[_group_col].unique()
     k = groups.size
     x_grouped = x.groupby(_group_col, observed=True)[_val_col]
-    samples = [x_grouped.get_group(group).dropna().to_numpy() for group in groups]
+    samples = [__numeric_array(x_grouped.get_group(group).dropna()) for group in groups]
     vs = np.zeros((k, k), dtype=float)
     combs = it.combinations(range(k), 2)
 
@@ -2061,7 +2079,7 @@ def posthoc_ttest(
     0.047356
     """
     x, _val_col, _group_col = __convert_to_df(a, val_col, group_col)
-    if not np.isfinite(x[_val_col].to_numpy(dtype=float)).all():
+    if not np.isfinite(x[_val_col].to_numpy(dtype=float, na_value=np.nan)).all():
         raise ValueError(
             "All observations must be finite; remove or impute missing values before calling "
             "posthoc_ttest"
@@ -2071,7 +2089,7 @@ def posthoc_ttest(
     groups = x[_group_col].unique()
     k = groups.size
     xg = x.groupby(by=_group_col, observed=True)[_val_col]
-    samples = [xg.get_group(group).to_numpy() for group in groups]
+    samples = [__numeric_array(xg.get_group(group)) for group in groups]
 
     vs = np.zeros((k, k), dtype=float)
     tri_upper = np.triu_indices(vs.shape[0], 1)
@@ -2158,7 +2176,7 @@ def posthoc_tukey_hsd(
     # Older pandas iterates ordered categoricals in category order even with
     # sort=False. Extract samples in the same order as the output labels.
     xg = x.groupby(_group_col, observed=True, sort=False)[_val_col]
-    results = ss.tukey_hsd(*[xg.get_group(group).to_numpy() for group in groups])
+    results = ss.tukey_hsd(*[__numeric_array(xg.get_group(group)) for group in groups])
 
     return DataFrame(results.pvalue, index=groups, columns=groups)
 
@@ -2240,7 +2258,7 @@ def posthoc_mannwhitney(
     x_len = groups.size
     vs = np.zeros((x_len, x_len))
     xg = x.groupby(_group_col, observed=True)[_val_col]
-    samples = [xg.get_group(group).to_numpy() for group in groups]
+    samples = [__numeric_array(xg.get_group(group)) for group in groups]
     tri_upper = np.triu_indices(vs.shape[0], 1)
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
@@ -2363,7 +2381,7 @@ def posthoc_wilcoxon(
     x_len = groups.size
     vs = np.zeros((x_len, x_len))
     xg = x.groupby(_group_col, observed=True)[_val_col]
-    samples = [xg.get_group(group).to_numpy() for group in groups]
+    samples = [__numeric_array(xg.get_group(group)) for group in groups]
     tri_upper = np.triu_indices(vs.shape[0], 1)
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
@@ -2461,8 +2479,8 @@ def posthoc_scheffe(
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
 
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
     dif = means[tri_upper[0]] - means[tri_upper[1]]
     A = sin * (1.0 / counts[tri_upper[0]] + 1.0 / counts[tri_upper[1]]) * (groups.size - 1.0)
     vs[tri_upper] = dif**2.0 / A
@@ -2549,9 +2567,9 @@ def posthoc_tamhane(
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
 
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
-    variances = si.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
+    variances = __numeric_array(si.reindex(groups))
     i, j = tri_upper
     scaled_variances = variances / counts
     A = scaled_variances[i] + scaled_variances[j]
@@ -2667,8 +2685,8 @@ def posthoc_tukey(
     vs = np.zeros((groups.size, groups.size), dtype=float)
     tri_upper = np.triu_indices(vs.shape[0], 1)
     tri_lower = np.tril_indices(vs.shape[0], -1)
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
     scale = sin * 0.5 * (1.0 / counts[tri_upper[0]] + 1.0 / counts[tri_upper[1]])
     q_values = (means[tri_upper[0]] - means[tri_upper[1]]) / np.sqrt(scale)
     vs[tri_upper] = ss.studentized_range.sf(np.abs(q_values), groups.size, n - groups.size)
@@ -2743,7 +2761,7 @@ def posthoc_dscf(
     x = x.sort_values(by=[_group_col], ascending=True) if sort else x
     groups = x[_group_col].unique()
     x_grouped = x.groupby(_group_col, observed=True)[_val_col]
-    samples = [x_grouped.get_group(group).dropna().to_numpy() for group in groups]
+    samples = [__numeric_array(x_grouped.get_group(group).dropna()) for group in groups]
     k = groups.size
 
     def compare(i, j):
@@ -2856,9 +2874,9 @@ def posthoc_games_howell(
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
 
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
-    variances = si.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
+    variances = __numeric_array(si.reindex(groups))
     i, j = tri_upper
     A = variances[i] / counts[i] + variances[j] / counts[j]
     q_values = (means[i] - means[j]) / np.sqrt(A) * np.sqrt(2.0)
@@ -2982,9 +3000,9 @@ def posthoc_dunnett_t3(
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
 
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
-    variances = si.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
+    variances = __numeric_array(si.reindex(groups))
     i, j = tri_upper
     A = variances[i] / counts[i] + variances[j] / counts[j]
     t_values = (means[i] - means[j]) / np.sqrt(A)
@@ -3076,8 +3094,8 @@ def posthoc_lsd(
     tri_lower = np.tril_indices(vs.shape[0], -1)
     vs[:, :] = 0
 
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
     i, j = tri_upper
     A = s2in * (1.0 / counts[i] + 1.0 / counts[j])
     t_values = (means[i] - means[j]) / np.sqrt(A)
@@ -3166,8 +3184,8 @@ def posthoc_snk(
     df = n - k
 
     # rank position of each group when means are sorted in decreasing order
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
+    counts = __numeric_array(ni.reindex(groups))
+    means = __numeric_array(xi.reindex(groups))
     order = means.argsort()[::-1]
     rank_pos = np.empty(k, dtype=int)
     rank_pos[order] = np.arange(k)
@@ -3257,8 +3275,7 @@ def posthoc_duncan(
     df = n - k
     r = k / np.sum(1.0 / ni)
 
-    counts = ni.reindex(groups).to_numpy()
-    means = xi.reindex(groups).to_numpy()
+    means = __numeric_array(xi.reindex(groups))
     order = means.argsort()[::-1]
     rank_pos = np.empty(k, dtype=int)
     rank_pos[order] = np.arange(k)

tests/test_input_conversion_contracts.py

@@ -70,19 +70,20 @@ def test_optimized_block_gate_rejects_inputs_requiring_conversion(data, melted,
     np.testing.assert_equal(result is None, True)
 
 
-def test_block_converter_uses_row_ids_without_losing_measurement_dtype():
+@pytest.mark.parametrize("dtype", ["float32", "float64", "Float32", "Float64", "Int64", "UInt64"])
+def test_block_converter_uses_row_ids_without_losing_measurement_dtype(dtype):
     data = pd.DataFrame(
         {
-            "y": pd.array([1, None, 3], dtype="Float64"),
-            "block_ids": pd.array([4, 5, 6], dtype="Float64"),
+            "y": pd.array([1, None, 3], dtype=dtype),
+            "block_ids": pd.array([4, 5, 6], dtype=dtype),
         }
     )
     data.index = ["subject", "subject", "other"]
     original = data.copy(deep=True)
 
     actual, value_col, group_col, block_col, block_id_col = sp.__convert_to_block_df(data)
 
-    np.testing.assert_equal(actual[value_col].dtype, pd.Float64Dtype())
+    np.testing.assert_equal(actual[value_col].dtype, pd.api.types.pandas_dtype(dtype))
     np.testing.assert_equal(len(actual), data.size)
     np.testing.assert_array_equal(actual[group_col], ["y"] * 3 + ["block_ids"] * 3)
     np.testing.assert_array_equal(actual[block_id_col], [0, 1, 2, 0, 1, 2])

tests/test_nullable_numeric_inputs.py

@@ -0,0 +1,88 @@
+import numpy as np
+import pandas as pd
+import pytest
+from scipy import stats
+
+import scikit_posthocs as sp
+import scikit_posthocs._posthocs as posthocs
+
+
+@pytest.mark.parametrize("dtype", ["Float32", "Float64", "Int64", "UInt64", "boolean"])
+@pytest.mark.parametrize("missing", [False, True])
+def test_nullable_numeric_array_has_native_dtype_and_preserves_missing_values(dtype, missing):
+    values = [True, False, True] if dtype == "boolean" else [1, 2, 3]
+    if missing:
+        values[1] = None
+    data = pd.Series(values, dtype=dtype)
+    original = data.copy(deep=True)
+    expected_dtype = data.dtype.numpy_dtype
+    if missing and expected_dtype.kind in "biu":
+        expected_dtype = np.dtype(float)
+    expected = np.array(
+        [np.nan if value is None else value for value in values], dtype=expected_dtype
+    )
+
+    actual = posthocs.__numeric_array(data)
+
+    np.testing.assert_equal(actual.dtype, expected_dtype)
+    np.testing.assert_array_equal(actual, expected)
+    pd.testing.assert_series_equal(data, original)
+
+
+@pytest.mark.parametrize("dtype", ["float32", "float64", "int64", "uint64", "bool"])
+def test_native_numeric_array_does_not_change_dtype(dtype):
+    data = pd.Series([0, 1, 0], dtype=dtype)
+
+    actual = posthocs.__numeric_array(data)
+
+    np.testing.assert_equal(actual.dtype, data.dtype)
+    np.testing.assert_array_equal(actual, data.to_numpy())
+
+
+@pytest.mark.parametrize("dtype", ["Int64", "UInt64"])
+@pytest.mark.parametrize("sort", [False, True])
+def test_nullable_large_integers_do_not_gain_ties_from_float_conversion(dtype, sort):
+    offset = 2**53
+    first = np.array([offset, offset + 1, offset + 4], dtype=dtype.lower())
+    second = np.array([offset + 2, offset + 3, offset + 5], dtype=dtype.lower())
+    data = pd.DataFrame(
+        {
+            "value": pd.array(np.concatenate([first, second]), dtype=dtype),
+            "group": ["z"] * 3 + ["a"] * 3,
+        }
+    )
+    original = data.copy(deep=True)
+
+    actual = sp.posthoc_mannwhitney(data, val_col="value", group_col="group", sort=sort)
+
+    np.testing.assert_allclose(actual.loc["z", "a"], stats.mannwhitneyu(first, second).pvalue)
+    native = posthocs.__numeric_array(data["value"])
+    np.testing.assert_equal(native.dtype, np.dtype(dtype.lower()))
+    np.testing.assert_array_equal(native, np.concatenate([first, second]))
+    pd.testing.assert_frame_equal(data, original)
+
+
+def test_nullable_missing_ttest_observations_keep_finite_input_error():
+    data = pd.DataFrame(
+        {"value": pd.array([1, None, 3, 5, 7, 9], dtype="Float64"), "group": ["a"] * 3 + ["b"] * 3}
+    )
+
+    with pytest.raises(ValueError, match="All observations must be finite"):
+        sp.posthoc_ttest(data, val_col="value", group_col="group")
+
+
+@pytest.mark.parametrize("dtype", ["Float64", "Int64"])
+def test_nullable_wilcoxon_matches_scipy_and_rejects_missing_pairs(dtype):
+    first, second = [1, 2, 4, 5, 8], [2, 4, 7, 8, 12]
+    data = pd.DataFrame(
+        {"value": pd.array(first + second, dtype=dtype), "group": ["a"] * 5 + ["b"] * 5}
+    )
+    original = data.copy(deep=True)
+
+    actual = sp.posthoc_wilcoxon(data, val_col="value", group_col="group")
+
+    np.testing.assert_allclose(actual.loc["a", "b"], stats.wilcoxon(first, second).pvalue)
+    pd.testing.assert_frame_equal(data, original)
+    data.loc[0, "value"] = pd.NA
+    with pytest.raises(ValueError, match="missing paired observations"):
+        sp.posthoc_wilcoxon(data, val_col="value", group_col="group")