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 <[email protected]> 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 <[email protected]> +# +# 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 <[email protected]> +# +# 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
