[ENH] Spatio-temporal cluster permutation testing: Define minimum "nodes" to define a cluster

[sappelhoff]

For spatio-temporal cluster permutation testing, it's important to consider how many "nodes" make up a cluster.

For example:

• In the time domain (considering 1 "channel"): Are two neighboring significant time points already a cluster? Or does it start with "3", "4", "N"?
• In the spatial domain: the same but for neighboring channels who have are significant at the same point in time

Currently I can only assume that in MNE-Python, any neighboring relationship (2 timepoints, or 2 channels) makes up a cluster.

It would be really nice to

1. document this
2. offer a parameter to customize this

In FieldTrip, this is done as so:

cfg.minnbchan        = 2;          % minimum number of neighborhood channels that is
                                   % required for a selected sample to be included
                                   % in the clustering algorithm (default=0).

https://www.fieldtriptoolbox.org/tutorial/cluster_permutation_timelock/

xref #4859

[larsoner]

Currently I can only assume that in MNE-Python, any neighboring relationship (2 timepoints, or 2 channels) makes up a cluster.

Actually a single time point and single spatial location will even be a "cluster". The only thing that matters is that it exceeds the threshold.

Do you need to set minimum spatio/temporal/spatiotemporal "extent"s to get reasonable results? For data that are sufficiently smooth I don't think it shouldn't matter much in practice because single-node "clusters" will have very small counts and stat values (which you get/use depends on t_power) so won't be the largest cluster(s) for each permutation used in the end for the maximal statistic.

[sappelhoff]

I think you are right with the assessment that "it won't matter" in most cases, and I think for my data presently it also doesn't matter too much.

However, I still remember this parameter from when I was using FieldTrip and would have liked to use it during exploration, as it also controls the number of "observed" clusters.

will have very small counts and stat values (which you get/use depends on t_power)

Now that we are talking about t_power: Is there actually a meaningful value for that param beyond 0 (for count of "extent") and 1 (for sum of cluster stats)?

[larsoner]

it also controls the number of "observed" clusters.

This is a pretty easy operation (almost/possibly a one-liner) to do on the observed clusters after the fact, though. So I'd rather not add a parameter for this, but instead maybe a line in a tutorial/example.

Now that we are talking about t_power: Is there actually a meaningful value for that param beyond 0 (for count of "extent") and 1 (for sum of cluster stats)?

I guess that continuous values between 0 and 1 give a corresponding weight to the count vs stat. Then values greater than 1 give more importance to larger stat values, which at least reminds me of the difference between minimizing squared error versus absolute error / L2 vs L1 norms. But I'd suggest looking at the original papers/refs to see what they say, this is just my intuition from thinking about it.

[sappelhoff]

This is a pretty easy operation (almost/possibly a one-liner) to do on the observed clusters after the fact, though. So I'd rather not add a parameter for this, but instead maybe a line in a tutorial/example.

Fair enough! But I also just thought of another way why such a parameter would be relevant -- I drew a quick sketch (you need to tilt your head to the right should, sorry 😉)

In the drawn scenario, this is 1 single big cluster that has a much larger extent simply because there are a few channel/time "bins" that extend like "fingers" from the "main cluster" ... this even has the potential to connect clusters that would actually be distinct.

With the proposed parameter, such "fingers" can be removed

I guess that continuous values between 0 and 1 give a corresponding weight to the count vs stat. Then values greater than 1 give more importance to larger stat values, which at least reminds me of the difference between minimizing squared error versus absolute error / L2 vs L1 norms. But I'd suggest looking at the original papers/refs to see what they say, this is just my intuition from thinking about it.

thanks for this intuition!

[hoechenberger]

I've run into this problem with a cluster permutation in the time domain, where like @larsoner said a single time point can become a cluster. I now simply filter those out manually according to a user-defined threshold.

[hoechenberger]

I've run into this problem with a cluster permutation in the time domain, where like @larsoner said a single time point can become a cluster. I now simply filter those out manually according to a user-defined threshold.

I'll make a PR for the study template soon that demonstrates this.

[larsoner]

In the drawn scenario, this is 1 single big cluster that has a much larger extent simply because there are a few channel/time "bins" that extend like "fingers" from the "main cluster" ... this even has the potential to connect clusters that would actually be distinct.

This is not the case for your example according to a standard lattice adjacency, there are four clusters:

>>> import numpy as np
>>> x = np.zeros((6, 5), bool)
>>> x[[0, 1, 2, 3, 3, 4, 4, 5], [0, 1, 4, 2, 3, 2, 3, 2]] = True
>>> 
>>> x
array([[ True, False, False, False, False],
       [False,  True, False, False, False],
       [False, False, False, False,  True],
       [False, False,  True,  True, False],
       [False, False,  True,  True, False],
       [False, False,  True, False, False]])
>>> from scipy import ndimage
>>> ndimage.label(x)
(array([[1, 0, 0, 0, 0],
       [0, 2, 0, 0, 0],
       [0, 0, 0, 0, 3],
       [0, 0, 4, 4, 0],
       [0, 0, 4, 4, 0],
       [0, 0, 4, 0, 0]], dtype=int32), 4)

This is because lattice adjacency (reasonably IMO) is a "cross" not a "filled 3x3 square" in the 2D case.

[sappelhoff]

This is not the case for your example according to a standard lattice adjacency, there are four clusters:

I see - but does that make my point invalid? 🙂

Because we could tune my example such that it would result in a single cluster, ... and still look like something that I'd rather control to not be a single cluster (by a parameter that tunes a "minimum extent").

Full disclosure: I won't have time to do this PR myself at this time. But I still don't think (yet?) that this issue needs to be closed.

[larsoner]

Because we could tune my example such that it would result in a single cluster, ... and still look like something that I'd rather control to not be a single cluster (by a parameter that tunes a "minimum extent").

Okay I'm not sure what you mean / how you want this "min" operation to operate. Can you provide an example of current behavior that would be bad, what you'd want the parameter to fix? For example, based on how I've been thinking about the problem, I would expect a min_spatial=2 or so to produce a cluster with 4 nodes for the thresholded map (time in rows, space in cols):

[[0, 1],
 [0, 1],
 [1, 1]]

But I'm starting to suspect you're saying that you'd want this to produce a cluster with just the last row, which is indeed different from how I've been thinking about the minimum extent... I'm thinking about a minimum after forming the cluster like an "any" operation, and you're thinking about a minimum during clustering like an "all" operation (i.e., "it has to maintain a minimum extent of X along dimension Y at all locations, rather than achieve it at any point")...?

Which way does the FieldTrip parameter actually operate?

[sappelhoff]

The fieldtrip parameter description says:

cfg.minnbchan        = 2;          % minimum number of neighborhood channels that is
                                   % required for a selected sample to be included
                                   % in the clustering algorithm (default=0).

So I think the default is identical to what we currently have in MNE.

with cfg.minnbchan=1 I would expect that for:

[[0, 1],
 [0, 1],
 [1, 1]]

... the first column would contain just a single channel significant at this timepoint 0, i.e., there are 0 neighbors -- so this does not go into the cluster, and -- as you say -- only the last column will be part of the cluster.

I think this can be helpful if one wants to avoid cases like this (just conceptually):

[[1, 0, 0, 1, 1],
 [1, 0, 0, 1, 1],
 [1, 1, 1, 1, 1]]

where 2 (or more) clusters are "bridged" by small connections (i.e., channels that are more sustainedly significant, when the rest of the cluster is already n.s. again).

[larsoner]

Okay that makes sense. It seems doable in principle given our existing code. I think we can achieve minimums for:

1. Time points spanned by a single spatial node by making a preliminary pass over the columns (time series) of thresholded array to just zero out any contiguous series of 1's that is not long enough.
2. Spatial nodes by working within the clustering step to just drop any "cluster" within a time slice that is not big enough.

Just to clarify, what should happen in this case for min_t_span=2, min_spatial_span=2?

[[0, 1, 0],
 [1, 1, 1],
 [0, 1, 0]]

Should there even be 1 cluster? I think no, since the singletons should be culled by the algorithm...

As far as the naming goes, we could consider min_t_span=1(default) rather than min_t_neighbors=0. I at least find thinking about the span more natural than the number of neighbors (even though they just differ by 1).

[sappelhoff]

Should there even be 1 cluster? I think no, since the singletons should be culled by the algorithm...

agreed

As far as the naming goes, we could consider min_t_span=1(default) rather than min_t_neighbors=0. I at least find thinking about the span more natural than the number of neighbors (even though they just differ by 1).

I think that most will agree with this as well 👍

[larsoner]

Okay then we just need a volunteer :)

[mmagnuski]

I agree that it is useful to have such parameters when doing cluster based permutation test. minnbchan can be very useful in avoiding birdges between clusters (ie two clusters merged into one because there is a relatively narrow path joining them).

I actually have some code for things like that in my borsar repository (which unfortunatelly, despite my best intentions, does not contain much good documentation apart from docstrings; but has good test coverage at least :) ). Unfortunatelly it does clustering differently to mne - uses separate functions (both numba and numpy) specialized for common scenarios (like 3d space with channel adjacency (TFR), 2d space with adjacency (ERP) etc.) - so making this code part of mne might not be straightforward.
If someone would like to try it out, here is a short example of finding clusters using minimum number of adjacent channels (min_adj_ch below):

import numpy as np
import borsar


data = np.array([[1, 0, 0, 1, 1],
                 [1, 0, 0, 1, 1],
                 [1, 1, 1, 1, 1]])
adjacency = np.array([[False, True,  False],
                      [True,  False, True],
                      [False, True,  False]])

# %% find clusters
clusters, cluster_stats = borsar.cluster.find_clusters(
    data.copy(), threshold=0.5, adjacency=adjacency, min_adj_ch=1)

this returns two boolean 2d arrays in clusters because there are two clusters when eliminating samples without at least one adjacent supra-threshold sample in the channels dimension. There are also higher level functions for performing the whole cluster based permutation test.