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
8 changes: 4 additions & 4 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]))
true = data.sum(axis=0)
nperms = 2 ** data.shape[0]
for perm in itertools.product([-1, 1], repeat=data.shape[0]):
Expand All @@ -98,7 +98,7 @@ 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 += randperm >= true
out = out / nperms
# control for cases where pval is 1
out[out == 1] = 1 - (1 / nperms)
Expand All @@ -121,15 +121,15 @@ 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.zeros((data.shape[1], data.shape[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.

In the case of exhaustive, this is not necessary.
But in the case of random, the initial statistic is not necessarily among those found within the permutation distribution.

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 += randperm >= true
out[out == nperms] = nperms - 1
return out / nperms
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Should it be out / (nperms + 1) instead?
So the correction apply to all p-value? And not to the extreme one only.


Expand Down
Loading