[med-svn] [python-mne] 45/52: ENH : adding support for FDR and Bonferonni p-value correction
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:49 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to annotated tag v0.2
in repository python-mne.
commit c17966ab91e9ba162ef34a2857783d4f7dde7cc8
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Mon Nov 7 18:00:20 2011 -0500
ENH : adding support for FDR and Bonferonni p-value correction
---
examples/stats/plot_fdr_stats_evoked.py | 80 +++++++++++++++++++++++++++
mne/stats/__init__.py | 1 +
mne/stats/multi_comp.py | 98 +++++++++++++++++++++++++++++++++
mne/stats/tests/test_multi_comp.py | 36 ++++++++++++
4 files changed, 215 insertions(+)
diff --git a/examples/stats/plot_fdr_stats_evoked.py b/examples/stats/plot_fdr_stats_evoked.py
new file mode 100644
index 0000000..0ff6d11
--- /dev/null
+++ b/examples/stats/plot_fdr_stats_evoked.py
@@ -0,0 +1,80 @@
+"""
+=======================================
+FDR correction on T-test on sensor data
+=======================================
+
+One tests if the evoked response significantly deviates from 0.
+Multiple comparison problem is adressed with
+False Discovery Rate (FDR) correction.
+
+"""
+
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+from scipy import stats
+import mne
+from mne import fiff
+from mne.datasets import sample
+from mne.stats import bonferroni_correction, fdr_correction
+
+###############################################################################
+# Set parameters
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
+event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
+event_id, tmin, tmax = 1, -0.2, 0.5
+
+# Setup for reading the raw data
+raw = fiff.Raw(raw_fname)
+events = mne.read_events(event_fname)[:30]
+
+channel = 'MEG 1332' # include only this channel in analysis
+include = [channel]
+
+###############################################################################
+# Read epochs for the channel of interest
+picks = fiff.pick_types(raw.info, meg=False, eog=True, include=include)
+event_id = 1
+reject = dict(grad=4000e-13, eog=150e-6)
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+ baseline=(None, 0), reject=reject)
+X = epochs.get_data() # as 3D matrix
+X = X[:, 0, :] # take only one channel to get a 2D array
+
+###############################################################################
+# Compute statistic
+T, pval = stats.ttest_1samp(X, 0)
+alpha = 0.05
+
+n_samples, n_tests = X.shape
+threshold_uncorrected = stats.t.ppf(1.0 - alpha, n_samples - 1)
+
+reject_bonferroni, pval_bonferroni = bonferroni_correction(pval, alpha=alpha)
+threshold_bonferroni = stats.t.ppf(1.0 - alpha / n_tests, n_samples - 1)
+
+reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='indep')
+threshold_fdr = np.min(np.abs(T)[reject_fdr])
+
+###############################################################################
+# Plot
+times = 1e3 * epochs.times
+
+import pylab as pl
+pl.close('all')
+pl.plot(times, T, 'k', label='T-stat')
+xmin, xmax = pl.xlim()
+pl.hlines(threshold_uncorrected, xmin, xmax, linestyle='--', colors='k',
+ label='p=0.05 (uncorrected)', linewidth=2)
+pl.hlines(threshold_bonferroni, xmin, xmax, linestyle='--', colors='r',
+ label='p=0.05 (Bonferroni)', linewidth=2)
+pl.hlines(threshold_fdr, xmin, xmax, linestyle='--', colors='b',
+ label='p=0.05 (FDR)', linewidth=2)
+pl.legend()
+pl.xlabel("Time (ms)")
+pl.ylabel("T-stat")
+pl.show()
diff --git a/mne/stats/__init__.py b/mne/stats/__init__.py
index 68488ea..4211030 100644
--- a/mne/stats/__init__.py
+++ b/mne/stats/__init__.py
@@ -1,3 +1,4 @@
from .permutations import permutation_t_test
from .cluster_level import permutation_cluster_test, \
permutation_cluster_1samp_test
+from .multi_comp import fdr_correction, bonferroni_correction
diff --git a/mne/stats/multi_comp.py b/mne/stats/multi_comp.py
new file mode 100644
index 0000000..3e508c6
--- /dev/null
+++ b/mne/stats/multi_comp.py
@@ -0,0 +1,98 @@
+# Authors: Josef Pktd and example from H Raja and rewrite from Vincent Davis
+# Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# Code borrowed from statsmodels
+#
+# License: BSD (3-clause)
+
+import numpy as np
+
+
+def _ecdf(x):
+ '''no frills empirical cdf used in fdrcorrection
+ '''
+ nobs = len(x)
+ return np.arange(1, nobs + 1) / float(nobs)
+
+
+def fdr_correction(pvals, alpha=0.05, method='indep'):
+ """P-value correction with False Discovery Rate (FDR)
+
+ Correction for multiple comparison using FDR.
+
+ This covers Benjamini/Hochberg for independent or positively correlated and
+ Benjamini/Yekutieli for general or negatively correlated tests.
+
+ Parameters
+ ----------
+ pvals : array_like
+ set of p-values of the individual tests.
+ alpha : float
+ error rate
+ method : {'indep', 'negcorr')
+ If 'interp' it implements Benjamini/Hochberg for independent or if
+ 'negcorr' it corresponds to Benjamini/Yekutieli.
+
+ Returns
+ -------
+ reject : array, bool
+ True if a hypothesis is rejected, False if not
+ pval_corrected : array
+ pvalues adjusted for multiple hypothesis testing to limit FDR
+
+ Notes
+ -----
+ Reference:
+ Genovese CR, Lazar NA, Nichols T.
+ Thresholding of statistical maps in functional neuroimaging using the false
+ discovery rate. Neuroimage. 2002 Apr;15(4):870-8.
+ """
+ pvals = np.asarray(pvals)
+
+ pvals_sortind = np.argsort(pvals)
+ pvals_sorted = pvals[pvals_sortind]
+ sortrevind = pvals_sortind.argsort()
+
+ if method in ['i', 'indep', 'p', 'poscorr']:
+ ecdffactor = _ecdf(pvals_sorted)
+ elif method in ['n', 'negcorr']:
+ cm = np.sum(1. / np.arange(1, len(pvals_sorted) + 1))
+ ecdffactor = _ecdf(pvals_sorted) / cm
+ else:
+ raise ValueError("Method should be 'indep' and 'negcorr'")
+
+ reject = pvals_sorted < ecdffactor * alpha
+ if reject.any():
+ rejectmax = max(np.nonzero(reject)[0])
+ else:
+ rejectmax = 0
+ reject[:rejectmax] = True
+
+ pvals_corrected_raw = pvals_sorted / ecdffactor
+ pvals_corrected = np.minimum.accumulate(pvals_corrected_raw[::-1])[::-1]
+ pvals_corrected[pvals_corrected > 1.0] = 1.0
+ return reject[sortrevind], pvals_corrected[sortrevind]
+
+
+def bonferroni_correction(pval, alpha=0.05):
+ """P-value correction with Bonferroni method
+
+ Parameters
+ ----------
+ pvals : array_like
+ set of p-values of the individual tests.
+ alpha : float
+ error rate
+
+ Returns
+ -------
+ reject : array, bool
+ True if a hypothesis is rejected, False if not
+ pval_corrected : array
+ pvalues adjusted for multiple hypothesis testing to limit FDR
+
+ """
+ pval = np.asarray(pval)
+ pval_corrected = pval / float(pval.size)
+ reject = pval < alpha
+ return reject, pval_corrected
diff --git a/mne/stats/tests/test_multi_comp.py b/mne/stats/tests/test_multi_comp.py
new file mode 100644
index 0000000..34eb1b5
--- /dev/null
+++ b/mne/stats/tests/test_multi_comp.py
@@ -0,0 +1,36 @@
+import numpy as np
+from numpy.testing import assert_almost_equal
+from nose.tools import assert_true
+from scipy import stats
+
+from mne.stats import fdr_correction, bonferroni_correction
+
+
+def test_multi_pval_correction():
+ """Test pval correction for multi comparison (FDR and Bonferroni)
+ """
+ rng = np.random.RandomState(0)
+ X = rng.randn(10, 10000)
+ X[:, :50] += 4.0 # 50 significant tests
+ alpha = 0.05
+
+ T, pval = stats.ttest_1samp(X, 0)
+
+ n_samples, n_tests = X.shape
+ 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)
+
+ fwer = np.mean(reject_bonferroni)
+ assert_almost_equal(fwer, alpha, 1)
+
+ reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='indep')
+ thresh_fdr = np.min(np.abs(T)[reject_fdr])
+ assert_true(0 <= (reject_fdr.sum() - 50) <= 50 * 1.05)
+ assert_true(thresh_uncorrected <= thresh_fdr <= thresh_bonferroni)
+
+ reject_fdr, pval_fdr = fdr_correction(pval, alpha=alpha, method='negcorr')
+ thresh_fdr = np.min(np.abs(T)[reject_fdr])
+ assert_true(0 <= (reject_fdr.sum() - 50) <= 50 * 1.05)
+ assert_true(thresh_uncorrected <= thresh_fdr <= thresh_bonferroni)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list