diff --git a/doc/changes/dev/14399.newfeature.rst b/doc/changes/dev/14399.newfeature.rst new file mode 100644 index 00000000000..f5efc3a6169 --- /dev/null +++ b/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`_. \ No newline at end of file diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index 33f2e9160e5..e0e3a4cb16a 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -2903,3 +2903,128 @@ def __repr__(self): # noqa: D105 f"" ) + + @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, + ) diff --git a/mne/stats/tests/test_cluster_level_modern.py b/mne/stats/tests/test_cluster_level_modern.py index f1c5c7cd92b..09d9b317acc 100644 --- a/mne/stats/tests/test_cluster_level_modern.py +++ b/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 diff --git a/tutorials/stats-sensor-space/76_new_cluster_test_api.py b/tutorials/stats-sensor-space/76_new_cluster_test_api.py index 3faccb8e6ae..8338403acff 100644 --- a/tutorials/stats-sensor-space/76_new_cluster_test_api.py +++ b/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")