[7f9fb8]: / mne / stats / tests / test_multi_comp.py

Download this file

57 lines (44 with data), 2.1 kB

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_almost_equal, assert_array_equal
from scipy import stats
from mne.stats import bonferroni_correction, fdr_correction
def test_bonferroni_pval_clip():
"""Test that p-values are never exceed 1.0."""
p = (0.2, 0.9)
_, p_corrected = bonferroni_correction(p)
assert p_corrected.max() <= 1.0
def test_multi_pval_correction():
"""Test pval correction for multi comparison (FDR and Bonferroni)."""
rng = np.random.RandomState(0)
X = rng.randn(10, 1000, 10)
X[:, :50, 0] += 4.0 # 50 significant tests
alpha = 0.05
T, pval = stats.ttest_1samp(X, 0)
n_samples = X.shape[0]
n_tests = X.size / n_samples
thresh_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)
reject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha)
thresh_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)
assert pval_bonferroni.ndim == 2
assert reject_bonferroni.ndim == 2
assert_allclose(pval_bonferroni, (pval * 10000).clip(max=1))
reject_expected = pval_bonferroni < alpha
assert_array_equal(reject_bonferroni, reject_expected)
fwer = np.mean(reject_bonferroni)
assert_almost_equal(fwer, alpha, 1)
reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method="indep")
assert pval_fdr.ndim == 2
assert reject_fdr.ndim == 2
thresh_fdr = np.min(np.abs(T)[reject_fdr])
assert 0 <= (reject_fdr.sum() - 50) <= 50 * 1.05
assert thresh_uncorrected <= thresh_fdr <= thresh_bonferroni
pytest.raises(ValueError, fdr_correction, pval, alpha, method="blah")
assert np.all(fdr_correction(pval, alpha=0)[0] == 0)
reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method="negcorr")
thresh_fdr = np.min(np.abs(T)[reject_fdr])
assert 0 <= (reject_fdr.sum() - 50) <= 50 * 1.05
assert thresh_uncorrected <= thresh_fdr <= thresh_bonferroni