diff --git a/doc/api/statistics.rst b/doc/api/statistics.rst index c6f174152b3..f55a57e7d18 100644 --- a/doc/api/statistics.rst +++ b/doc/api/statistics.rst @@ -39,6 +39,8 @@ Non-parametric (clustering) resampling methods: .. autosummary:: :toctree: ../generated/ + ClusterResult + cluster_test combine_adjacency permutation_cluster_test permutation_cluster_1samp_test diff --git a/doc/changes/dev/12663.newfeature.rst b/doc/changes/dev/12663.newfeature.rst new file mode 100644 index 00000000000..83ca873d24e --- /dev/null +++ b/doc/changes/dev/12663.newfeature.rst @@ -0,0 +1 @@ +Add new API for cluster permutation statistics of sensor data: :func:`mne.stats.cluster_test`, see the tutorial :ref:`tut-new-cluster-test-api`, by `Carina Forster`_, `Sophie Herbst`_, :newcontrib:`Maximilien Chaumon`, and `Scott Huberty`_. \ No newline at end of file diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 8a10bc6d8bb..774f35c2916 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -297,6 +297,7 @@ .. _Matti Toivonen: https://github.com/mattitoi .. _Maureen Shader: https://github.com/mshader .. _Mauricio Cespedes Tenorio: https://github.com/mcespedes99 +.. _Maximilien Chaumon: https://github.com/dnacombo .. _Melih Yayli: https://github.com/yaylim .. _Michael Krause: https://github.com/octomike .. _Michael Straube: https://github.com/mistraube diff --git a/doc/conf.py b/doc/conf.py index f4eb8eda8f6..7931108b0eb 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -326,6 +326,7 @@ "RawPersyst": "mne.io.Raw", "RawSNIRF": "mne.io.Raw", "Calibration": "mne.preprocessing.eyetracking.Calibration", + "ClusterResult": "mne.stats.ClusterResult", # dipy "dipy.align.AffineMap": "dipy.align.imaffine.AffineMap", "dipy.align.DiffeomorphicMap": "dipy.align.imwarp.DiffeomorphicMap", diff --git a/environment.yml b/environment.yml index cfd278d4599..158d8b97f96 100644 --- a/environment.yml +++ b/environment.yml @@ -15,6 +15,7 @@ dependencies: - eeglabio - ffmpeg ==8.1.2 - filelock >=3.18.0 + - formulaic - h5io >=0.2.4 - h5py >=2.4 - imageio >=2.6.1 diff --git a/mne/datasets/config.py b/mne/datasets/config.py index 293c29bcbb3..e71452d6f9c 100644 --- a/mne/datasets/config.py +++ b/mne/datasets/config.py @@ -88,7 +88,7 @@ # here: ↓↓↓↓↓↓↓↓ RELEASES = dict( testing="0.176", - misc="0.27", + misc="0.30", phantom_kit="0.2", ucl_opm_auditory="0.2", ) @@ -129,7 +129,7 @@ ) MNE_DATASETS["misc"] = dict( archive_name=f"{MISC_VERSIONED}.tar.gz", # 'mne-misc-data', - hash="md5:e343d3a00cb49f8a2f719d14f4758afe", + hash="md5:201d35531d3c03701cf50e38bb73481f", url=( f"https://codeload.github.com/mne-tools/mne-misc-data/tar.gz/{RELEASES['misc']}" ), diff --git a/mne/stats/__init__.pyi b/mne/stats/__init__.pyi index a206a608eab..d52719e9613 100644 --- a/mne/stats/__init__.pyi +++ b/mne/stats/__init__.pyi @@ -1,9 +1,11 @@ __all__ = [ + "ClusterResult", "_ci", "_parametric_ci", "_st_mask_from_s_inds", "bonferroni_correction", "bootstrap_confidence_interval", + "cluster_test", "combine_adjacency", "erp", "f_mway_rm", @@ -24,7 +26,9 @@ __all__ = [ from . import erp from ._adjacency import combine_adjacency from .cluster_level import ( + ClusterResult, _st_mask_from_s_inds, + cluster_test, permutation_cluster_1samp_test, permutation_cluster_test, spatio_temporal_cluster_1samp_test, diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index 858d259d7d5..36363ccad62 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -4,23 +4,43 @@ # License: BSD-3-Clause # Copyright the MNE-Python contributors. +from __future__ import annotations + +from functools import partial +from string import ascii_uppercase +from typing import TYPE_CHECKING, Literal + import numpy as np +from ..epochs import BaseEpochs +from ..evoked import Evoked, combine_evoked from ..parallel import parallel_func from ..source_estimate import MixedSourceEstimate, SourceEstimate, VolSourceEstimate from ..source_space import SourceSpaces +from ..time_frequency import BaseTFR, EpochsTFR from ..utils import ( + GetEpochsMixin, ProgressBar, _check_option, + _check_rng, _legacy_rng, _pl, + _soft_import, _validate_type, + legacy, logger, split_list, verbose, warn, ) -from .parametric import f_oneway, ttest_1samp_no_p +from .parametric import f_mway_rm, f_oneway, f_threshold_mway_rm, ttest_1samp_no_p + +if TYPE_CHECKING: + from scipy import sparse # Used in type hints for cluster_test + +# need this at top-level of file due to type hints +pd = _soft_import("pandas", purpose="DataFrame integration", strict=False) +DataFrame = getattr(pd, "DataFrame", None) def _get_labels_st(x_in, adjacency, max_step): @@ -840,12 +860,17 @@ def _permutation_cluster_test( out_type, check_disjoint, buffer_size, + within_subject=False, ): """Aux Function. Note. X is required to be a list. Depending on the length of X either a 1 sample t-test or an F test / more sample permutation scheme is elicited. + + ``within_subject=True`` restricts multi-group permutations to swapping each + subject's observations across the groups (repeated-measures designs); rows + of each element of X must then be aligned by subject. """ _check_option("out_type", out_type, ["mask", "indices"]) _check_option("tail", tail, [-1, 0, 1]) @@ -871,7 +896,7 @@ def _permutation_cluster_test( sample_shape = X[0].shape[1:] for x in X: if x.shape[1:] != sample_shape: - raise ValueError("All samples mush have the same size") + raise ValueError("All samples must have the same size") # flatten the last dimensions in case the data is high dimensional X = [np.reshape(x, (x.shape[0], -1)) for x in X] @@ -985,7 +1010,28 @@ def _permutation_cluster_test( n_samples_per_condition = [x.shape[0] for x in X] splits_idx = np.append([0], np.cumsum(n_samples_per_condition)) slices = [slice(splits_idx[k], splits_idx[k + 1]) for k in range(len(X))] - orders = [rng.permutation(len(X_full)) for _ in range(n_permutations - 1)] + if within_subject: + # Repeated-measures design: permute each subject's observations + # only across the conditions (cells), never across subjects -- the + # exchangeability assumption for repeated measures (FieldTrip's + # depsamples* statistics permute the same way). + n_cells, n_subjects = len(X), len(X[0]) + assert all(len(x) == n_subjects for x in X) # checked by callers + # a random permutation of the cells per (permutation, subject) + cell_orders = np.argsort( + rng.uniform(size=(n_permutations - 1, n_subjects, n_cells)), axis=-1 + ) + # the row index of (cell j, subject s) in X_full is + # j * n_subjects + s, so position (j, s) draws from row + # (cell_orders[:, s, j], s): + orders = list( + ( + cell_orders.transpose(0, 2, 1) * n_subjects + + np.arange(n_subjects)[np.newaxis, np.newaxis] + ).reshape(n_permutations - 1, -1) + ) + else: + orders = [rng.permutation(len(X_full)) for _ in range(n_permutations - 1)] del rng parallel, my_do_perm_func, n_jobs = parallel_func( do_perm_func, n_jobs, verbose=False @@ -1073,7 +1119,23 @@ def _permutation_cluster_test( return t_obs, clusters, cluster_pv, H0 -def _check_fun(X, stat_fun, threshold, tail=0, kind="within"): +def _rm_anova_stat_fun(*X, factor_levels, effects): + """Wrap `f_mway_rm` for use as a cluster-test ``stat_fun``. + + ``X`` arrives as one 2D array (replications x flattened locations) per cell of + the design, ordered so that the first factor varies slowest (matching how + :func:`pandas.DataFrame.groupby` orders a multi-column group-by, and what + :func:`~mne.stats.f_mway_rm` expects). + """ + data = np.stack(X, axis=1) # subjects x conditions x locations + return f_mway_rm( + data, factor_levels=factor_levels, effects=effects, return_pvals=False + )[0] + + +def _check_fun( + X, stat_fun, threshold, tail=0, kind="within", factor_levels=None, effects=None +): """Check the stat_fun and threshold values.""" from scipy.stats import f as fstat from scipy.stats import t as tstat @@ -1092,6 +1154,22 @@ def _check_fun(X, stat_fun, threshold, tail=0, kind="within"): threshold = -threshold logger.info(f"Using a threshold of {threshold:.6f}") stat_fun = ttest_1samp_no_p if stat_fun is None else stat_fun + elif kind == "within_rm": + n_subjects = len(X[0]) + if threshold is None: + if stat_fun is not None: + warn( + "Automatic threshold is only valid for stat_fun=None " + f"(uses f_mway_rm internally), got {stat_fun}" + ) + elif tail != 1: + warn('Ignoring argument "tail", performing 1-tailed F-test') + threshold = f_threshold_mway_rm(n_subjects, factor_levels, effects=effects) + logger.info(f"Using a threshold of {threshold:.6f}") + if stat_fun is None: + stat_fun = partial( + _rm_anova_stat_fun, factor_levels=factor_levels, effects=effects + ) else: assert kind == "between" if threshold is None: @@ -1111,6 +1189,7 @@ def _check_fun(X, stat_fun, threshold, tail=0, kind="within"): return stat_fun, threshold +@legacy(alt="mne.stats.cluster_test(...)") @_legacy_rng("seed") @verbose def permutation_cluster_test( @@ -1688,3 +1767,529 @@ def summarize_clusters_stc( data_summary[:, 0] = np.sum(data_summary, axis=1) return klass(data_summary, vertices, tmin, tstep, subject) + + +def _validate_cluster_df(df: DataFrame, dv_name: str, iv_names: list[str]): + """Validate the input DataFrame for cluster tests.""" + # check if all necessary columns are present + missing = ({dv_name} | set(iv_names)) - set(df.columns) # should be empty + sep = '", "' + if missing: # if not empty, there are missing columns + raise ValueError( + f"DataFrame must contain a column named for each term in `formula`. " + f"Column{_pl(missing)} missing for term{_pl(missing)} " # _pl = pluralize + f'"{sep.join(missing)}".' + ) + # check if the data column contains valid (and consistent) instance types + inst = df[dv_name].iloc[0] + valid_types = ( + Evoked, + BaseEpochs, + BaseTFR, + np.ndarray, + ) # Base covers all Epochs and TFRs + _validate_type(inst, valid_types, f"Data in dependent variable column '{dv_name}'") + all_types = set(df[dv_name].map(type)) + all_type_names = ", ".join([type(x).__name__ for x in all_types]) + prologue = f"Data in dependent variable column '{dv_name}' must all have " + if len(all_types) > 1: + raise ValueError( + f"{prologue} the same type, but found types {{{all_type_names}}}." + ) + # check if the shape of the data is consistent + if isinstance(inst, np.ndarray): + all_shapes = set( + df[dv_name].map(lambda x: x.shape[1:]) + ) # first dim may vary (participants or epochs) + elif isinstance(inst, (BaseEpochs | EpochsTFR)): + all_shapes = set(df[dv_name].map(lambda x: x.get_data().shape[1:])) + else: + all_shapes = set(df[dv_name].map(lambda x: x.get_data().shape)) + if len(all_shapes) > 1: + raise ValueError( + f"{prologue} consistent shape, but {len(all_shapes)} different " + f"shapes were found: {'; '.join(all_shapes)}." + ) + obj_type = all_types.pop() + is_epo = GetEpochsMixin in obj_type.__mro__ + is_tfr = BaseTFR in obj_type.__mro__ + is_arr = np.ndarray in obj_type.__mro__ + return is_epo, is_tfr, is_arr + + +# TODO: design/analysis features FieldTrip's cluster stats support that +# cluster_test does not (yet): +# - continuous predictors / regression & correlation designs +# (ft_statfun_indepsamplesregrT, _depsamplesregrT, _correlationT); the +# formula right-hand side currently must be categorical +# - multivariate within-subject F across conditions +# (ft_statfun_depsamplesFmultivariate) +# - activation-versus-baseline tests (ft_statfun_actvsblT) +# - control variables / stratified or blocked resampling (cfg.cvar, cfg.wvar) +# - requiring a minimum number of neighboring channels for cluster membership +# (cfg.minnbchan) +# - the weighted cluster mass statistic (cfg.clusterstatistic='wcm'); +# ``t_power`` covers maxsum (t_power=1) and maxsize (t_power=0) only +@verbose +def cluster_test( + df: DataFrame, + formula: str, + *, # end of positional-only parameters + within_id: str | None = None, + reference: str | None = None, + stat_fun: callable | None = None, + tail: Literal[-1, 0, 1] = 0, + threshold=None, + n_permutations: str | int = 1024, + adjacency: sparse.spmatrix + | None + | Literal[False] = None, # should be None (default) + max_step: int = 1, # TODO may need to provide `max_step_time` and `max_step_freq` + exclude: list | None = None, # TODO needs rethink because user passes MNE objects + step_down_p: float = 0.0, + t_power: float = 1.0, + check_disjoint: bool = False, + out_type: Literal["indices", "mask"] = "indices", + rng: None | int | np.random.Generator | np.random.RandomState = None, + buffer_size: int | None = None, + n_jobs: int = 1, + verbose=None, +): + """Run a cluster permutation test from a DataFrame and a formula. + + Parameters + ---------- + df : pandas.DataFrame + Dataframe containing the data, dependent and independent variables. + formula : str + Wilkinson notation formula naming the dependent variable and either a single + independent variable (e.g. ``"data ~ condition"``) or a single interaction + term between two or more independent variables (e.g. ``"data ~ a:b"``, tested + with a repeated-measures ANOVA; see ``within_id``). All names must match + columns in ``df``. Testing several effects (e.g. two main effects, or a main + effect and an interaction) requires calling :func:`cluster_test` once per + effect. + within_id : None | str + Name of column in ``df`` to use in identifying within-group contrasts. + + - If ``within_id`` is not ``None``: + ``within_id`` must match a column name in ``df``, e.g. ``"subject_index"`` + (a name not in ``df.columns`` will result in an error). If the independent + variable has 1 level per participant, the data will be treated as + already subtracted (e.g., condition A - condition B) and a paired t-test + against zero will be performed (using + :func:`mne.stats.ttest_1samp_no_p`). If the independent + variable has 2 levels, the data will be subtracted for each participant + (e.g., condition A - condition B) first. If it has more than 2 levels, + a one-way repeated-measures ANOVA is performed (using + :func:`mne.stats.f_mway_rm`), with permutations swapping each + subject's observations across the levels (never across subjects). + + - If ``within_id`` is ``None``: + Will perform a between-group test (using :func:`mne.stats.f_oneway`; This + works for 2 levels or more). + + - This parameter is required if: + ``formula``'s right-hand side is an interaction term (e.g. + ``"data ~ a:b"``), in which case each combination of ``within_id`` and the + factors must appear exactly once (a fully balanced repeated-measures + design). + reference : str | None + Level of the independent variable to treat as the reference, i.e. the level + that is *subtracted*. The test statistic is then computed on + ``other_level - reference``, so positive values mean the other level is + larger. Only valid for paired two-level contrasts (a single factor with 2 + levels, with ``within_id`` given); for F-tests and repeated-measures ANOVAs + the statistic is sign-invariant and passing ``reference`` raises an error. + If ``None`` (default), levels are taken in sorted order (or in category + order if the column is a :class:`pandas.Categorical`) and the second one is + the reference. + %(stat_fun_clust_both)s + %(tail_clust)s + %(threshold_clust_both)s + %(n_permutations_clust_all)s + %(adjacency_clust_both)s + max_step : int + Maximum distance between samples (time points). Default is 1. + exclude : array-like of bool | None + Mask to apply to the data to exclude certain points from clustering + (e.g., medial wall vertices). Should be the same shape as the channels/vertices + dimension of the data objects. If ``None``, no points are excluded. + %(step_down_p_clust)s + %(t_power_clust)s + check_disjoint : bool + Whether to check if the ``adjacency`` matrix can be separated into disjoint + sets before clustering. This may lead to faster clustering, especially if + the "time" and/or "frequency" dimensions are large. + out_type : 'mask' | 'indices' + Format used to represent each cluster in the list of clusters stored in + the ``clusters`` attribute of :class:`mne.stats.ClusterResult`: + + - ``'mask'``:s + Each cluster is represented by a boolean array of the same shape as + the ``stat_obs`` attribute array of :class:`mne.stats.ClusterResult`, + with ``True`` values indicating locations that are part of a cluster. Note + that MNE-Python's legacy API + (e.g. :func:`mne.stats.permutation_cluster_test`) would return slices if the + shape is 1D and adjacency is ``None``, whereas ``cluster_test`` will always + return a boolean array. + + - ``'indices'``: + Each cluster is represented by a tuple of 1D integer arrays, one array per + dimension of the array in the ``stat_obs`` attribute of + :class:`mne.stats.ClusterResult`. The arrays + together give the coordinates of all locations belonging to the cluster and + can be used to index ``stat_obs``. + Note that for large datasets, ``'indices'`` may use far less memory than + ``'mask'``. + %(rng)s + buffer_size : int | None + Block size to use when computing test statistics. This can significantly + reduce memory usage when ``n_jobs > 1`` and memory sharing between + processes is enabled (see :func:`mne.set_cache_dir`), because the data will be + shared between processes and each process only needs to allocate space for + a small block of locations at a time. + %(n_jobs)s + %(verbose)s + + Returns + ------- + mne.stats.ClusterResult + Object containing the results of the cluster permutation test. + + Notes + ----- + %(threshold_clust_t_or_f_notes)s + + .. versionadded:: 1.13 + """ + # parse formula + formulaic = _soft_import("formulaic", purpose="parse formula for clustering") + parser = formulaic.parser.DefaultFormulaParser(include_intercept=False) + rng = _check_rng(rng) + + formula_str = formula + formula = formulaic.Formula(formula, _parser=parser) + # extract the dependent variable name + dv_name = str(formula.lhs) + # the right-hand side must be a single term: either one factor (main effect, + # e.g. "a") or a single interaction between factors (e.g. "a:b") + rhs_terms = list(formula.rhs) + if len(rhs_terms) != 1: + raise ValueError( + "the right-hand side of `formula` must be a single term: either one " + 'factor (e.g. "data ~ a") or a single interaction (e.g. "data ~ a:b"). ' + f'Got "{formula.rhs}", which has {len(rhs_terms)} terms. To test ' + "several effects, call `cluster_test` once per effect." + ) + factor_names = [str(factor) for factor in rhs_terms[0].factors] + is_interaction = len(factor_names) > 1 + iv_name = factor_names[0] if not is_interaction else ":".join(factor_names) + + # validate the input dataframe and return the type of the data column entries + is_epo, is_tfr, is_arr = _validate_cluster_df(df, dv_name, factor_names) + + _validate_type(within_id, (str, None), "within_id") + if within_id is not None and within_id not in df.columns: + raise ValueError( + f"within_id must be one of {list(df.columns)}, got {within_id!r}" + ) + + # check if within_id has 1 or 2 levels to do paired t-test (within) + if is_interaction and within_id is None: + raise ValueError( + f'testing the interaction "{iv_name}" requires repeated-measures data; ' + "pass `within_id` naming the column that identifies each subject/" + "replication." + ) + # for within-subject designs, check that each subject has one observation per + # combination of factor(s) (2 for a simple paired test; more for a one-way + # repeated-measures ANOVA or an interaction) + n_groups = df[factor_names].drop_duplicates().shape[0] + if within_id and (is_interaction or n_groups >= 2): + df = df.copy(deep=False) # Don't mutate input dataframe row order! + df.sort_values([*factor_names, within_id], inplace=True) + counts = df[within_id].value_counts() + + iv_names = iv_name.split(":") + groups = df[[dv_name, *iv_names, within_id]].groupby([*iv_names, within_id]) + elem = df[dv_name].iloc[0] + # TODO: Support this for other input types e.g. array, epochs, TFR, etc. + if isinstance(elem, Evoked): + reduce = set(df.columns) - set([*iv_names, within_id, dv_name]) + if reduce: + logger.info( + f"To test '{formula_str}', reducing along column(s): {reduce}" + ) + func = {dv_name: lambda evs: combine_evoked(evs.tolist(), weights="nave")} + df = groups.agg(func).reset_index() + + else: + if any(counts != n_groups): + raise ValueError( + f"for a within-subject test, each subject (column {within_id!r}) " + f"must have exactly {n_groups} observations, one per combination " + f"of {factor_names}." + ) + # extract the data from the dataframe + outer_func = np.concatenate if is_epo else np.array + axes = (-3, -1) if is_tfr else (-2, -1) + + def func_arr(series): + return np.concatenate(series.values) + + def func_mne(series): + return outer_func( + series.map(lambda inst: inst.get_data().swapaxes(*axes)).to_list() + ) + + func = func_arr if is_arr else func_mne + + # convert to a list-like X for clustering. Grouping by multiple columns sorts + # lexicographically (first factor varies slowest), which is what f_mway_rm + # expects for interaction effects. + grouped = df.groupby(factor_names, observed=True).agg({dv_name: func})[dv_name] + levels = grouped.index.to_list() # parallel to X by construction + X = grouped.to_list() + contrast = None # set below if a subtraction is performed + + _validate_type(reference, (str, None), "reference") + if reference is not None: + if is_interaction or within_id is None or len(levels) != 2: + raise ValueError( + "`reference` only applies to paired two-level contrasts (a single " + "factor with 2 levels, with `within_id` given); for F-tests and " + "repeated-measures ANOVAs the statistic is sign-invariant." + ) + if reference not in levels: + raise ValueError( + f"reference must be one of the levels of {iv_name!r} ({levels}), " + f"got {reference!r}" + ) + if levels.index(reference) == 0: # reference is subtracted → put it last + levels, X = levels[::-1], X[::-1] + + # determine test type. NOTE: branches that set kind="within" also collapse X + # from a list of groups to an ndarray of shape (n_subjects, ...), so `len(X)` + # below means "number of groups" only until that happens. + if is_interaction: + kind = "within_rm" + factor_levels = [df[name].nunique() for name in factor_names] + # f_mway_rm/f_threshold_mway_rm only understand generic "A", "B", ... + # factor labels (in the order given in `formula`), not the actual column + # names, so translate the interaction accordingly. + rm_effects = ":".join(ascii_uppercase[: len(factor_names)]) + elif len(X) == 1: + kind = "within" # single group -- e.g. already-subtracted paired data + X = X[0] + elif within_id is not None and len(X) > 2: + # one within-subject factor with 3+ levels: one-way repeated-measures + # ANOVA (each subject contributes one observation per level) + kind = "within_rm" + factor_levels = [len(X)] + rm_effects = "A" + elif len(X) > 2: + kind = "between" + elif ( + len(set(x.shape for x in X)) > 1 + ): # unequal number of observations in each group + if within_id is not None: + raise ValueError( + "for a within-subject test, all groups must have the same number " + "of observations; check that every subject has data for every " + f"level of {iv_name!r}." + ) + kind = "between" + # by now we know there are exactly 2 elements in X, and their shapes match + elif within_id in df: + kind = "within" + assert len(X) == 2 + contrast = (levels[0], levels[1]) + logger.info( + f"Subtracting ({levels[0]} - {levels[1]}) of column {iv_name!r} before " + "computing cluster statistics." + ) + X = X[0] - X[1] + else: # 2 elements in X but no within_id provided → unpaired test + kind = "between" + + # define stat function and threshold + if kind == "within_rm": + stat_fun, threshold = _check_fun( + X=X, + stat_fun=stat_fun, + threshold=threshold, + tail=tail, + kind=kind, + factor_levels=factor_levels, + effects=rm_effects, + ) + else: + stat_fun, threshold = _check_fun( + X=X, stat_fun=stat_fun, threshold=threshold, tail=tail, kind=kind + ) + + # check_fun doesn't work with list input` + if kind == "within": # will this create an issue for already subtracted data? + X = [X] + + kind_descs = { + "between": "between-groups F-test", + "within": "one-sample T-test", + "within_rm": "M-way repeated measures ANOVA", + } + func_name = stat_fun.__name__ if "__name__" in dir(stat_fun) else str(stat_fun) + logger.info(f"Chosen statistic: {kind_descs[kind]} -- {func_name}") + + # Run the cluster-based permutation test + stat_obs, clusters, cluster_p_values, H0 = _permutation_cluster_test( + X, + n_permutations=n_permutations, + threshold=threshold, + stat_fun=stat_fun, + tail=tail, + n_jobs=n_jobs, + adjacency=adjacency, + max_step=max_step, # maximum distance between samples (time points) + exclude=exclude, # exclude no time points or channels + step_down_p=step_down_p, # step down in jumps test + t_power=t_power, # weigh each location by its stats score + out_type=out_type, + check_disjoint=check_disjoint, + buffer_size=buffer_size, # block size for chunking the data + rng=rng, + # repeated-measures ANOVA: permute within subjects only + within_subject=kind == "within_rm", + ) + + stat_obs = stat_obs.T + if out_type == "mask": + if isinstance(clusters[0], np.ndarray) and clusters[0].dtype == "bool": + clusters = [cl.T for cl in clusters] + elif isinstance(clusters[0], tuple) and isinstance(clusters[0][0], slice): + clusters = [tuple(reversed(cluster)) for cluster in clusters] + # Convert from old form of slices to mask, make users life easier. + new_clusters = list() + for clust in clusters: + new_clust = np.zeros(stat_obs.shape, bool) + new_clust[clust] = True + new_clusters.append(new_clust) + clusters = new_clusters + elif out_type == "indices": + clusters = [tuple(reversed(cluster)) for cluster in clusters] + return ClusterResult( + stat_obs=stat_obs, + clusters=clusters, + cluster_p_values=cluster_p_values, + H0=H0, + stat_fun=stat_fun, + n_permutations=n_permutations, + t_power=t_power, + contrast=contrast, + ) + + +def _cluster_mass(stat_obs, cluster, t_power): + """Compute a cluster's mass, matching _find_clusters_1dir's own formula.""" + vals = stat_obs[cluster] + if t_power == 1: + return vals.sum() + return (np.sign(vals) * np.abs(vals) ** t_power).sum() + + +class ClusterResult: + """Object containing the results of the cluster permutation test. + + .. note:: + This class is not meant to be instantiated directly, but rather returned + by :func:`~mne.stats.cluster_test`. + + Parameters + ---------- + stat_obs : np.ndarray + The observed test statistic. + clusters : list + List of clusters. + cluster_p_values : np.ndarray + P-values for each cluster. + H0 : np.ndarray + Max cluster level stats observed under permutation. + stat_fun : callable | None + Function called to calculate the test statistic. Must accept 1D-array as + input and return a 1D array. If ``None`` (the default), uses + :func:`mne.stats.ttest_1samp_no_p` for paired tests and + :func:`mne.stats.f_oneway` for unpaired tests or tests of more than 2 groups. + n_permutations : int + The number of permutations that were taken to compute the test statistic. + t_power : float + Power to which the observed statistic was raised (sign retained) before + summing within a cluster to obtain its mass (see ``cluster_masses``). + Should match whatever ``t_power`` was passed to :func:`cluster_test`. + contrast : tuple of str | None + The two levels that were contrasted, as ``(positive, reference)``; the data + were computed as the first minus the second. ``None`` when no subtraction + was performed. + + Attributes + ---------- + cluster_masses : np.ndarray + The mass of each cluster, i.e. the sum (optionally ``t_power``-weighted) + of ``stat_obs`` within that cluster. This is the same per-cluster + statistic that is compared against the permutation distribution (``H0``) + to obtain ``cluster_p_values``, so it is a natural way to rank clusters by + how extreme they are, independent of the resulting p-value. + reference : str | None + The level that was subtracted, i.e. ``contrast[1]``, or ``None``. + + Notes + ----- + .. versionadded:: 1.13 + """ + + def __init__( + self, + *, + stat_obs: np.typing.NDArray, + clusters: list, + cluster_p_values: np.typing.NDArray, + H0: np.typing.NDArray, + stat_fun: callable, + n_permutations: int, + t_power: float = 1.0, + contrast: tuple | None = None, + ): + self.stat_obs = stat_obs + self.clusters = clusters + self.cluster_p_values = cluster_p_values + self.H0 = H0 + self.stat_fun = stat_fun + self.t_power = t_power + self.cluster_masses = np.array( + [_cluster_mass(stat_obs, c, t_power) for c in clusters] + ) + self.n_permutations = n_permutations + self.contrast = contrast + self.reference = None if contrast is None else contrast[1] + + # unpaired t-test equivalent to f_oneway w/ 2 groups + if stat_fun is f_oneway: + self.stat_name = "F-statistic" + elif stat_fun is ttest_1samp_no_p: + self.stat_name = "paired T-statistic" + if contrast is not None: + self.stat_name += f" ({contrast[0]} - {contrast[1]})" + elif isinstance(stat_fun, partial) and stat_fun.func is _rm_anova_stat_fun: + self.stat_name = "F-statistic (repeated-measures ANOVA)" + else: + self.stat_name = "test statistic" + + def __repr__(self): # noqa: D105 + contrast = ( + "" + if self.contrast is None + else f", {self.contrast[0]} - {self.contrast[1]}" + ) + return ( + f"" + ) diff --git a/mne/stats/tests/conftest.py b/mne/stats/tests/conftest.py new file mode 100644 index 00000000000..6389f862781 --- /dev/null +++ b/mne/stats/tests/conftest.py @@ -0,0 +1,36 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import numpy as np +import pytest + + +def _get_conditions(): + n_space = 50 + noise_level = 20 + n_time_1 = 20 + n_time_2 = 13 + normfactor = np.hanning(20).sum() + rng = np.random.default_rng(42) + condition1_1d = rng.normal(scale=noise_level, size=(n_time_1, n_space)) + for c in condition1_1d: + c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor + + condition2_1d = rng.normal(scale=noise_level, size=(n_time_2, n_space)) + for c in condition2_1d: + c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor + + pseudoekp = 10 * np.hanning(25)[None, :] + condition1_1d[:, 25:] += pseudoekp + condition2_1d[:, 25:] -= pseudoekp + + condition1_2d = condition1_1d[:, :, np.newaxis] + condition2_2d = condition2_1d[:, :, np.newaxis] + return condition1_1d, condition2_1d, condition1_2d, condition2_2d + + +@pytest.fixture(scope="session") +def stat_conditions(): + """Get data for mne.stats.tests.""" + return _get_conditions() diff --git a/mne/stats/tests/test_cluster_equiv.py b/mne/stats/tests/test_cluster_equiv.py new file mode 100644 index 00000000000..59679a821b0 --- /dev/null +++ b/mne/stats/tests/test_cluster_equiv.py @@ -0,0 +1,365 @@ +"""FieldTrip equivalence tests for :func:`mne.stats.cluster_test`. + +Reference values were computed with FieldTrip (fieldtrip-20260812) in MATLAB R2026a: +https://gist.github.com/larsoner/5c99b464bccf67f5641c1a2babc2c84e +""" + +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import numpy as np +import pytest +from scipy import sparse, stats + +import mne +from mne.stats import cluster_test + +pd = pytest.importorskip("pandas") +pytest.importorskip("formulaic") + + +def _bump(n_times, sl, amp): + """Return a length-n_times signal with a Hann bump of peak ``amp`` in ``sl``.""" + out = np.zeros(n_times) + width = sl.stop - sl.start + out[sl] = amp * np.hanning(width + 2)[1:-1] # no zero-amplitude edge samples + return out + + +def _within_rows(rng, cond_signals, n_subjects, cols=("condition",), offsets=True): + """Make one df row per subject and condition (with random intercepts).""" + offs = rng.standard_normal(n_subjects) if offsets else np.zeros(n_subjects) + rows = [] + for si in range(n_subjects): + for cond, sig in cond_signals.items(): + data = sig + offs[si] + rng.standard_normal(sig.shape) + row = dict(data=data[np.newaxis], subject=si) + row.update(zip(cols, (cond,) if isinstance(cond, str) else cond)) + rows.append(row) + return pd.DataFrame(rows) + + +def _between_rows(rng, group_signals, group_ns): + """Make ``group_ns[group]`` df rows per independent group.""" + rows = [ + dict(data=(sig + rng.standard_normal(sig.shape))[np.newaxis], group=group) + for group, sig in group_signals.items() + for _ in range(group_ns[group]) + ] + return pd.DataFrame(rows) + + +def _kwargs(**overrides): + """Paired-test defaults; n_permutations > 2**(n-1) -> exact sign flips.""" + out = dict( + formula="data ~ condition", + within_id="subject", + tail=0, + n_permutations=4096, + rng=0, + ) + out.update(overrides) + return out + + +# Scenario builders return (df, kwargs-for-cluster_test); the gist's .mat +# exporter imports these same builders (via SCENARIOS) for the FieldTrip runs. + + +def scenario_one_sample(): + """Pre-subtracted diffs, one pos bump; FT: depsamplesT vs. all-zero partner.""" + rng = np.random.default_rng(1) + n_subjects, n_times = 12, 30 + signals = {"diff": _bump(n_times, slice(8, 15), 1.2)} + df = _within_rows(rng, signals, n_subjects, offsets=False) + return df, _kwargs(threshold=stats.t.ppf(1 - 0.025, n_subjects - 1)) + + +def scenario_paired(): + """Paired t, 2 conditions, one pos + one neg bump; FT: depsamplesT.""" + rng = np.random.default_rng(2) + n_subjects, n_times = 12, 30 + signals = { + "a": _bump(n_times, slice(5, 12), 1.5), # a > b early + "b": _bump(n_times, slice(18, 25), 1.5), # a < b late + } + df = _within_rows(rng, signals, n_subjects) + return df, _kwargs(threshold=stats.t.ppf(1 - 0.025, n_subjects - 1)) + + +def scenario_between_t(): + """2 independent groups, unequal n, opposite-sign group differences. + + FT: indepsamplesF, and indepsamplesT (MNE f_oneway F == FT t**2). + """ + rng = np.random.default_rng(3) + n_times, group_ns = 30, {"ctrl": 10, "pat": 12} + signals = { + "ctrl": _bump(n_times, slice(4, 11), 1.8), # ctrl > pat early + "pat": _bump(n_times, slice(17, 24), 1.8), # pat > ctrl late + } + df = _between_rows(rng, signals, group_ns) + threshold = stats.f.ppf(1 - 0.05, 1, sum(group_ns.values()) - 2) + return df, _kwargs( + formula="data ~ group", + within_id=None, + tail=1, + threshold=threshold, + n_permutations=5000, + ) + + +def scenario_between_anova(): + """3 independent groups, unequal n, 2 bumps with different orderings. + + FT: indepsamplesF. + """ + rng = np.random.default_rng(4) + n_times, group_ns = 30, {"g1": 8, "g2": 9, "g3": 10} + means_1 = {"g1": 1.8, "g2": 0.0, "g3": -1.8} # bump 1 group means + means_2 = {"g1": 0.0, "g2": -1.8, "g3": 1.8} # bump 2 group means + signals = { + g: _bump(n_times, slice(3, 10), means_1[g]) + + _bump(n_times, slice(17, 24), means_2[g]) + for g in group_ns + } + df = _between_rows(rng, signals, group_ns) + threshold = stats.f.ppf(1 - 0.05, 2, sum(group_ns.values()) - 3) + return df, _kwargs( + formula="data ~ group", + within_id=None, + tail=1, + threshold=threshold, + n_permutations=5000, + ) + + +def scenario_rm_anova_interaction(): + """2x2 repeated-measures ANOVA interaction, two opposite-sign bumps. + + FT: depsamplesT on the double difference (a1b1 - a1b2) - (a2b1 - a2b2); + the interaction F equals that t**2. + """ + rng = np.random.default_rng(5) + n_subjects, n_times, c = 12, 30, 0.8 + # pure interaction pattern: sign of the bumps per (a, b) cell + signals = { + (a, b): _bump(n_times, slice(6, 13), sgn) + _bump(n_times, slice(19, 26), -sgn) + for (a, b), sgn in { + ("a1", "b1"): +c, + ("a1", "b2"): -c, + ("a2", "b1"): -c, + ("a2", "b2"): +c, + }.items() + } + df = _within_rows(rng, signals, n_subjects, cols=("a", "b")) + threshold = stats.f.ppf(1 - 0.05, 1, n_subjects - 1) + return df, _kwargs( + formula="data ~ a:b", tail=1, threshold=threshold, n_permutations=1024 + ) + + +def scenario_spatiotemporal_paired(): + """Paired t on 4-channel Evoked data with chain adjacency A1-A2-A3-A4. + + Three bumps in the a-minus-b difference: positive on A1+A2 (adjacent -> + must merge) and on A4 (same time window, but not adjacent to A2 because A3 + is clean -> must stay separate), negative on A3 only, later. + FT: depsamplesT with cfg.neighbours encoding the same chain. + """ + rng = np.random.default_rng(6) + n_subjects, n_channels, n_times = 12, 4, 20 + sig_a = np.zeros((n_channels, n_times)) + for ch in (0, 1, 3): + sig_a[ch] += _bump(n_times, slice(3, 9), 1.5) + sig_b = np.zeros((n_channels, n_times)) + sig_b[2] += _bump(n_times, slice(12, 18), 1.5) # a - b negative on A3 + info = mne.create_info([f"A{n + 1}" for n in range(n_channels)], 1000.0, "eeg") + df = _within_rows(rng, {"a": sig_a, "b": sig_b}, n_subjects) + df["data"] = df["data"].map(lambda d: mne.EvokedArray(d[0], info, tmin=0.0)) + adjacency = sparse.coo_array( + np.diag(np.ones(n_channels - 1), 1) + np.diag(np.ones(n_channels - 1), -1) + ) + threshold = stats.t.ppf(1 - 0.025, n_subjects - 1) + return df, _kwargs(threshold=threshold, adjacency=adjacency) + + +def scenario_rm_3level(): + """3-level within factor: one-way rm ANOVA; FT: depsamplesFunivariate.""" + rng = np.random.default_rng(7) + n_subjects, n_times = 12, 30 + means = {"c1": 1.0, "c2": 0.0, "c3": -1.0} + signals = {c: _bump(n_times, slice(10, 17), amp) for c, amp in means.items()} + df = _within_rows(rng, signals, n_subjects) + threshold = stats.f.ppf(1 - 0.05, 2, 2 * (n_subjects - 1)) + return df, _kwargs(tail=1, threshold=threshold, n_permutations=1024) + + +SCENARIOS = { + "one_sample": scenario_one_sample, + "paired": scenario_paired, + "between_t": scenario_between_t, + "between_anova": scenario_between_anova, + "rm_anova_interaction": scenario_rm_anova_interaction, + "spatiotemporal_paired": scenario_spatiotemporal_paired, + "rm_3level": scenario_rm_3level, +} + +# FieldTrip reference values (see gist for the full cfg): ``critval`` is FT's +# cluster-forming threshold (== ``threshold``; squared when ``square=True``, +# i.e. FT ran a signed t where MNE runs the equivalent F = t**2 -- masses are +# then sums of the squared FT t map). ``clusters`` holds (member indices, +# cluster mass, FT prob); ``p_slack`` is tight for exhaustive sign-flip tests, +# loose for Monte Carlo and for differing permutation groups (2x2 interaction). +FT_REF = { + "one_sample": dict( # cfg.statistic='depsamplesT' (data vs. zeros) + stat_name="paired T-statistic", + critval=2.20098516009, + p_slack=0.001, + stat_max=8.30997573613, + stat_sum=15.2992735697, + clusters=[ + ([9, 10, 11, 12], 24.215180039, 0), + ([8], -2.80238763033, 0.218017578125), + ([0], 2.42485110412, 0.376220703125), + ([4], 2.39754436956, 0.38916015625), + ([2], -2.39390450562, 0.39013671875), + ([17], -2.29386021535, 0.4423828125), + ], + ), + "paired": dict( # cfg.statistic='depsamplesT' + stat_name="paired T-statistic", + critval=2.20098516009, + p_slack=0.001, + stat_max=4.17741541515, + stat_sum=-10.4491353518, + clusters=[ + ([19, 20, 21, 22, 23], -20.5133161849, 0), + ([7, 8], 6.48082672584, 0.00830078125), + ], + ), + "between_t": dict( # cfg.statistic='indepsamplesF' + stat_name="F-statistic", + critval=4.35124350333, + p_slack=0.03, + stat_max=21.5844596996, + stat_sum=106.323938106, + clusters=[ + ([6, 7, 8, 9], 44.1166473563, 0.0001999900005), + ([19, 20, 21], 32.99796735, 0.00109994500275), + ([17], 9.91559364543, 0.183390830458), + ], + ), + "between_anova": dict( # cfg.statistic='indepsamplesF' + stat_name="F-statistic", + critval=3.40282610535, + p_slack=0.03, + stat_max=42.9850709813, + stat_sum=251.00244401, + clusters=[ + ([18, 19, 20, 21, 22], 113.413482696, 4.9997500125e-05), + ([5, 6, 7, 8], 111.934649709, 4.9997500125e-05), + ([12], 4.11422783286, 0.566721663917), + ], + ), + # cfg.statistic='depsamplesT' on the per-subject double difference + # (a1b1 - a1b2) - (a2b1 - a2b2) vs. zeros; interaction F = t**2 + "rm_anova_interaction": dict( + stat_name="F-statistic (repeated-measures ANOVA)", + critval=2.20098516009, + square=True, + p_slack=0.05, + stat_max=35.751136081, + stat_sum=227.529952225, + clusters=[ + ([20, 21, 22, 23], 104.752717043, 0), + ([8, 9, 10, 11], 83.9317992007, 0.00048828125), + ([6], 9.38758796849, 0.144287109375), + ], + ), + # cfg.statistic='depsamplesT', cfg.neighbours = chain A1-A2-A3-A4 + "spatiotemporal_paired": dict( + stat_name="paired T-statistic", + critval=2.20098516009, + p_slack=0.001, + stat_max=5.28186348184, + stat_sum=27.7385987319, + clusters=[ + ([(4, 0), (5, 0), (5, 1), (6, 0), (6, 1), (7, 1)], 18.6062044319, 0), + ([(14, 2), (15, 2), (16, 2)], -11.5475200021, 0.001708984375), + ([(5, 3), (6, 3)], 7.58354029468, 0.013671875), + ([(1, 2)], -2.91124038589, 0.4404296875), + ([(16, 0)], 2.83280102942, 0.477783203125), + ([(10, 3)], 2.38278255024, 0.7353515625), + ], + ), + "rm_3level": dict( # cfg.statistic='depsamplesFunivariate' + stat_name="F-statistic (repeated-measures ANOVA)", + critval=3.44335677937, + p_slack=0.03, + stat_max=24.7722284367, + stat_sum=83.3268364113, + clusters=[ + ([12, 13, 14], 50.6420452705, 4.9997500125e-05), + ], + ), +} +# cfg.statistic='indepsamplesT' on the between_t data: same clusters and (via +# F = t**2) masses as indepsamplesF, but a t critval and per-tail probs +_t_probs = (4.9997500125e-05, 0.00029998500075, 0.0896955152242) +FT_REF["between_t_T"] = dict( + FT_REF["between_t"], + critval=2.08596344727, + square=True, + clusters=[(*c[:2], p) for c, p in zip(FT_REF["between_t"]["clusters"], _t_probs)], +) + + +def _mne_clusters(result): + """Map frozenset of cluster members -> (mass, p) for a ClusterResult.""" + out = {} + for ci, cl in enumerate(result.clusters): + cl = cl if isinstance(cl, tuple) else (cl,) + if len(cl) == 1: + members = frozenset(int(i) for i in cl[0]) + else: + # reveersal here to deal with transpoe added after FT reulsts were generated + members = frozenset(zip(*cl[::-1])) + out[members] = (result.cluster_masses[ci], result.cluster_p_values[ci]) + return out + + +def _assert_ft_equiv(result, ref, threshold): + """Assert a ClusterResult matches one FieldTrip reference run.""" + square = ref.get("square", False) + assert result.stat_name == ref["stat_name"] # check test-type routing + # the threshold each scenario passes must equal FT's parametric + # cfg.clustercritval (from cfg.clusteralpha=0.05), or its square + np.testing.assert_allclose(threshold, ref["critval"] ** (1 + square), rtol=1e-9) + np.testing.assert_allclose(result.stat_obs.max(), ref["stat_max"], rtol=1e-6) + np.testing.assert_allclose(result.stat_obs.sum(), ref["stat_sum"], rtol=1e-6) + got = _mne_clusters(result) + expected = { + frozenset(members): (mass, prob) for members, mass, prob in ref["clusters"] + } + assert set(got) == set(expected), (sorted(map(sorted, got)),) + ft_is_t = square or "T-statistic" in ref["stat_name"] + slack = ref["p_slack"] + for members, (ft_mass, ft_prob) in expected.items(): + mass, p = got[members] + np.testing.assert_allclose(mass, ft_mass, rtol=1e-6) + if ft_is_t: # MNE null pools both tails; FT prob is per-tail + assert ft_prob - slack <= p <= 2 * ft_prob + slack, (p, ft_prob) + else: # same one-sided F null on both sides + assert abs(p - ft_prob) <= slack, (p, ft_prob) + + +@pytest.mark.parametrize("name", list(SCENARIOS)) +def test_fieldtrip_equivalence(name): + """Compare cluster_test output against FieldTrip reference values.""" + df, kwargs = SCENARIOS[name]() + result = cluster_test(df, **kwargs, verbose="error") + _assert_ft_equiv(result, FT_REF[name], kwargs["threshold"]) + if name == "between_t": # also equivalent to FT's independent-samples t + _assert_ft_equiv(result, FT_REF["between_t_T"], kwargs["threshold"]) diff --git a/mne/stats/tests/test_cluster_level.py b/mne/stats/tests/test_cluster_level_legacy.py similarity index 96% rename from mne/stats/tests/test_cluster_level.py rename to mne/stats/tests/test_cluster_level_legacy.py index 7019b2703ce..5b148a1439e 100644 --- a/mne/stats/tests/test_cluster_level.py +++ b/mne/stats/tests/test_cluster_level_legacy.py @@ -15,7 +15,12 @@ ) from scipy import linalg, sparse, stats -from mne import MixedSourceEstimate, SourceEstimate, SourceSpaces, VolSourceEstimate +from mne import ( + MixedSourceEstimate, + SourceEstimate, + SourceSpaces, + VolSourceEstimate, +) from mne.stats import combine_adjacency, ttest_ind_no_p from mne.stats.cluster_level import ( _find_clusters, @@ -31,31 +36,6 @@ ) from mne.utils import _record_warnings, catch_logging -n_space = 50 - - -def _get_conditions(): - noise_level = 20 - n_time_1 = 20 - n_time_2 = 13 - normfactor = np.hanning(20).sum() - rng = np.random.default_rng(42) - condition1_1d = rng.normal(scale=noise_level, size=(n_time_1, n_space)) - for c in condition1_1d: - c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor - - condition2_1d = rng.normal(scale=noise_level, size=(n_time_2, n_space)) - for c in condition2_1d: - c[:] = np.convolve(c, np.hanning(20), mode="same") / normfactor - - pseudoekp = 10 * np.hanning(25)[None, :] - condition1_1d[:, 25:] += pseudoekp - condition2_1d[:, 25:] -= pseudoekp - - condition1_2d = condition1_1d[:, :, np.newaxis] - condition2_2d = condition2_1d[:, :, np.newaxis] - return condition1_1d, condition2_1d, condition1_2d, condition2_2d - def test_thresholds(numba_conditional): """Test automatic threshold calculations.""" @@ -212,9 +192,9 @@ def test_permutation_step_down_p(numba_conditional): assert_allclose(p_next, 0.015625, atol=1e-6) -def test_cluster_permutation_test(numba_conditional): +def test_cluster_permutation_test(numba_conditional, stat_conditions): """Test cluster level permutations tests.""" - condition1_1d, condition2_1d, condition1_2d, condition2_2d = _get_conditions() + condition1_1d, condition2_1d, condition1_2d, condition2_2d = stat_conditions for condition1, condition2 in zip( (condition1_1d, condition1_2d), (condition2_1d, condition2_2d) ): @@ -258,9 +238,9 @@ def stat_fun(X, Y): @pytest.mark.parametrize( "stat_fun", [ttest_1samp_no_p, partial(ttest_1samp_no_p, sigma=1e-1)] ) -def test_cluster_permutation_t_test(numba_conditional, stat_fun): +def test_cluster_permutation_t_test(numba_conditional, stat_conditions, stat_fun): """Test cluster level permutations T-test.""" - condition1_1d, _, condition1_2d, _ = _get_conditions() + condition1_1d, _, condition1_2d, _ = stat_conditions # use a very large sigma to make sure Ts are not independent for condition1, p in ((condition1_1d, 0.01), (condition1_2d, 0.01)): @@ -338,14 +318,17 @@ def test_cluster_permutation_t_test(numba_conditional, stat_fun): ) -def test_cluster_permutation_with_adjacency(numba_conditional, monkeypatch): +def test_cluster_permutation_with_adjacency( + numba_conditional, monkeypatch, stat_conditions +): """Test cluster level permutations with adjacency matrix.""" pytest.importorskip("sklearn") from sklearn.feature_extraction.image import grid_to_graph - condition1_1d, condition2_1d, _, _ = _get_conditions() + condition1_1d, condition2_1d, _, _ = stat_conditions n_pts = condition1_1d.shape[1] + n_space = 50 # we don't care about p-values in any of these, so do fewer permutations args = dict( rng=None, @@ -751,12 +734,12 @@ def test_labels_to_clusters(): assert_array_equal(got[0], active) -def test_spatio_temporal_cluster_adjacency(numba_conditional): +def test_spatio_temporal_cluster_adjacency(numba_conditional, stat_conditions): """Test spatio-temporal cluster permutations.""" pytest.importorskip("sklearn") from sklearn.feature_extraction.image import grid_to_graph - condition1_1d, condition2_1d, condition1_2d, condition2_2d = _get_conditions() + condition1_1d, condition2_1d, condition1_2d, condition2_2d = stat_conditions rng = np.random.default_rng(0) noise1_2d = rng.standard_normal( diff --git a/mne/stats/tests/test_cluster_level_modern.py b/mne/stats/tests/test_cluster_level_modern.py new file mode 100644 index 00000000000..f809e012e74 --- /dev/null +++ b/mne/stats/tests/test_cluster_level_modern.py @@ -0,0 +1,306 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal, assert_array_equal + +from mne import EpochsArray, EvokedArray, create_info +from mne.stats import ( + cluster_test, + f_mway_rm, + f_threshold_mway_rm, + permutation_cluster_1samp_test, + permutation_cluster_test, +) +from mne.time_frequency import AverageTFRArray, BaseTFR, EpochsTFRArray +from mne.utils import GetEpochsMixin + +pd = pytest.importorskip("pandas") +pytest.importorskip("formulaic") # required for cluster_test API + + +def _convert_cluster_slices_to_arrays(clusters, stat_obs_shape): + """Old API sometimes returns slices, we always want masked arrays.""" + cluster_masks = list() + for clust in clusters: + clust_mask = np.zeros(stat_obs_shape, bool) + clust_mask[clust] = True + cluster_masks.append(clust_mask) + return cluster_masks + + +def test_cluster_test_one_sample(stat_conditions): + """Test cluster_test with a single-group (1-sample) design.""" + condition1_1d, _, _, _ = stat_conditions + df = pd.DataFrame(dict(data=[condition1_1d], group=["only"])) + kwargs = dict(n_permutations=100, tail=0, rng=1, buffer_size=None) + T_obs, clusters, cluster_pvals, H0 = permutation_cluster_1samp_test( + condition1_1d, **kwargs + ) + result = cluster_test(df, "data ~ group", **kwargs) + assert result.stat_name == "paired T-statistic" + assert_array_equal(result.H0, H0) + assert_array_equal(result.stat_obs, T_obs) + assert_array_equal(result.cluster_p_values, cluster_pvals) + assert len(result.clusters) == len(clusters) + for clu1, clu2 in zip(result.clusters, clusters): + assert_array_equal(clu1, clu2) + + +def test_compare_old_and_new_cluster_api(stat_conditions): + """Test for same results from old and new APIs.""" + condition1_1d, condition2_1d, condition1_2d, condition2_2d = stat_conditions + df_1d = pd.DataFrame( + dict( + data=[condition1_1d, condition2_1d], + condition=["a", "b"], + ) + ) + kwargs = dict(n_permutations=100, tail=1, rng=1, buffer_size=None, out_type="mask") + F_obs, clusters, cluster_pvals, H0 = permutation_cluster_test( + [condition1_1d, condition2_1d], **kwargs + ) + formula = "data ~ condition" + cluster_result = cluster_test(df_1d, formula, **kwargs) + + for clust in cluster_result.clusters: + assert clust.shape == cluster_result.stat_obs.shape + assert_array_equal(cluster_result.H0, H0) + assert_array_equal(cluster_result.stat_obs, F_obs) + assert_array_equal(cluster_result.cluster_p_values, cluster_pvals) + + assert len(clusters) == len(cluster_result.clusters) + # Convert slices to masked arrays + cluster_masks = _convert_cluster_slices_to_arrays(clusters, F_obs.shape) + for cluster, res_clust in zip(cluster_masks, cluster_result.clusters): + # bool_clust = np.zeros(F_obs.shape, bool) + # bool_clust[cluster] = True + np.testing.assert_array_equal(cluster.T, res_clust) + + +@pytest.mark.parametrize( + "Inst", (EpochsArray, EvokedArray, EpochsTFRArray, AverageTFRArray) +) +@pytest.mark.filterwarnings('ignore:Ignoring argument "tail":RuntimeWarning') +def test_new_cluster_api(Inst): + """Test handling different MNE objects in the cluster API.""" + rng = np.random.default_rng(seed=8675309) + is_epo = GetEpochsMixin in Inst.__mro__ + is_tfr = BaseTFR in Inst.__mro__ + + n_epo, n_chan, n_freq, n_times = 6, 3, 4, 5 + + # prepare the dimensions of the simulated data, then simulate + size = (n_chan,) + if is_epo: + size = (n_epo, *size) + if is_tfr: + size = (*size, n_freq) + size = (*size, n_times) + data = rng.normal(size=size) + + # construct the instance + info = create_info(ch_names=n_chan, sfreq=1000, ch_types="eeg") + kw = dict(times=np.arange(n_times), freqs=np.arange(n_freq)) if is_tfr else dict() + cond_a = Inst(data=data, info=info, **kw) + cond_b = cond_a.copy() + # introduce a significant difference in a specific region, time, and frequency + ch_start, ch_end = 0, 2 # 2 channels + t_start, t_end = 2, 4 # 2 times + f_start, f_end = 2, 4 # 2 freqs + if is_tfr: + cond_b._data[..., ch_start:ch_end, f_start:f_end, t_start:t_end] += 2 + else: + cond_b._data[..., ch_start:ch_end, t_start:t_end] += 2 + # for Evokeds/AverageTFRs, we create fake "subjects" as our observations within each + # condition. We add a bit of noise while we do so. + if not is_epo: + insts = list() + for cond in cond_a, cond_b: + for _n in range(n_epo): + if not _n: + insts.append(cond) + continue + _cond = cond.copy() + _cond.data += rng.normal(scale=0.1, size=_cond.data.shape) + insts.append(_cond) + conds = np.repeat(["a", "b"], n_epo).tolist() + else: + # For Epochs(TFR)Array, each epoch is an observation and they're already + # noisy/non-identical, so no duplication / noise-addition necessary. + insts = [cond_a, cond_b] + conds = ["a", "b"] + + # run new clustering API + df = pd.DataFrame(dict(data=insts, condition=conds)) + kwargs = dict(n_permutations=100, rng=42, tail=1, buffer_size=None, out_type="mask") + result_new_api = cluster_test(df, "data~condition", **kwargs) + # make sure channels are last dimension for old API + if is_epo: + assert result_new_api.stat_obs.shape == df["data"][0].get_data()[0, ...].shape + axes = (0, 3, 2, 1) if is_tfr else (0, 2, 1) + X = [cond_a.get_data().transpose(*axes), cond_b.get_data().transpose(*axes)] + else: + assert result_new_api.stat_obs.shape == df["data"][0].get_data().shape + axes = (2, 1, 0) if is_tfr else (1, 0) + Xa = list() + Xb = list() + for inst, cond in zip(insts, conds): + container = Xa if cond == "a" else Xb + container.append(inst.get_data().transpose(*axes)) + X = [np.stack(Xa), np.stack(Xb)] + + F_obs, clusters, cluster_pvals, H0 = permutation_cluster_test(X, **kwargs) + + for clust in result_new_api.clusters: + assert clust.shape == result_new_api.stat_obs.shape + + assert_array_almost_equal(result_new_api.H0, H0) + assert_array_almost_equal(result_new_api.stat_obs, F_obs.T) + assert_array_almost_equal(result_new_api.cluster_p_values, cluster_pvals) + assert len(result_new_api.clusters) == len(clusters) + for clu1, clu2 in zip(result_new_api.clusters, clusters): + assert_array_equal(clu1, clu2.T) + + +@pytest.mark.filterwarnings('ignore:Ignoring argument "tail":RuntimeWarning') +def test_cluster_test_rm_anova(): + """Test the interaction-formula (repeated-measures ANOVA) branch of cluster_test.""" + rng = np.random.default_rng(seed=0) + n_subjects, n_channels, n_times = 8, 3, 6 + info = create_info(n_channels, sfreq=100.0, ch_types="eeg") + factor_levels = [2, 2] + conditions = ["a1b1", "a1b2", "a2b1", "a2b2"] + data = { + cond: rng.normal(size=(n_subjects, n_channels, n_times)) for cond in conditions + } + # inject an interaction effect (crossover pattern) in the first 2 channels + data["a1b1"][:, :2] += 3 + data["a2b2"][:, :2] += 3 + data["a1b2"][:, :2] -= 3 + data["a2b1"][:, :2] -= 3 + + # reference: old-style call with a hand-rolled f_mway_rm stat_fun, exactly as + # done in tutorials/stats-sensor-space/70_cluster_rmANOVA_time_freq.py + def stat_fun(*args): + return f_mway_rm( + np.swapaxes(np.asarray(args), 1, 0), + factor_levels=factor_levels, + effects="A:B", + return_pvals=False, + )[0] + + f_thresh = f_threshold_mway_rm( + n_subjects, factor_levels, effects="A:B", pvalue=0.001 + ) + # channels last, as required by permutation_cluster_test + X_old = [data[cond].transpose(0, 2, 1) for cond in conditions] + kwargs = dict( + n_permutations=100, + tail=1, + rng=3, + buffer_size=None, + out_type="indices", + threshold=f_thresh, + ) + F_obs, clusters, cluster_pvals, H0 = permutation_cluster_test( + X_old, stat_fun=stat_fun, **kwargs + ) + + # new API: one row per (subject, condition), with an EvokedArray holding that + # subject's data for that condition + rows = list() + for cond in conditions: + for subj in range(n_subjects): + rows.append( + dict( + data=EvokedArray(data[cond][subj], info, tmin=0), + modality=cond[1], + location=cond[3], + subject=subj, + ) + ) + df = pd.DataFrame(rows) + result = cluster_test(df, "data ~ modality:location", within_id="subject", **kwargs) + + assert result.stat_name == "F-statistic (repeated-measures ANOVA)" + assert_array_almost_equal(result.stat_obs, F_obs.T) + assert len(result.clusters) == len(clusters) + for clu1, clu2 in zip(result.clusters, clusters): + assert_array_equal(clu1, tuple(reversed(clu2))) + # the observed stat and clusters match the legacy API, but the null differs + # by design: cluster_test permutes repeated-measures data within subject + # only, whereas the legacy API shuffles rows across the whole design + assert result.H0.shape == H0.shape + assert not np.allclose(result.H0, H0) + assert result.cluster_p_values.shape == cluster_pvals.shape + + +def test_cluster_test_formula_validation(stat_conditions): + """Test that cluster_test raises clear errors for unsupported formulas.""" + condition1_1d, condition2_1d, _, _ = stat_conditions + df = pd.DataFrame(dict(data=[condition1_1d, condition2_1d], a=["x", "y"])) + df["b"] = "z" + + # multi-term right-hand side ("a+b") is not a single effect + with pytest.raises(ValueError, match="single term"): + cluster_test(df, "data ~ a+b") + + # interaction effect requires within_id + with pytest.raises(ValueError, match="repeated-measures"): + cluster_test(df, "data ~ a:b") + + # unbalanced repeated-measures design (subject missing an observation) + rows = [ + dict(data=condition1_1d, a="x", b="p", subject=0), + dict(data=condition2_1d, a="x", b="q", subject=0), + dict(data=condition1_1d, a="y", b="p", subject=0), + # subject 0 is missing the "y"/"q" combination + dict(data=condition1_1d, a="x", b="p", subject=1), + dict(data=condition2_1d, a="x", b="q", subject=1), + dict(data=condition1_1d, a="y", b="p", subject=1), + dict(data=condition2_1d, a="y", b="q", subject=1), + ] + df_unbalanced = pd.DataFrame(rows) + with pytest.raises(ValueError, match="must have exactly"): + cluster_test(df_unbalanced, "data ~ a:b", within_id="subject") + + +@pytest.mark.filterwarnings('ignore:Ignoring argument "tail":RuntimeWarning') +@pytest.mark.filterwarnings("ignore:divide by zero:RuntimeWarning") +@pytest.mark.filterwarnings("ignore:invalid value encountered:RuntimeWarning") +@pytest.mark.filterwarnings("ignore:No clusters found:RuntimeWarning") +def test_cluster_test_reduce(stat_conditions): + """Reduce multiple observations for paired t-test.""" + # TODO: parametrize this test for Epochs, AveragedTFR etc. + + condition1_1d, _, _, _ = stat_conditions + # For this test we need equal sized arrays + condition2_1d = condition1_1d.copy() + rng = np.random.default_rng(0) + rng.shuffle(condition2_1d) + + info = create_info( + ch_names=[f"ch_{ii}" for ii in range(condition1_1d.shape[0])], + sfreq=10, + ch_types="eeg", + ) + data = [EvokedArray(arr, info) for arr in [condition1_1d, condition2_1d]] + df = pd.DataFrame(dict(data=data, a=["x", "y"])) + df["b"] = 1 + + df_2 = df.copy() + df_2["b"] = 2 + df = pd.concat([df, df_2]) + del df_2 + + df["c"] = "foo" + + df_2 = df.copy() + df_2["c"] = "bar" + df = pd.concat([df, df_2]) + del df_2 + # This should not raise + cluster_test(df, formula="data ~ a", within_id="c") diff --git a/mne/utils/docs.py b/mne/utils/docs.py index 92182d13885..328c12db657 100644 --- a/mne/utils/docs.py +++ b/mne/utils/docs.py @@ -140,61 +140,54 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): formatting. This can add overhead so is meant only for debugging. """ -docdict["adjacency_clust"] = """ -adjacency : scipy.sparse.spmatrix | None | False +_adjacency_clust_template = """ +adjacency : scipy.sparse.spmatrix | {param_none}False Defines adjacency between locations in the data, where "locations" can be spatial vertices, frequency bins, time points, etc. For spatial vertices (i.e. sensor space data), see :func:`mne.channels.find_ch_adjacency` or :func:`mne.spatial_inter_hemi_adjacency`. For source space data, see - :func:`mne.spatial_src_adjacency` or - :func:`mne.spatio_temporal_src_adjacency`. If ``False``, assumes - no adjacency (each location is treated as independent and unconnected). - If ``None``, a regular lattice adjacency is assumed, connecting - each {sp} location to its neighbor(s) along the last dimension - of {{eachgrp}} ``{{x}}``{lastdim}. + :func:`mne.spatial_src_adjacency` or :func:`mne.spatio_temporal_src_adjacency`. + If ``False``, assumes no adjacency (each location is treated as independent and + unconnected).{if_none} If ``adjacency`` is a matrix, it is assumed to be symmetric (only the upper triangular half is used) and must be square with dimension equal to - ``{{x}}.shape[-1]`` {parone} or ``{{x}}.shape[-1] * {{x}}.shape[-2]`` - {partwo} or (optionally) - ``{{x}}.shape[-1] * {{x}}.shape[-2] * {{x}}.shape[-3]`` - {parthree}.{memory} + the product of the last 1, 2, or 3 data dimensions (e.g., for time-frequency data: + n_channels, n_channels * n_freqs, or n_channels * n_freqs * n_times).{memory} +""" +_if_none = """ If ``None``, a regular lattice adjacency is assumed, connecting + each {spatial}location to its neighbor(s) along the last dimension + of {the_data}. """ - -mem = ( - " If spatial adjacency is uniform in time, it is recommended to use " - "a square matrix with dimension ``{x}.shape[-1]`` (n_vertices) to save " - "memory and computation, and to use ``max_step`` to define the extent " - "of temporal adjacency to consider when clustering." -) -comb = " The function `mne.stats.combine_adjacency` may be useful for 4D data." st = dict( - sp="spatial", - lastdim="", - parone="(n_vertices)", - partwo="(n_times * n_vertices)", - parthree="(n_times * n_freqs * n_vertices)", - memory=mem, + param_none="None | ", + if_none=_if_none.format(spatial="spatial ", the_data="{eachgrp} ``{x}``"), + memory=""" + If spatial adjacency is uniform in time, it is recommended to use a square matrix + with dimension ``{x}.shape[-1]`` (n_vertices) to save memory and computation, + and to use ``max_step`` to define the extent of temporal adjacency to consider when + clustering. +""", ) tf = dict( - sp="", - lastdim=" (or the last two dimensions if ``{x}`` is 2D)", - parone="(for 2D data)", - partwo="(for 3D data)", - parthree="(for 4D data)", - memory=comb, + param_none="None | ", + if_none=_if_none.format( + spatial="", + the_data="{eachgrp} ``{x}`` (or the last two dimensions if ``{x}`` is 2D)", + ), + memory=""" + The function `mne.stats.combine_adjacency` may be useful for 4D data. +""", ) -nogroups = dict(eachgrp="", x="X") +nogrps = dict(eachgrp="", x="X") groups = dict(eachgrp="each group ", x="X[k]") -docdict["adjacency_clust_1"] = ( - docdict["adjacency_clust"].format(**tf).format(**nogroups) -) -docdict["adjacency_clust_n"] = docdict["adjacency_clust"].format(**tf).format(**groups) -docdict["adjacency_clust_st1"] = ( - docdict["adjacency_clust"].format(**st).format(**nogroups) -) -docdict["adjacency_clust_stn"] = ( - docdict["adjacency_clust"].format(**st).format(**groups) + +docdict["adjacency_clust_1"] = _adjacency_clust_template.format(**tf).format(**nogrps) +docdict["adjacency_clust_both"] = _adjacency_clust_template.format( + param_none="", if_none="", memory="" ) +docdict["adjacency_clust_n"] = _adjacency_clust_template.format(**tf).format(**groups) +docdict["adjacency_clust_st1"] = _adjacency_clust_template.format(**st).format(**nogrps) +docdict["adjacency_clust_stn"] = _adjacency_clust_template.format(**st).format(**groups) docdict["adjust_dig_chpi"] = """ adjust_dig : bool @@ -785,7 +778,7 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): docdict["check_disjoint_clust"] = """ check_disjoint : bool - Whether to check if the connectivity matrix can be separated into disjoint + Whether to check if the ``adjacency`` matrix can be separated into disjoint sets before clustering. This may lead to faster clustering, especially if the second dimension of ``X`` (usually the "time" dimension) is large. """ @@ -1526,7 +1519,7 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): """ docdict["exclude_clust"] = """ -exclude : bool array or None +exclude : array-like of bool | None Mask to apply to the data to exclude certain points from clustering (e.g., medial wall vertices). Should be the same shape as ``X``. If ``None``, no points are excluded. @@ -4551,16 +4544,23 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): channel names in the file will be used when possible. """ -_stat_fun_clust_base = """ +_stat_fun_template = """ stat_fun : callable | None Function called to calculate the test statistic. Must accept 1D-array as - input and return a 1D array. If ``None`` (the default), uses - `mne.stats.{}`. + input and return a 1D array. If ``None`` (the default), uses {}. """ -docdict["stat_fun_clust_f"] = _stat_fun_clust_base.format("f_oneway") +docdict["stat_fun_clust_both"] = _stat_fun_template.format( + """:func:`mne.stats.ttest_1samp_no_p` + for paired tests and :func:`mne.stats.f_oneway` for unpaired tests or tests of + more than 2 groups.""" +) + +docdict["stat_fun_clust_f"] = _stat_fun_template.format(":func:`mne.stats.f_oneway`") -docdict["stat_fun_clust_t"] = _stat_fun_clust_base.format("ttest_1samp_no_p") +docdict["stat_fun_clust_t"] = _stat_fun_template.format( + ":func:`mne.stats.ttest_1samp_no_p`" +) docdict["static"] = """ static : instance of SpatialImage @@ -4773,10 +4773,10 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): threshold : float | dict | None The so-called "cluster forming threshold" in the form of a test statistic (note: this is not an alpha level / "p-value"). - If numeric, vertices with data values more extreme than ``threshold`` will - be used to form clusters. If ``None``, {} will be chosen + If numeric, vertices with stat values more extreme than ``threshold`` will + be used to form clusters. If ``None``, {which_thresh} will be chosen automatically that corresponds to a p-value of 0.05 for the given number of - observations (only valid when using {}). If ``threshold`` is a + observations (only valid when using {which_stat}). If ``threshold`` is a :class:`dict` (with keys ``'start'`` and ``'step'``) then threshold-free cluster enhancement (TFCE) will be used (see the :ref:`TFCE example ` and :footcite:`SmithNichols2009`). @@ -4784,8 +4784,14 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): a particular p-value for one-tailed or two-tailed tests. """ -f_test = ("an F-threshold", "an F-statistic") -docdict["threshold_clust_f"] = _threshold_clust_base.format(*f_test) +docdict["threshold_clust_both"] = _threshold_clust_base.format( + which_thresh="a t- or F-threshold", + which_stat="``stat_fun=None``, i.e., a paired t-test or one-way F-test", +) + +docdict["threshold_clust_f"] = _threshold_clust_base.format( + which_thresh="an F-threshold", which_stat="an F-statistic" +) docdict["threshold_clust_f_notes"] = """ For computing a ``threshold`` based on a p-value, use the conversion @@ -4797,8 +4803,9 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): thresh = scipy.stats.f.ppf(1 - pval, dfn=dfn, dfd=dfd) # F distribution """ -t_test = ("a t-threshold", "a t-statistic") -docdict["threshold_clust_t"] = _threshold_clust_base.format(*t_test) +docdict["threshold_clust_t"] = _threshold_clust_base.format( + which_thresh="a t-threshold", which_stat="a t-statistic" +) docdict["threshold_clust_t_notes"] = """ For computing a ``threshold`` based on a p-value, use the conversion @@ -4812,6 +4819,23 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): For testing the lower tail (``tail=-1``), don't subtract ``pval`` from 1. """ +docdict["threshold_clust_t_or_f_notes"] = """ +For computing a ``threshold`` based on a p-value, use the conversion +from :meth:`scipy.stats.rv_continuous.ppf`:: + + pval = 0.001 # arbitrary + # for t-statistic + df = n_observations - 1 # degrees of freedom for the t-test + thresh = scipy.stats.t.ppf(1 - pval / 2, df) # two-tailed, t distribution + # for f-statistic + dfn = n_conditions - 1 # degrees of freedom numerator + dfd = n_observations - n_conditions # degrees of freedom denominator + thresh = scipy.stats.f.ppf(1 - pval, dfn=dfn, dfd=dfd) # F distribution + +For a one-tailed test (``tail=1``), don't divide the p-value by 2. +For testing the lower tail (``tail=-1``), don't subtract ``pval`` from 1. +""" + docdict["time_bandwidth_tfr"] = """ time_bandwidth : float ``≥ 2.0`` Product between the temporal window length (in seconds) and the *full* diff --git a/pyproject.toml b/pyproject.toml index 2ea1ff516dc..0303ddf8b57 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -171,6 +171,7 @@ full-no-qt = [ "edfio >= 0.4.10", "eeglabio", "filelock >= 3.18.0", + "formulaic", "h5py >= 2.4", "imageio >= 2.6.1", "imageio-ffmpeg >= 0.4.1", diff --git a/tools/vulture_allowlist.py b/tools/vulture_allowlist.py index 0ded28e9a33..0f678ba4f28 100644 --- a/tools/vulture_allowlist.py +++ b/tools/vulture_allowlist.py @@ -177,6 +177,8 @@ _qt_raise_window _qt_disable_paint _qt_get_stylesheet + +# used in tutorial, not sure why shows up _show_help_fig # Called by Qt, or only from a subprocess (mne/viz/backends/tests/test_utils.py) diff --git a/tutorials/stats-sensor-space/76_new_cluster_test_api.py b/tutorials/stats-sensor-space/76_new_cluster_test_api.py new file mode 100644 index 00000000000..1c30c579a2a --- /dev/null +++ b/tutorials/stats-sensor-space/76_new_cluster_test_api.py @@ -0,0 +1,221 @@ +""" +.. _tut-new-cluster-test-api: + +=============================================================== +Group-level cluster permutation testing with formula contrasts +=============================================================== + +Run a cluster-based permutation test on evoked data from several subjects, +specifying the contrast with a Wilkinson (R-style) formula. By the end you +will have run a paired *t*-test across subjects and inspected the cluster +permutation results. + +You will: + +- load evoked data from multiple subjects +- build a long-format dataframe with one row per subject and condition +- run :func:`mne.stats.cluster_test` with a Wilkinson-notation formula +- inspect the cluster permutation results +""" +# Author: Carina Forster +# +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +# %% +# Load the required packages +# -------------------------- + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +import mne + +# %% +# Load evoked data for multiple subjects +# -------------------------------------- +# +# We use the P3b data from the freely available ERP CORE dataset: a visual +# oddball task contrasting rare *target* stimuli with frequent *non-target* +# stimuli. Each subject has an evoked response for both conditions; five +# subjects are included here. + +# path to the P3 dataset +path_to_p3 = mne.datasets.misc.data_path() / "ERP_CORE" / "P3" + +# participant IDs available in this dataset (15 to 19) +participant_ids = range(15, 20) + +# load each subject's evoked data into a list +evokeds_allsubs = [] +for pid in participant_ids: + # filename with the ID zero-padded to three digits + filename_p3 = f"sub-{pid:03d}_ses-P3_task-P3_ave.fif" + p3_file_path = Path(path_to_p3) / filename_p3 + evokeds = mne.read_evokeds(p3_file_path) + evokeds_allsubs.append(evokeds) + +# split the two conditions into separate per-subject lists +target_only = [evoked[0] for evoked in evokeds_allsubs] +non_target_only = [evoked[1] for evoked in evokeds_allsubs] + +# %% +# Inspect the contrast before testing +# ----------------------------------- +# +# Before running any statistics, look at the effect you are about to test. +# We form the per-subject difference (target minus non-target) and plot its +# grand average. A positive deflection means targets evoke the stronger +# response. + +diff_evoked = [ + mne.combine_evoked([evoked_target, evoked_non_target], weights=[1, -1]) + for evoked_target, evoked_non_target in zip(target_only, non_target_only) +] + +grand_avg_diff = mne.grand_average(diff_evoked) +grand_avg_diff.plot() +grand_avg_diff.plot_topomap() + +# %% +# You should see the largest difference around 400 ms over central-parietal +# channels -- the expected P3b effect, stronger for target stimuli. This is the +# contrast the cluster test will evaluate formally. +# +# ``diff_evoked`` is used only for this visualization. The cluster test below +# works from a dataframe holding *both* conditions and forms the contrast from +# the formula. + +# %% +# Build the dataframe for the cluster test +# ---------------------------------------- +# +# The formula interface takes a long-format :class:`pandas.DataFrame` with one +# row per observation. Each row holds one subject's evoked response for one +# condition, so every subject contributes two rows (target and non-target). +# Every subject must contribute the same set of conditions. The columns are: +# +# - ``evoked``: the single-subject :class:`~mne.Evoked` object +# - ``condition``: the condition label, referenced by the formula +# - ``subject_index``: identifies which observations are paired within a subject + +evokeds_conditions = target_only + non_target_only +conditions = ["target"] * len(target_only) + ["non-target"] * len(non_target_only) +subject_index = list(participant_ids) * 2 + +df = pd.DataFrame( + { + "evoked": evokeds_conditions, + "condition": conditions, + "subject_index": subject_index, + } +) +df + +# %% +# You should see the largest difference around 400 ms over central-parietal +# channels -- the expected P3b effect, stronger for target stimuli. This is the +# contrast the cluster test will evaluate formally. +# +# ``diff_evoked`` is used only for this visualization. The cluster test below +# works from a dataframe holding *both* conditions and forms the contrast from +# the formula. + +# %% +# The sign of the contrast follows the order of the condition levels. We set +# "target" as the first level so the difference is formed as target minus +# non-target (positive = stronger response to targets), matching the grand +# average we plotted above. + +# TODO: do this within cluster test? +df["condition"] = pd.Categorical( + df["condition"], categories=["target", "non-target"], ordered=True +) + +df + +# %% +# Run the cluster test with a formula +# ----------------------------------- +# +# The contrast is written as a Wilkinson (R-style) formula, the same notation +# used by R's ``lmer``/``glmer``. Here ``"evoked ~ condition"`` models the +# evoked response as a function of condition: ``condition`` is categorical and +# is dummy-coded automatically, and an intercept is included implicitly. +# Passing ``within_id="subject_index"`` makes this a within-subject (paired) +# test: the two conditions are subtracted within each subject and the resulting +# differences are tested against zero (a one-sample t-test), with the null +# distribution built by sign-flipping those per-subject differences. Because we +# set ``target`` as the first condition level above, the difference is formed as +# target minus non-target. TODO: should be a parameter in cluster_test? + +formula = "evoked ~ condition" + +cluster_result = mne.stats.cluster_test( + df=df, formula=formula, within_id="subject_index" +) + +print(f"Smallest cluster p-value: {cluster_result.cluster_p_values.min():.4f}") + +# %% +# The smallest cluster p-value is about 0.06, so no cluster is significant at +# alpha = 0.05 -- and with five subjects none ever could be. Here is why. +# +# The null distribution is built by sign-flipping the five per-subject +# difference scores. There are ``2 ** 5 = 32`` ways to assign signs, but +# flipping every sign only mirrors the partition, so just ``2 ** (5 - 1) = 16`` +# are distinct; excluding the observed arrangement leaves +# ``2 ** (5 - 1) - 1 = 15`` permutations. Because the test is exact, all 15 are +# evaluated (you will see ``15/15`` in the progress log) rather than sampled at +# random. +# +# The finest p-value this can resolve is ``1 / (15 + 1) = 0.0625``, so even the +# most extreme possible cluster lands just above 0.05. The near-0.0625 result +# means the observed cluster *was* the most extreme one -- there is simply not +# enough data to reach significance. Detecting an effect here would need more +# subjects; that, not the specific p-value, is the takeaway. + +# %% +# Inspect the results +# ------------------- +# +# 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. + +print( + f"Number of permutations run: {cluster_result.n_permutations}" +) # TODO: fix this in separate PR + +# 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; ensure it is arranged as (channels, times) +stat_obs = cluster_result.stat_obs +if stat_obs.shape != (len(ch_names), len(times)): + stat_obs = stat_obs.T + +# 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, +) +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")