[med-svn] [python-mne] 132/376: ENH : adding 1 sample cluster level stat
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:21 UTC 2015
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 <alexandre.gramfort at inria.fr>
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 <gramfort at nmr.mgh.harvard.edu>
#
@@ -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 <gramfort at nmr.mgh.harvard.edu>
#
@@ -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
More information about the debian-med-commit
mailing list