Correct permutation t-tests #684

gcattan · 2025-02-07T17:19:19Z

I think the out array should be initiated with zeros.
If not, then out can take nperms+1 as values, and the corrective statements (out[out == nperms] = nperms - 1) apply on the wrong sample.

@bruAristimunha what do you think?

Signed-off-by: gcattan <[email protected]>

bruAristimunha · 2025-02-07T17:46:55Z

hey @gcattan,

I am not a super expert on t-tests, I am pinging @sylvchev here.

gcattan · 2025-02-07T19:35:23Z

Ok thanks ! Just to give more context: we wanted to display the summary plot with many datasets. When computing the dataset statistics, we spot this warning saying there were some NaN values out of the stouffer method. When looking at the p-values some were greater or equals to 1. Envoyé à partir de Outlook pour Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Bru ***@***.***> Sent: Friday, February 7, 2025 6:47:18 PM To: NeuroTechX/moabb ***@***.***> Cc: gcattan ***@***.***>; Mention ***@***.***> Subject: Re: [NeuroTechX/moabb] Correct permutation t-tests (PR #684) hey @gcattan<https://github.com/gcattan>, I am not a super expert on t-tests, I am pinging @sylvchev<https://github.com/sylvchev> here. — Reply to this email directly, view it on GitHub<#684 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABPQYJ2C6QG7QBJ7KBH45GL2OTWSNAVCNFSM6AAAAABWWNWHOGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNBTGU4TKNRQGQ>. You are receiving this because you were mentioned.Message ID: ***@***.***>

gcattan · 2025-02-10T13:17:14Z

moabb/analysis/meta_analysis.py

@@ -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]))


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.

gcattan · 2025-02-10T13:18:58Z

moabb/analysis/meta_analysis.py

    out[out == nperms] = nperms - 1
    return out / nperms


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

Switch back to 1 for the random t-test to have at least 1 in the result. And changed check for p=1 values for _pairedttest_exhaustive to be the same as for _pairedttest_random.

Added comment

toncho11 · 2025-02-10T15:17:25Z

This last version was consulted with @Marco-Congedo.

qbarthelemy

Unit tests could be completed to test that p-values are never equal to 0.

qbarthelemy · 2025-02-11T10:49:54Z

moabb/analysis/meta_analysis.py

-    return out
+        out += randperm >= true
+
+    out[out >= nperms] = nperms - 1


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

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.

qbarthelemy · 2025-02-11T10:49:58Z

moabb/analysis/meta_analysis.py

-    out[out == nperms] = nperms - 1
+        out += randperm >= true
+
+    out[out >= nperms] = nperms - 1


Suggested change

out[out >= nperms] = nperms - 1

out[out >= nperms] = nperms - 1

out[out == 0] = 1

Here also, I will add a comment.

qbarthelemy · 2025-02-11T10:57:24Z

moabb/analysis/meta_analysis.py

@@ -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)


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)

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

moabb/analysis/meta_analysis.py

toncho11 · 2025-02-11T12:10:27Z

The main problem was that there were p values produced equal to 1 or bigger than 1.

Co-authored-by: Quentin Barthélemy <[email protected]> Signed-off-by: gcattan <[email protected]>

moabb/analysis/meta_analysis.py

Co-authored-by: Quentin Barthélemy <[email protected]> Signed-off-by: gcattan <[email protected]>

PierreGtch · 2025-02-11T16:25:31Z

Hi @gcattan, thank you for investing time into this issue!

As we can see here, custom re-implementations in our codebase are prone to errors. Do you think it would be possible to reimplement these using functions from common libraries, such as scipy.stats.permutation_test, which have been validated by a larger community?

If yes, I think it could be done in two steps:

Test if the current implementation or your fix is consistent with the SciPy version.
Replace the custom implementation with the SciPy version.

These two steps could be done in two PRs eventually.

Marco-Congedo · 2025-02-11T16:43:43Z

Hi,

the problem is not with the implementation, but with the fact that a permutation test can give a p-value in [1/nperm, 1] and that any p-value equal to 1 will screw up the Liptak (Stouffer) combination function as z(1)=Inf. The proposed solution is to allow p to be in [1/nperm, 1-1/nperm] only. This is not very straight, but i don't see any other solution if the Liptak combination function is to be used. In any case, using any other implementation, one would still need to constraint the p-values in [1/nperm, 1-1/nperm].

gcattan · 2025-02-12T10:23:11Z

@PierreGtch Then the question is rather if you want to change the current implementation of combine_pvalues:

What I can do as part of this PR is to validate the latest modification against the datasets we used and push some unit testing as suggested by @qbarthelemy.

bruAristimunha · 2025-02-12T10:38:01Z

I think it's a great way to go, the less we keep the better.

…

On Wed, 12 Feb 2025, 11:23 gcattan, ***@***.***> wrote: @PierreGtch <https://github.com/PierreGtch> Then the question is rather if you want to change the current implementation of combine_pvalues: image.png (view on web) <https://github.com/user-attachments/assets/cb2822e6-512b-4fa4-a8d6-6a6c558ef642> What I can do as part of this PR is to validate the latest modification against the datasets we used and push some unit testing as suggested by @qbarthelemy <https://github.com/qbarthelemy>. — Reply to this email directly, view it on GitHub <#684 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKFZNATSJ56423F5BEG76IL2PMOKTAVCNFSM6AAAAABWWNWHOGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJTGI4DSMZUGQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

bruAristimunha · 2025-02-12T10:38:43Z

Regarding the discussion about 1 and 0, I'll let Sylvain have the final word and merge On Wed, 12 Feb 2025, 11:37 Bruno Aristimunha Pinto, ***@***.***> wrote:

…

I think it's a great way to go, the less we keep the better. On Wed, 12 Feb 2025, 11:23 gcattan, ***@***.***> wrote: > @PierreGtch <https://github.com/PierreGtch> Then the question is rather > if you want to change the current implementation of combine_pvalues: > > image.png (view on web) > <https://github.com/user-attachments/assets/cb2822e6-512b-4fa4-a8d6-6a6c558ef642> > > What I can do as part of this PR is to validate the latest modification > against the datasets we used and push some unit testing as suggested by > @qbarthelemy <https://github.com/qbarthelemy>. > > — > Reply to this email directly, view it on GitHub > <#684 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AKFZNATSJ56423F5BEG76IL2PMOKTAVCNFSM6AAAAABWWNWHOGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJTGI4DSMZUGQ> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

gcattan · 2025-02-12T10:55:30Z

Sure, but let's wait for testing before merging, this is not covered by the current test suite. I will try to tackle this during the week.

PierreGtch · 2025-02-12T11:29:15Z

Hi,

the problem is not with the implementation, but with the fact that a permutation test can give a p-value in [1/nperm, 1] and that any p-value equal to 1 will screw up the Liptak (Stouffer) combination function as z(1)=Inf. The proposed solution is to allow p to be in [1/nperm, 1-1/nperm] only. This is not very straight, but i don't see any other solution if the Liptak combination function is to be used. In any case, using any other implementation, one would still need to constraint the p-values in [1/nperm, 1-1/nperm].

I never used this code from moabb, so your guess is better than mine :)
My remark was more general : it is a better habit to use code that is already validated than re-writing custom implementations

gcattan · 2025-02-12T11:49:47Z

The implementation is custom because you need to correct your p-value for the Stouffer method of combining p-value. If you really want to move away from custom implementation to a more standardized one, you need to use another method of combining p-value (in MoABB). Marco can probably advise what it is best to do in this situation. The trade off is that previously reported analysis may have slightly different result. Envoyé à partir de Outlook pour Android<https://aka.ms/AAb9ysg>

…

________________________________ From: Pierre Guetschel ***@***.***> Sent: Wednesday, February 12, 2025 12:29:39 PM To: NeuroTechX/moabb ***@***.***> Cc: gcattan ***@***.***>; Mention ***@***.***> Subject: Re: [NeuroTechX/moabb] Correct permutation t-tests (PR #684) Hi, the problem is not with the implementation, but with the fact that a permutation test can give a p-value in [1/nperm, 1] and that any p-value equal to 1 will screw up the Liptak (Stouffer) combination function as z(1)=Inf. The proposed solution is to allow p to be in [1/nperm, 1-1/nperm] only. This is not very straight, but i don't see any other solution if the Liptak combination function is to be used. In any case, using any other implementation, one would still need to constraint the p-values in [1/nperm, 1-1/nperm]. I never used this code from moabb, so your guess is better than mine :) My remark was more general : it is a better habit to use code that is already validated than re-writing custom implementations — Reply to this email directly, view it on GitHub<#684 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABPQYJZV5KNNHJ53XFW3POL2PMWCHAVCNFSM6AAAAABWWNWHOGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJTGQZTSNJQHA>. You are receiving this because you were mentioned.Message ID: ***@***.***>

Marco-Congedo · 2025-02-12T12:48:17Z

The implementation is custom because you need to correct your p-value for the Stouffer method of combining p-value. If you really want to move away from custom implementation to a more standardized one, you need to use another method of combining p-value (in MoABB). Marco can probably advise what it is best to do in this situation. The trade off is that previously reported analysis may have slightly different result. Envoyé à partir de Outlook pour Androidhttps://aka.ms/AAb9ysg
…
________________________________ From: Pierre Guetschel @.> Sent: Wednesday, February 12, 2025 12:29:39 PM To: NeuroTechX/moabb @.> Cc: gcattan @.>; Mention @.> Subject: Re: [NeuroTechX/moabb] Correct permutation t-tests (PR #684) Hi, the problem is not with the implementation, but with the fact that a permutation test can give a p-value in [1/nperm, 1] and that any p-value equal to 1 will screw up the Liptak (Stouffer) combination function as z(1)=Inf. The proposed solution is to allow p to be in [1/nperm, 1-1/nperm] only. This is not very straight, but i don't see any other solution if the Liptak combination function is to be used. In any case, using any other implementation, one would still need to constraint the p-values in [1/nperm, 1-1/nperm]. I never used this code from moabb, so your guess is better than mine :) My remark was more general : it is a better habit to use code that is already validated than re-writing custom implementations — Reply to this email directly, view it on GitHub<#684 (comment)>, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABPQYJZV5KNNHJ53XFW3POL2PMWCHAVCNFSM6AAAAABWWNWHOGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJTGQZTSNJQHA. You are receiving this because you were mentioned.Message ID: @.***>

Yeah, and who knows how the problem has been handled before, as it is not uncommon that a pipeline gives worse results as compared to another for all subjects, which by definition results in a p-value of 1 when running a paired t-test. The Fisher combination function does not have this limitation (and probably that is why Fisher did not propose the Liptak combination function) since it is based on the sum of logarithms of the p-values; now, for a permutation test those are in [1/nperm, 1], thus the sum is always defined. However, using the Fisher combination function weights cannot be used anymore and if i understand well currently MOABB weights the combinations to account for the different numerosity in different databases. That is why, all in all, if you wish to stick to the current methodology for comparing pipelines across databases in MOABB, i don't see for the moment a better soluation than artificially restricting the p-values in [1/nperm, 1-1/nperm].

toncho11 · 2025-02-12T13:27:16Z

One important thing for the future test is that we had undeterministic error. That's it. For the same code and results data it will sometimes produce p values equal or bigger than 1. This is probably because of the random test. So it was actually hard to detect. So if you make a test ... maybe the test should be run 20 times in order to be sure that the handling of the p values is correct. Because even before there were guards for the p values, but they were not working. Or this was the idea - to warn you that these p values are strange.

Correct permutation t-tests

25d1022

Signed-off-by: gcattan <[email protected]>

gcattan and others added 2 commits February 10, 2025 14:09

Update meta_analysis.py

0667dc6

Merge branch 'develop' into patch-2

c7114de

gcattan commented Feb 10, 2025

View reviewed changes

toncho11 and others added 4 commits February 10, 2025 16:08

Update meta_analysis.py

4d546f7

Switch back to 1 for the random t-test to have at least 1 in the result. And changed check for p=1 values for _pairedttest_exhaustive to be the same as for _pairedttest_random.

[pre-commit.ci] auto fixes from pre-commit.com hooks

7a7c4f4

Update meta_analysis.py

a5e6160

Added comment

[pre-commit.ci] auto fixes from pre-commit.com hooks

9495ae0

qbarthelemy reviewed Feb 11, 2025

View reviewed changes

gcattan and others added 5 commits February 11, 2025 14:18

Apply suggestions from code review

507d028

Co-authored-by: Quentin Barthélemy <[email protected]> Signed-off-by: gcattan <[email protected]>

Update meta_analysis.py

14c04a6

[pre-commit.ci] auto fixes from pre-commit.com hooks

d199d2b

Merge branch 'develop' into patch-2

8b94e0d

Update meta_analysis.py

acb565a

qbarthelemy reviewed Feb 11, 2025

View reviewed changes

moabb/analysis/meta_analysis.py Outdated Show resolved Hide resolved

Update moabb/analysis/meta_analysis.py

d2518e2

Co-authored-by: Quentin Barthélemy <[email protected]> Signed-off-by: gcattan <[email protected]>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Correct permutation t-tests #684

Correct permutation t-tests #684

gcattan commented Feb 7, 2025 •

edited

Loading

bruAristimunha commented Feb 7, 2025

gcattan commented Feb 7, 2025 via email

gcattan Feb 10, 2025

gcattan Feb 10, 2025

toncho11 commented Feb 10, 2025

qbarthelemy left a comment

qbarthelemy Feb 11, 2025

gcattan Feb 11, 2025

qbarthelemy Feb 11, 2025

gcattan Feb 11, 2025

qbarthelemy Feb 11, 2025

gcattan Feb 11, 2025 •

edited

Loading

toncho11 commented Feb 11, 2025 •

edited

Loading

PierreGtch commented Feb 11, 2025

Marco-Congedo commented Feb 11, 2025

gcattan commented Feb 12, 2025

bruAristimunha commented Feb 12, 2025 via email

bruAristimunha commented Feb 12, 2025 via email

gcattan commented Feb 12, 2025 •

edited

Loading

PierreGtch commented Feb 12, 2025

gcattan commented Feb 12, 2025 via email

Marco-Congedo commented Feb 12, 2025

toncho11 commented Feb 12, 2025

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

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

Correct permutation t-tests #684

Are you sure you want to change the base?

Correct permutation t-tests #684

Conversation

gcattan commented Feb 7, 2025 • edited Loading

bruAristimunha commented Feb 7, 2025

gcattan commented Feb 7, 2025 via email

gcattan Feb 10, 2025

Choose a reason for hiding this comment

gcattan Feb 10, 2025

Choose a reason for hiding this comment

toncho11 commented Feb 10, 2025

qbarthelemy left a comment

Choose a reason for hiding this comment

qbarthelemy Feb 11, 2025

Choose a reason for hiding this comment

gcattan Feb 11, 2025

Choose a reason for hiding this comment

qbarthelemy Feb 11, 2025

Choose a reason for hiding this comment

gcattan Feb 11, 2025

Choose a reason for hiding this comment

qbarthelemy Feb 11, 2025

Choose a reason for hiding this comment

gcattan Feb 11, 2025 • edited Loading

Choose a reason for hiding this comment

toncho11 commented Feb 11, 2025 • edited Loading

PierreGtch commented Feb 11, 2025

Marco-Congedo commented Feb 11, 2025

gcattan commented Feb 12, 2025

bruAristimunha commented Feb 12, 2025 via email

bruAristimunha commented Feb 12, 2025 via email

gcattan commented Feb 12, 2025 • edited Loading

PierreGtch commented Feb 12, 2025

gcattan commented Feb 12, 2025 via email

Marco-Congedo commented Feb 12, 2025

toncho11 commented Feb 12, 2025

gcattan commented Feb 7, 2025 •

edited

Loading

gcattan Feb 11, 2025 •

edited

Loading

toncho11 commented Feb 11, 2025 •

edited

Loading

gcattan commented Feb 12, 2025 •

edited

Loading