Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correct permutation t-tests #684

Draft
wants to merge 13 commits into
base: develop
Choose a base branch
from
20 changes: 11 additions & 9 deletions moabb/analysis/meta_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ def _pairedttest_exhaustive(data):
pvals: ndarray of shape (n_pipelines, n_pipelines)
array of pvalues
"""
out = np.ones((data.shape[1], data.shape[1]))
out = np.zeros((data.shape[1], data.shape[1]), dtype=np.int32)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Initialization of out should be similar between _pairedttest_exhaustive and _pairedttest_random:

Suggested change
out = np.zeros((data.shape[1], data.shape[1]), dtype=np.int32)
out = np.ones(data.shape[1:], dtype=np.int32)

Copy link
Collaborator Author

@gcattan gcattan Feb 11, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm ok. But then we will take into account the original statistic two times, as we have random >= true

true = data.sum(axis=0)
nperms = 2 ** data.shape[0]
for perm in itertools.product([-1, 1], repeat=data.shape[0]):
Expand All @@ -98,11 +98,10 @@ def _pairedttest_exhaustive(data):
# multiply permutation by subject dimension and sum over subjects
randperm = (data * perm[:, None, None]).sum(axis=0)
# compare to true difference (numpy autocasts bool to 0/1)
out += randperm > true
out = out / nperms
# control for cases where pval is 1
out[out == 1] = 1 - (1 / nperms)
return out
out += randperm >= true

out[out >= nperms] = nperms - 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could add another control to avoid a p-value equal to 0:

Suggested change
out[out >= nperms] = nperms - 1
out[out >= nperms] = nperms - 1
out[out == 0] = 1

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree this may deserve a comment in the code. If out is initialized to 1, or if it is initialized to 0 and we have random >= true, then this check should not be necessary.

return out / nperms


def _pairedttest_random(data, nperms):
Expand All @@ -121,16 +120,19 @@ def _pairedttest_random(data, nperms):
pvals: ndarray of shape (n_pipelines, n_pipelines)
array of pvalues
"""
out = np.ones((data.shape[1], data.shape[1]))
out = np.ones(
(data.shape[1], data.shape[1]), dtype=np.int32
) # we use ones() so that out is never 0
gcattan marked this conversation as resolved.
Show resolved Hide resolved
true = data.sum(axis=0)
for _ in range(nperms):
perm = np.random.randint(2, size=(data.shape[0],))
perm[perm == 0] = -1
# multiply permutation by subject dimension and sum over subjects
randperm = (data * perm[:, None, None]).sum(axis=0)
# compare to true difference (numpy autocasts bool to 0/1)
out += randperm > true
out[out == nperms] = nperms - 1
out += randperm >= true

out[out >= nperms] = nperms - 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
out[out >= nperms] = nperms - 1
out[out >= nperms] = nperms - 1
out[out == 0] = 1

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here also, I will add a comment.

return out / nperms


Expand Down
Loading