diff --git a/xskillscore/core/stattests.py b/xskillscore/core/stattests.py
index ef1238c4..5dbcc820 100644
--- a/xskillscore/core/stattests.py
+++ b/xskillscore/core/stattests.py
@@ -1,14 +1,26 @@
 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 +111,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 +123,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 +173,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()