This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit 33aa24c201044947f332dd50340b11f7ef508643 Author: Alexandre Gramfort <[email protected]> Date: Mon Mar 14 11:24:12 2011 -0400 ENH : adding 1 sample cluster level stat --- ...y => plot_cluster_1samp_test_time_frequency.py} | 54 ++++--- .../stats/plot_cluster_stats_time_frequency.py | 31 ++-- mne/stats/__init__.py | 2 +- mne/stats/cluster_level.py | 168 ++++++++++++++++----- mne/stats/tests/test_cluster_level.py | 11 +- mne/time_frequency/tfr.py | 2 +- 6 files changed, 187 insertions(+), 81 deletions(-) diff --git a/examples/stats/plot_cluster_stats_time_frequency.py b/examples/stats/plot_cluster_1samp_test_time_frequency.py similarity index 64% copy from examples/stats/plot_cluster_stats_time_frequency.py copy to examples/stats/plot_cluster_1samp_test_time_frequency.py index 04f497b..711c187 100644 --- a/examples/stats/plot_cluster_stats_time_frequency.py +++ b/examples/stats/plot_cluster_1samp_test_time_frequency.py @@ -1,13 +1,19 @@ """ -====================================================== -Non-parametric cluster statistic on single trial power -====================================================== +=============================================================== +Non-parametric 1 sample cluster statistic on single trial power +=============================================================== This script shows how to estimate significant clusters in time-frequency power estimates. It uses a non-parametric statistical procedure based on permutations and cluster level statistics. +The procedure consists in : +- extracting epochs +- compute single trial power estimates +- baseline line correct the power estimates (power ratios) +- compute stats to see if ratio deviates from 1. + """ # Authors: Alexandre Gramfort <[email protected]> # @@ -20,7 +26,7 @@ import numpy as np import mne from mne import fiff from mne.time_frequency import single_trial_power -from mne.stats import permutation_cluster_test +from mne.stats import permutation_cluster_t_test from mne.datasets import sample ############################################################################### @@ -48,38 +54,30 @@ ch_name = raw.info['ch_names'][picks[0]] # Load condition 1 event_id = 1 -epochs_condition_1 = mne.Epochs(raw, events, event_id, - tmin, tmax, picks=picks, baseline=(None, 0)) -data_condition_1 = epochs_condition_1.get_data() # as 3D matrix -data_condition_1 *= 1e13 # change unit to fT / cm - -# Load condition 1 -event_id = 2 -epochs_condition_2 = mne.Epochs(raw, events, event_id, +epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0)) -data_condition_2 = epochs_condition_2.get_data() # as 3D matrix -data_condition_2 *= 1e13 # change unit to fT / cm +data = epochs.get_data() # as 3D matrix +data *= 1e13 # change unit to fT / cm # Time vector -times = 1e3 * epochs_condition_1.times # change unit to ms +times = 1e3 * epochs.times # change unit to ms frequencies = np.arange(7, 30, 3) # define frequencies of interest Fs = raw.info['sfreq'] # sampling in Hz -epochs_power_1 = single_trial_power(data_condition_1, Fs=Fs, +epochs_power = single_trial_power(data, Fs=Fs, frequencies=frequencies, - n_cycles=2, use_fft=False) + n_cycles=3, use_fft=False) -epochs_power_2 = single_trial_power(data_condition_2, Fs=Fs, - frequencies=frequencies, - n_cycles=2, use_fft=False) - -epochs_power_1 = epochs_power_1[:,0,:,:] # only 1 channel to get a 3D matrix -epochs_power_2 = epochs_power_2[:,0,:,:] # only 1 channel to get a 3D matrix +epochs_power = epochs_power[:,0,:,:] # only 1 channel to get a 3D matrix +# do ratio with baseline power: +epochs_power /= np.mean(epochs_power[:,:,times < 0], axis=2)[:,:,None] +# null hypothesis is that the ratio is 1 (set it to 0 for stats below) +epochs_power -= 1.0 ############################################################################### # Compute statistic threshold = 6.0 T_obs, clusters, cluster_p_values, H0 = \ - permutation_cluster_test([epochs_power_1, epochs_power_2], + permutation_cluster_t_test(epochs_power, n_permutations=100, threshold=threshold, tail=0) ############################################################################### @@ -88,9 +86,9 @@ import pylab as pl pl.clf() pl.subplots_adjust(0.12, 0.08, 0.96, 0.94, 0.2, 0.43) pl.subplot(2, 1, 1) -evoked_contrast = np.mean(data_condition_1, 0) - np.mean(data_condition_2, 0) -pl.plot(times, evoked_contrast.T) -pl.title('Contrast of evoked response (%s)' % ch_name) +evoked_data = np.mean(data, 0) +pl.plot(times, evoked_data.T) +pl.title('Evoked response (%s)' % ch_name) pl.xlabel('time (ms)') pl.ylabel('Magnetic Field (fT/cm)') pl.xlim(times[0], times[-1]) @@ -108,7 +106,7 @@ pl.imshow(T_obs, cmap=pl.cm.gray, extent=[times[0], times[-1], frequencies[0], frequencies[-1]], aspect='auto', origin='lower') pl.imshow(T_obs_plot, cmap=pl.cm.jet, extent=[times[0], times[-1], - frequencies[0], frequencies[-1]], + frequencies[0], frequencies[-1]], aspect='auto', origin='lower') pl.xlabel('time (ms)') diff --git a/examples/stats/plot_cluster_stats_time_frequency.py b/examples/stats/plot_cluster_stats_time_frequency.py index 04f497b..095bef4 100644 --- a/examples/stats/plot_cluster_stats_time_frequency.py +++ b/examples/stats/plot_cluster_stats_time_frequency.py @@ -1,13 +1,20 @@ """ -====================================================== -Non-parametric cluster statistic on single trial power -====================================================== +========================================================================= +Non-parametric between conditions cluster statistic on single trial power +========================================================================= -This script shows how to estimate significant clusters -in time-frequency power estimates. It uses a non-parametric +This script shows how to compare clusters in time-frequency +power estimates between conditions. It uses a non-parametric statistical procedure based on permutations and cluster level statistics. +The procedure consists in : +- extracting epochs for 2 conditions +- compute single trial power estimates +- baseline line correct the power estimates (power ratios) +- compute stats to see if the power estimates are significantly different + between conditions. + """ # Authors: Alexandre Gramfort <[email protected]> # @@ -53,7 +60,7 @@ epochs_condition_1 = mne.Epochs(raw, events, event_id, data_condition_1 = epochs_condition_1.get_data() # as 3D matrix data_condition_1 *= 1e13 # change unit to fT / cm -# Load condition 1 +# Load condition 2 event_id = 2 epochs_condition_2 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0)) @@ -65,17 +72,23 @@ times = 1e3 * epochs_condition_1.times # change unit to ms frequencies = np.arange(7, 30, 3) # define frequencies of interest Fs = raw.info['sfreq'] # sampling in Hz +n_cycles = 1.5 epochs_power_1 = single_trial_power(data_condition_1, Fs=Fs, frequencies=frequencies, - n_cycles=2, use_fft=False) + n_cycles=n_cycles, use_fft=False) epochs_power_2 = single_trial_power(data_condition_2, Fs=Fs, frequencies=frequencies, - n_cycles=2, use_fft=False) + n_cycles=n_cycles, use_fft=False) epochs_power_1 = epochs_power_1[:,0,:,:] # only 1 channel to get a 3D matrix epochs_power_2 = epochs_power_2[:,0,:,:] # only 1 channel to get a 3D matrix - ############################################################################### + +# do ratio with baseline power: +epochs_power_1 /= np.mean(epochs_power_1[:,:,times < 0], axis=2)[:,:,None] +epochs_power_2 /= np.mean(epochs_power_2[:,:,times < 0], axis=2)[:,:,None] + +############################################################################### # Compute statistic threshold = 6.0 T_obs, clusters, cluster_p_values, H0 = \ diff --git a/mne/stats/__init__.py b/mne/stats/__init__.py index 2d810de..d61f12e 100644 --- a/mne/stats/__init__.py +++ b/mne/stats/__init__.py @@ -1,2 +1,2 @@ from .permutations import permutation_t_test -from .cluster_level import permutation_cluster_test +from .cluster_level import permutation_cluster_test, permutation_cluster_t_test diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index f1a5bbb..3a8c19e 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -8,7 +8,7 @@ import numpy as np from scipy import ndimage -from scipy.stats import percentileofscore +from scipy.stats import percentileofscore, ttest_1samp from .parametric import f_oneway @@ -119,7 +119,7 @@ def permutation_cluster_test(X, stat_fun=f_oneway, threshold=1.67, slices = [slice(splits_idx[k], splits_idx[k+1]) for k in range(len(X))] - # Step 2: If we have some clusters, repeat process on permutated data + # Step 2: If we have some clusters, repeat process on permuted data # ------------------------------------------------------------------- if len(clusters) > 0: cluster_pv = np.ones(len(clusters), dtype=np.float) @@ -148,6 +148,93 @@ def permutation_cluster_test(X, stat_fun=f_oneway, threshold=1.67, permutation_cluster_test.__test__ = False + +def permutation_cluster_t_test(X, threshold=1.67, n_permutations=1000, tail=0): + """Non-parametric cluster-level 1 sample T-test + + From a array of observations, e.g. signal amplitudes or power spectrum + estimates etc., calculate if the observed mean significantly deviates + from 0. The procedure uses a cluster analysis with permutation test + for calculating corrected p-values. Randomized data are generated with + random sign flips. + + Parameters + ---------- + X: array + Array where the first dimension corresponds to the + samples (observations). X[k] can be a 1D or 2D array (time series + or TF image) associated to the kth observation. + threshold: float + The threshold for the statistic. + n_permutations: int + The number of permutations to compute. + tail : -1 or 0 or 1 (default = 0) + If tail is 1, the statistic is thresholded above threshold. + If tail is -1, the statistic is thresholded below threshold. + If tail is 0, the statistic is thresholded on both sides of + the distribution. + + Returns + ------- + T_obs : array of shape [n_tests] + T-statistic observerd for all variables + clusters: list of tuples + Each tuple is a pair of indices (begin/end of cluster) + cluster_pv: array + P-value for each cluster + H0 : array of shape [n_permutations] + Max cluster level stats observed under permutation. + + Notes + ----- + Reference: + Cluster permutation algorithm as described in + Maris/Oostenveld (2007), + "Nonparametric statistical testing of EEG- and MEG-data" + Journal of Neuroscience Methods, Vol. 164, No. 1., pp. 177-190. + doi:10.1016/j.jneumeth.2007.03.024 + """ + X_copy = X.copy() + n_samples = X.shape[0] + shape_ones = tuple([1] * X[0].ndim) + # Step 1: Calculate T-stat for original data + # ------------------------------------------------------------- + T_obs, _ = ttest_1samp(X, 0) + + clusters, cluster_stats = _find_clusters(T_obs, threshold, tail) + + # Step 2: If we have some clusters, repeat process on permuted data + # ------------------------------------------------------------------- + if len(clusters) > 0: + cluster_pv = np.ones(len(clusters), dtype=np.float) + H0 = np.zeros(n_permutations) # histogram + for i_s in range(n_permutations): + # new surrogate data with random sign flip + signs = np.sign(0.5 - np.random.rand(n_samples, *shape_ones)) + X_copy *= signs + + # Recompute statistic on randomized data + T_obs_surr, _ = ttest_1samp(X_copy, 0) + _, perm_clusters_sums = _find_clusters(T_obs_surr, threshold, tail) + + if len(perm_clusters_sums) > 0: + H0[i_s] = np.max(perm_clusters_sums) + else: + H0[i_s] = 0 + + # for each cluster in original data, calculate p-value as percentile + # of its cluster statistics within all cluster statistics in surrogate + # data + cluster_pv[:] = [percentileofscore(H0, cluster_stats[i_cl]) + for i_cl in range(len(clusters))] + cluster_pv[:] = (100.0 - cluster_pv[:]) / 100.0 # from pct to fraction + return T_obs, clusters, cluster_pv, H0 + else: + return T_obs, np.array([]), np.array([]), np.array([]) + + +permutation_cluster_t_test.__test__ = False + # if __name__ == "__main__": # noiselevel = 30 # np.random.seed(0) @@ -176,41 +263,44 @@ permutation_cluster_test.__test__ = False # condition1[..., -3:] = np.random.randn(*shape1) * noiselevel # condition2[..., -3:] = np.random.randn(*shape2) * noiselevel # -# fs, clusters, cluster_p_values, histogram = permutation_cluster_test( -# [condition1, condition2], n_permutations=1000) +# fs, clusters, cluster_p_values, histogram = permutation_cluster_t_test( +# condition1, n_permutations=100) +# +# # fs, clusters, cluster_p_values, histogram = permutation_cluster_test( +# # [condition1, condition2], n_permutations=1000) +# +# # Plotting for a better understanding +# import pylab as pl +# pl.close('all') +# +# if condition1.ndim == 2: +# pl.subplot(211) +# pl.plot(condition1.mean(axis=0), label="Condition 1") +# pl.plot(condition2.mean(axis=0), label="Condition 2") +# pl.ylabel("signal [a.u.]") +# pl.subplot(212) +# for i_c, c in enumerate(clusters): +# c = c[0] +# if cluster_p_values[i_c] <= 0.05: +# h = pl.axvspan(c.start, c.stop-1, color='r', alpha=0.3) +# else: +# pl.axvspan(c.start, c.stop-1, color=(0.3, 0.3, 0.3), alpha=0.3) +# hf = pl.plot(fs, 'g') +# pl.legend((h, ), ('cluster p-value < 0.05', )) +# pl.xlabel("time (ms)") +# pl.ylabel("f-values") +# else: +# fs_plot = np.nan * np.ones_like(fs) +# for c, p_val in zip(clusters, cluster_p_values): +# if p_val <= 0.05: +# fs_plot[c] = fs[c] +# +# pl.imshow(fs.T, cmap=pl.cm.gray) +# pl.imshow(fs_plot.T, cmap=pl.cm.jet) +# # pl.imshow(fs.T, cmap=pl.cm.gray, alpha=0.6) +# # pl.imshow(fs_plot.T, cmap=pl.cm.jet, alpha=0.6) +# pl.xlabel('time') +# pl.ylabel('Freq') +# pl.colorbar() # -# # # Plotting for a better understanding -# # import pylab as pl -# # pl.close('all') -# # -# # if condition1.ndim == 2: -# # pl.subplot(211) -# # pl.plot(condition1.mean(axis=0), label="Condition 1") -# # pl.plot(condition2.mean(axis=0), label="Condition 2") -# # pl.ylabel("signal [a.u.]") -# # pl.subplot(212) -# # for i_c, c in enumerate(clusters): -# # c = c[0] -# # if cluster_p_values[i_c] <= 0.05: -# # h = pl.axvspan(c.start, c.stop-1, color='r', alpha=0.3) -# # else: -# # pl.axvspan(c.start, c.stop-1, color=(0.3, 0.3, 0.3), alpha=0.3) -# # hf = pl.plot(fs, 'g') -# # pl.legend((h, ), ('cluster p-value < 0.05', )) -# # pl.xlabel("time (ms)") -# # pl.ylabel("f-values") -# # else: -# # fs_plot = np.nan * np.ones_like(fs) -# # for c, p_val in zip(clusters, cluster_p_values): -# # if p_val <= 0.05: -# # fs_plot[c] = fs[c] -# # -# # pl.imshow(fs.T, cmap=pl.cm.gray) -# # pl.imshow(fs_plot.T, cmap=pl.cm.jet) -# # # pl.imshow(fs.T, cmap=pl.cm.gray, alpha=0.6) -# # # pl.imshow(fs_plot.T, cmap=pl.cm.jet, alpha=0.6) -# # pl.xlabel('time') -# # pl.ylabel('Freq') -# # pl.colorbar() -# # -# # pl.show() +# pl.show() diff --git a/mne/stats/tests/test_cluster_level.py b/mne/stats/tests/test_cluster_level.py index 3b9e3bf..8e2f8e2 100644 --- a/mne/stats/tests/test_cluster_level.py +++ b/mne/stats/tests/test_cluster_level.py @@ -1,11 +1,11 @@ import numpy as np from numpy.testing import assert_equal -from ..cluster_level import permutation_cluster_test +from ..cluster_level import permutation_cluster_test, permutation_cluster_t_test -def test_permutation_t_test(): - """Test T-test based on permutations +def test_cluster_permutation_test(): + """Test cluster level permutations tests. """ noiselevel = 20 @@ -37,3 +37,8 @@ def test_permutation_t_test(): [condition1, condition2], n_permutations=500, tail=-1) assert_equal(np.sum(cluster_p_values < 0.05), 0) + + condition1 = condition1[:,:,None] # to test 2D also + T_obs, clusters, cluster_p_values, hist = permutation_cluster_t_test( + condition1, n_permutations=500, tail=0) + assert_equal(np.sum(cluster_p_values < 0.05), 1) diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py index cd1aa67..5b933e3 100644 --- a/mne/time_frequency/tfr.py +++ b/mne/time_frequency/tfr.py @@ -216,7 +216,7 @@ def single_trial_power(epochs, Fs, frequencies, use_fft=True, n_cycles=7): return power -def induced_power(data, Fs, frequencies, use_fft=True, n_cycles=25, +def induced_power(data, Fs, frequencies, use_fft=True, n_cycles=7, n_jobs=1): """Compute time induced power and inter-trial phase-locking factor -- 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
