add plotting function to new cluster test api

doc/changes/dev/14399.newfeature.rst

@@ -0,0 +1 @@
+Add :meth:`mne.stats.ClusterResult.plot_joint` to visualize the observed cluster test statistic as a joint plot, by `Carina Forster`_.

mne/stats/cluster_level.py

@@ -2903,3 +2903,128 @@ def __repr__(self):  # noqa: D105
             f"<ClusterResult | p={self.cluster_p_values.min()}, "
             f"{len(self.clusters)} clusters{contrast}>"
         )
+
+    @staticmethod
+    def _get_stat_unit(stat_fun, default="stat"):
+        """Return the unit label for a test statistic (e.g. "t" or "F")."""
+        if stat_fun is f_oneway:
+            return "F"
+        if stat_fun is ttest_1samp_no_p:
+            return "t"
+        if isinstance(stat_fun, partial) and stat_fun.func is _rm_anova_stat_fun:
+            return "F"
+        return default
+
+    def plot_joint(
+        self,
+        info,
+        times_sec,
+        times="peaks",
+        title=None,
+        picks=None,
+        exclude="bads",
+        show=True,
+        ts_args=None,
+        topomap_args=None,
+    ):
+        """Plot the observed cluster test statistic as a joint plot.
+
+        Visualizes test statistics as a time series (butterfly plot) with
+        topographic maps at selected time points, preserving spatial relationships
+        between channels.
+
+        .. note::
+            Only results of shape ``(n_channels, n_times)`` (e.g. from
+            :class:`~mne.Evoked` input) are supported. Time-frequency results
+            are not supported yet.
+
+        Parameters
+        ----------
+        info : instance of Info
+            The Info object from the evoked data, containing channel information.
+        times_sec : array-like
+            Must have length stat_obs.shape[1] and spacing 1 / info["sfreq"].
+        times : float | array of float | "auto" | "peaks"
+            The time point(s) at which to show topomaps. If ``"auto"``, 5 evenly
+            spaced topographies will be shown. If ``"peaks"``, automatically finds
+            time points by checking for local maxima. Defaults to ``"peaks"``.
+        title : str | None
+            The title for the plot. If ``None``, a default title is created based
+            on the test statistic name and contrast. If an empty string, no title
+            is shown. Defaults to ``None``.
+        picks : str | array-like | slice | None
+            Channels to include. Defaults to ``None`` (all channels).
+        exclude : list of str | "bads"
+            Channels to exclude. Defaults to ``"bads"``.
+        show : bool
+            Show figure if ``True``. Defaults to ``True``.
+        ts_args : dict | None
+            Arguments forwarded to :func:`mne.Evoked.plot`.
+            Defaults to ``None``.
+        topomap_args : dict | None
+            Arguments forwarded to :func:`mne.Evoked.plot_topomap`.
+            Defaults to ``None``.
+
+        Returns
+        -------
+        fig : instance of matplotlib.figure.Figure | list
+            The figure object(s) containing the plot.
+        """
+        from mne import EvokedArray
+        from mne.viz.evoked import plot_evoked_joint
+
+        # for now only supports 2D clusters
+        if self.stat_obs.ndim != 2:
+            raise ValueError(
+                "plot_joint supports only (n_channels, n_times) statistics, got "
+                f"stat_obs with shape {self.stat_obs.shape}. Time-frequency "
+                "results are not supported yet."
+            )
+        if len(info["ch_names"]) != self.stat_obs.shape[0]:
+            raise ValueError(
+                f"info has {len(info['ch_names'])} channels, but stat_obs has "
+                f"{self.stat_obs.shape[0]}."
+            )
+
+        times_sec = np.asarray(times_sec, float)
+
+        stat_evoked = EvokedArray(self.stat_obs, info, tmin=times_sec[0])
+
+        if times_sec.shape != stat_evoked.times.shape or not np.allclose(
+            stat_evoked.times, times_sec
+        ):
+            raise ValueError(
+                "times_sec does not match the time axis of stat_obs given "
+                f"info['sfreq']={info['sfreq']}."
+            )
+
+        # Create default title if not provided
+        if title is None:
+            if self.contrast is not None:
+                title = (
+                    f"Observed cluster {self.stat_name} "
+                    f"({self.contrast[0]} − {self.contrast[1]})"
+                )
+            else:
+                title = f"Observed cluster {self.stat_name}"
+
+        # Test statistics are not in Volts: disable the V -> uV scaling
+        # and label the axes with t-values instead.
+        unit = self._get_stat_unit(self.stat_fun)
+        ch_types = set(stat_evoked.get_channel_types())
+        scalings = {ch_type: 1.0 for ch_type in ch_types}
+        units = {ch_type: unit for ch_type in ch_types}
+
+        ts_args = {"scalings": scalings, "units": units, **(ts_args or {})}
+        topomap_args = {"scalings": scalings, "units": units, **(topomap_args or {})}
+
+        return plot_evoked_joint(
+            stat_evoked,
+            times=times,
+            title=title,
+            picks=picks,
+            exclude=exclude,
+            show=show,
+            ts_args=ts_args,
+            topomap_args=topomap_args,
+        )

mne/stats/tests/test_cluster_level_modern.py

@@ -2,6 +2,9 @@
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
+from functools import partial
+from unittest import mock
+
 import numpy as np
 import pytest
 from numpy.testing import assert_array_almost_equal, assert_array_equal
@@ -12,10 +15,13 @@
     cluster_test,
     combine_adjacency,
     f_mway_rm,
+    f_oneway,
     f_threshold_mway_rm,
     permutation_cluster_1samp_test,
     permutation_cluster_test,
+    ttest_1samp_no_p,
 )
+from mne.stats.cluster_level import ClusterResult, _rm_anova_stat_fun
 from mne.time_frequency import AverageTFRArray, BaseTFR, EpochsTFRArray
 from mne.utils import GetEpochsMixin
 
@@ -419,3 +425,57 @@ def test_cluster_test_adjacency():
     )
     with pytest.raises(ValueError, match='"auto" could not infer sensor adjacency'):
         cluster_test(df_no_montage, "data ~ condition", **kwargs)
+
+
+@pytest.mark.parametrize(
+    ("stat_fun", "unit"),
+    [
+        (ttest_1samp_no_p, "t"),
+        (f_oneway, "F"),
+        pytest.param(
+            partial(_rm_anova_stat_fun, factor_levels=[2, 2], effects="A:B"),
+            "F",
+            id="rm_anova-F",
+        ),
+    ],
+)
+def test_cluster_result_plot_joint_units(stat_fun, unit):
+    """Test plot_joint labels the statistic and does not rescale it to µV."""
+    info = create_info(["Fz", "Cz", "Pz", "Oz"], 100, "eeg")
+    result = ClusterResult(
+        stat_obs=np.zeros((4, 10)),  # required keywords
+        clusters=[],
+        cluster_p_values=np.array([]),
+        H0=np.array([]),
+        stat_fun=stat_fun,
+        n_permutations=1,
+    )
+    # mock a figure
+    with mock.patch("mne.viz.evoked.plot_evoked_joint") as plot:
+        result.plot_joint(info, np.arange(10) / 100, times=[0.02])
+    for key in ("ts_args", "topomap_args"):
+        assert plot.call_args.kwargs[key]["units"]["eeg"] == unit
+        assert plot.call_args.kwargs[key]["scalings"] == dict(eeg=1.0)
+
+
+def test_cluster_result_plot_joint_errors():
+    """Test plot_joint rejects input it cannot plot."""
+    info = create_info(["Fz", "Cz", "Pz", "Oz"], 100, "eeg")
+    kwargs = dict(
+        clusters=[],
+        cluster_p_values=np.array([]),
+        H0=np.array([]),
+        stat_fun=ttest_1samp_no_p,
+        n_permutations=1,
+    )
+    times = np.arange(10) / 100.0
+    tfr = ClusterResult(stat_obs=np.zeros((4, 3, 10)), **kwargs)
+    with pytest.raises(ValueError, match="Time-frequency results"):
+        tfr.plot_joint(info, times)
+    result = ClusterResult(stat_obs=np.zeros((4, 10)), **kwargs)
+    with pytest.raises(ValueError, match="info has 2 channels"):
+        result.plot_joint(create_info(["Fz", "Cz"], 100, "eeg"), times)
+    with pytest.raises(ValueError, match="times_sec does not match"):
+        result.plot_joint(info, times[:5])
+    with pytest.raises(ValueError, match="times_sec does not match"):
+        result.plot_joint(info, np.arange(10) / 50.0)  # wrong sfreq

tutorials/stats-sensor-space/76_new_cluster_test_api.py

@@ -28,8 +28,6 @@
 
 from pathlib import Path
 
-import matplotlib.pyplot as plt
-import numpy as np
 import pandas as pd
 
 import mne
@@ -172,35 +170,18 @@
 # -------------------
 #
 # The result object carries the observed cluster-level statistics. We plot the
-# observed t-values as a heatmap with time on the x-axis and channel names on
-# the y-axis. Because the contrast is target minus non-target, positive t-values
-# mean a stronger response to targets, matching the difference plotted earlier.
+# observed t-values as a joint plot: a butterfly plot of the t-value time course
+# of every channel, with topographic maps at the time points of the largest
+# peaks. Because the contrast is target minus non-target, positive t-values mean
+# a stronger response to targets, matching the difference plotted earlier.
+# Note that this shows the test statistic itself, not which channels and time
+# points belong to a significant cluster.
 
 print(f"Number of permutations run: {cluster_result.n_permutations}")
 
-# times (in seconds) and channel names come from the evoked data
-times = grand_avg_diff.times
-ch_names = grand_avg_diff.ch_names
-
-# stat_obs holds the observed t-values, arranged as (channels, times)
-stat_obs = cluster_result.stat_obs
-
-# symmetric colour limits so the diverging colormap is centred on zero
-vlim = np.abs(stat_obs).max()
-
-fig, ax = plt.subplots(layout="constrained")
-im = ax.imshow(
-    stat_obs,
-    aspect="auto",
-    origin="lower",
-    extent=[times[0], times[-1], 0, len(ch_names)],
-    cmap="RdBu_r",
-    vmin=-vlim,
-    vmax=vlim,
+# Plot cluster result statistic (for now without the actual cluster highlighted)
+fig = cluster_result.plot_joint(
+    info=df["evoked"][0].info,
+    times_sec=df["evoked"][0].times,
+    show=False,
 )
-ax.set_yticks(np.arange(len(ch_names)) + 0.5)
-ax.set_yticklabels(ch_names)
-ax.set_xlabel("time (s)")
-ax.set_ylabel("channel")
-ax.set_title("Observed cluster statistic (target - non-target)")
-fig.colorbar(im, ax=ax, label="t-value")