TFCE Calculation Error

[charfisher]

Description of the problem

The outputs of mne.stats.permutation_cluster_1samp_test() do not report the correct clusters when TFCE is being used. The returned "t_obs" variable does not appear to be the t statistics, rather the TFCE values (although these are incorrect). Additionally, the reported number of clusters is equal to the number of time samples. I believe I have found part of the issue. However, since this function is being used with other functions, I do not want to push this to git in case it affects other related functions. All of these issues are in _find_clusters() within mne.stats.cluster_level.py

Issue 1: Extent of cluster is incorrectly calculated

Lines 472-482,

for c in clusters:
      if isinstance(c, slice):
          len_c = c.stop - c.start
      elif isinstance(c, tuple):
          len_c = len(c)
      elif c.dtype == np.dtype(bool):
          len_c = np.sum(c)
      else:
          len_c = len(c)
      scores[c] += h * (len_c**e_power)

Each c is a tuple, but the tuple contains a singular slice. Therefore, len_c is 1 regardless of the length of the cluster.

Issue 2: Height of cluster is incorrectly calculated

Lines 465-471,

# the score of each point is the sum of the h^H * e^E for each
# supporting section "rectangle" h x e.
if ti == 0:
    h = abs(thresh)
else:
    h = abs(thresh - thresholds[ti - 1])
h = h**h_power 

I believe this should be:

# the score of each point is the sum of the h^H * e^E for each
# supporting section "rectangle" h x e.
if ti == 0:
    dh = abs(thresh)
else:
    dh = abs(thresh - thresholds[ti - 1])
h = (abs(thresh)**h_power) * dh

Issue 3: Incorrect number of clusters reported

This issue still remains. I didn't look further into this.

Steps to reproduce

import numpy as np
import matplotlib.pyplot as plt
import scipy
import mne
from mne.io import concatenate_raws, read_raw_edf
from mne.datasets import eegbci

subjects = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
runs = [6, 10, 14] 
sfreq = 160
tmin = -1
tmax = 4

# gets an ERP for each subject from -1 to 4s 
erp = np.zeros((len(subjects), int((tmax-tmin)*sfreq) + 1))

for i in range(0, len(subjects)):
    raw_fnames = eegbci.load_data(subjects[i], runs)
    raws = [read_raw_edf(f, preload=True) for f in raw_fnames]
    # concatenate runs from subject
    raw = concatenate_raws(raws)

    raw.annotations.rename(dict(T1="hands", T2="feet"))  # as documented on PhysioNet
    raw.set_eeg_reference(projection=True)

    # Apply band-pass filter
    raw.filter(1., 10.0, fir_design="firwin", skip_by_annotation="edge")

    picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False, exclude="bads")
    
    # Read epochs (train will be done only between 1 and 2s)
    # Testing will be done with a running classifier
    epochs = mne.Epochs(
        raw,
        event_id=["hands", "feet"],
        tmin=tmin,
        tmax=tmax,
        proj=True,
        picks=picks,
        baseline=None,
        preload=True,
    )
    
    erp[i,:] = np.mean(epochs['hands'].get_data(), axis=(0,1)) # average over epochs and channels for "hands" trial

thresh = {'start': 0, 'step': 0.1}

t_obs, clusters, p_vals, h0 = mne.stats.permutation_cluster_1samp_test(erp, threshold=thresh, tail=0, adjacency=None, max_step=1)
time = np.linspace(tmin, tmax, int((tmax-tmin)*sfreq) + 1)

plt.figure()
plt.plot(time, erp.T)
plt.ylabel("Ampltiude")
plt.xlabel("Time (s)")
plt.title("ERP for each subject")

plt.figure()
plt.plot(time, scipy.stats.ttest_1samp(erp, popmean=0).statistic)
plt.ylabel("t-values from scipy")
plt.xlabel("Time (s)")

plt.figure()
plt.plot(time, t_obs)
plt.ylabel("'t-values' from MNE with step = 0.1")
plt.xlabel("Time (s)")

thresh = {'start': 0, 'step': 1}
t_obs, clusters, p_vals, h0 = mne.stats.permutation_cluster_1samp_test(erp, threshold=thresh, tail=0, adjacency=None, max_step=1)

plt.figure()
plt.plot(time, t_obs)
plt.ylabel("'t-values' from MNE with step = 0.5")
plt.xlabel("Time (s)")

Link to data

No response

Expected results

There should be a cluster around 500 ms. The TFCE values should look something like this:

Actual results

This is what the TFCE values (or "t_obs") looks like:

Additional information

Platform macOS-26.6.2-arm64-arm-64bit
Python 3.10.13 | packaged by conda-forge | (main, Dec 23 2023, 15:35:25) [Clang 16.0.6 ]
Executable /Users/charlief/miniconda3/bin/python
CPU Apple M3 Ultra (28 cores)
Memory 96.0 GiB

Core
├☑ mne 1.12.1 (latest release)
├☑ numpy 2.2.6 (unknown linalg bindings (threadpoolctl module not found: No module named 'threadpoolctl'))
├☑ scipy 1.15.3
└☑ matplotlib 3.10.9 (backend=module://matplotlib_inline.backend_inline)

Numerical (optional)
└☐ unavailable sklearn, numba, nibabel, nilearn, dipy, openmeeg, cupy, pandas, h5io, h5py

Visualization (optional)
├☑ qtpy 2.4.3 (PyQt5=5.15.15)
└☐ unavailable pyvista, pyvistaqt, vtk, ipympl, pyqtgraph, mne-qt-browser, ipywidgets, trame_client, trame_server, trame_vtk, trame_vuetify

Ecosystem (optional)
├☑ defusedxml 0.7.1
└☐ unavailable mne-bids, mne-nirs, mne-features, mne-connectivity, mne-icalabel, mne-bids-pipeline, neo, eeglabio, edfio, curryreader, mffpy, pybv, pymef, antio

[larsoner]

These do sound like genuine bugs, but I haven't verified them. Would you be up for attempting a fix? If so, I'd suggest starting from some synthetic ground truth data where the result is known and verifiable. Turn this into a little unit test. It should fail on main, and pass on the PR. Some version of the data from this tutorial might make sense as a good test case to start from.

[charfisher]

Sure! I can take a crack at it.

[scott-huberty]

Each c is a tuple, but the tuple contains a singular slice. Therefore, len_c is 1 regardless of the length of the cluster.

---

Just a hint, I don't think this is the bug. Per the code you referenced, when isinstance(c, slice) then len_c = c.stop - c.start, which wouldnt' be 1 unless c.stop is somehow guaranteed to be c.start + 1 (?)

[charfisher]

That is correct when c is a slice, but c is a tuple that contains a slice. Therefore, to get the correct length in an inelegant way, it would be len_c = c[0].stop - c[0].start.