From 47a137baebadaa650ff9f59460c8f9a2116f1930 Mon Sep 17 00:00:00 2001
From: Aaron Spring <draaronspring@proton.me>
Date: Thu, 11 Jun 2026 13:09:20 +0200
Subject: [PATCH 1/2] Fix multipletests compatibility with NumPy >= 2.0

Replace per-element vectorize=True approach with a proper core-dims call
that passes the entire flattened p-value array to statsmodels at once.
The old approach relied on numpy.vectorize implicitly coercing 1-element
arrays to scalars, which numpy 2.0 no longer allows.

Also fix test_multipletests_alpha to use explicit p-values with known
small values, since proper BH correction on random Pearson-r p-values
may not reject any tests at either alpha level.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 xskillscore/core/stattests.py       | 64 ++++++++++++++++-------------
 xskillscore/tests/test_stattests.py | 13 ++++--
 2 files changed, 46 insertions(+), 31 deletions(-)

diff --git a/xskillscore/core/stattests.py b/xskillscore/core/stattests.py
index ef1238c4..e0676b8f 100644
--- a/xskillscore/core/stattests.py
+++ b/xskillscore/core/stattests.py
@@ -1,14 +1,28 @@
 from __future__ import annotations
 
-import warnings
 from typing import Literal, Mapping, Optional, Tuple, Union
 
+import numpy as np
 import xarray as xr
 from statsmodels.stats.multitest import multipletests as statsmodels_multipletests
 
 from .types import XArray
 
 
+def _multipletests_numpy(pvals, **kwargs):
+    """Call statsmodels multipletests and broadcast scalar outputs to match pvals shape."""
+    reject, pvals_corrected, alphacSidak, alphacBonf = statsmodels_multipletests(
+        pvals, **kwargs
+    )
+    n = len(pvals)
+    return (
+        reject.astype(float),
+        pvals_corrected,
+        np.full(n, float(alphacSidak)),
+        np.full(n, float(alphacBonf)),
+    )
+
+
 def multipletests(
     p: XArray,
     alpha: float = 0.05,
@@ -99,6 +113,7 @@ def multipletests(
 
     Examples
     --------
+    >>> np.random.seed(42)
     >>> p = xr.DataArray(
     ...     np.random.normal(size=(3, 3)),
     ...     coords=[("x", np.arange(3)), ("y", np.arange(3))],
@@ -110,17 +125,17 @@ def multipletests(
             [ 0.        ,  1.        ,  1.        ],
             [ 0.        ,  0.        ,  1.        ]],
     <BLANKLINE>
-           [[ 0.49671415, -0.1382643 ,  0.64768854],
-            [ 1.        , -0.23415337, -0.23413696],
-            [ 1.        ,  0.76743473, -0.46947439]],
+           [[ 0.89408548, -0.31109468,  0.97153281],
+            [ 1.        , -1.05369019, -0.70241087],
+            [ 1.        ,  0.98670179, -4.22526947]],
     <BLANKLINE>
-           [[ 0.1       ,  0.1       ,  0.1       ],
-            [ 0.1       ,  0.1       ,  0.1       ],
-            [ 0.1       ,  0.1       ,  0.1       ]],
+           [[ 0.01163847,  0.01163847,  0.01163847],
+            [ 0.01163847,  0.01163847,  0.01163847],
+            [ 0.01163847,  0.01163847,  0.01163847]],
     <BLANKLINE>
-           [[ 0.1       ,  0.1       ,  0.1       ],
-            [ 0.1       ,  0.1       ,  0.1       ],
-            [ 0.1       ,  0.1       ,  0.1       ]]])
+           [[ 0.01111111,  0.01111111,  0.01111111],
+            [ 0.01111111,  0.01111111,  0.01111111],
+            [ 0.01111111,  0.01111111,  0.01111111]]])
     Coordinates:
       * result                (result) <U15 240B 'reject' ... 'alphacBonf'
       * x                     (x) int64 24B 0 1 2
@@ -160,24 +175,17 @@ def multipletests(
             f"Expected `return_results` from {allowed_return_results}, found {return_results}"
         )
 
-    # Suppress NumPy scalar conversion deprecation warning from internal numpy operations
-    with warnings.catch_warnings():
-        warnings.filterwarnings(
-            "ignore",
-            message="Conversion of an array with ndim > 0 to a scalar is deprecated",
-            category=DeprecationWarning,
-        )
-        ret = xr.apply_ufunc(
-            statsmodels_multipletests,
-            p.stack(s=p.dims),
-            input_core_dims=[[]],
-            vectorize=True,
-            output_core_dims=[[]] * 4,
-            output_dtypes=[bool, float, float, float],
-            kwargs=dict(method=method, alpha=alpha, **multipletests_kwargs),
-            dask="parallelized",
-            keep_attrs=keep_attrs,
-        )
+    ret = xr.apply_ufunc(
+        _multipletests_numpy,
+        p.stack(s=p.dims),
+        input_core_dims=[["s"]],
+        vectorize=False,
+        output_core_dims=[["s"], ["s"], ["s"], ["s"]],
+        output_dtypes=[float, float, float, float],
+        kwargs=dict(method=method, alpha=alpha, **multipletests_kwargs),
+        dask="parallelized",
+        keep_attrs=keep_attrs,
+    )
 
     ret = tuple(r.unstack("s").transpose(*p.dims, ...) for r in ret)
 
diff --git a/xskillscore/tests/test_stattests.py b/xskillscore/tests/test_stattests.py
index f306a68a..b797af98 100644
--- a/xskillscore/tests/test_stattests.py
+++ b/xskillscore/tests/test_stattests.py
@@ -1,4 +1,6 @@
+import numpy as np
 import pytest
+import xarray as xr
 from dask import is_dask_collection
 
 import xskillscore as xs
@@ -31,11 +33,16 @@ def test_multipletests_inputs(r_p, input, chunk):
     assert ret.coords["multipletests_alpha"] == alpha
 
 
-def test_multipletests_alpha(r_p):
+def test_multipletests_alpha():
     """Test that larger alpha leads to more rejected in multipletests."""
     method = "fdr_bh"
-    reject = multipletests(r_p, alpha=0.05, method=method).sel(result="reject")
-    reject_larger_alpha = multipletests(r_p, alpha=0.5, method=method).sel(result="reject")
+    # Use explicit p-values spanning a wide range so BH correction can distinguish alphas
+    p = xr.DataArray(
+        np.array([0.001, 0.003, 0.01, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0]),
+        dims=["x"],
+    )
+    reject = multipletests(p, alpha=0.05, method=method).sel(result="reject")
+    reject_larger_alpha = multipletests(p, alpha=0.5, method=method).sel(result="reject")
     # check more reject
     assert reject_larger_alpha.sum() > reject.sum()
 

From 755f8f20895bfd9af24df6569d028690eab1ae4b Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 11 Jun 2026 11:09:55 +0000
Subject: [PATCH 2/2] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 xskillscore/core/stattests.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/xskillscore/core/stattests.py b/xskillscore/core/stattests.py
index e0676b8f..5dbcc820 100644
--- a/xskillscore/core/stattests.py
+++ b/xskillscore/core/stattests.py
@@ -11,9 +11,7 @@
 
 def _multipletests_numpy(pvals, **kwargs):
     """Call statsmodels multipletests and broadcast scalar outputs to match pvals shape."""
-    reject, pvals_corrected, alphacSidak, alphacBonf = statsmodels_multipletests(
-        pvals, **kwargs
-    )
+    reject, pvals_corrected, alphacSidak, alphacBonf = statsmodels_multipletests(pvals, **kwargs)
     n = len(pvals)
     return (
         reject.astype(float),