plot_cluster_stats_spatio_temporal_tris.py

"""
=================================================================
Permutation t-test on source data with spatio-temporal clustering
=================================================================

In this version of the example we're trying an analysis based on
within-label connectivity.

Watchout. This examples was written against the background of
MNE-Python PR #1503. The code and the functionality might change.
"""

# Authors: Denis Engemann <denis.engemann@gmail.com>
# License: BSD (3-clause)

print(__doc__)

import numpy as np
from numpy.random import randn
from scipy import stats as stats

import mne
from mne import io, spatial_tris_connectivity
from mne.epochs import equalize_epoch_counts
from mne.stats import spatio_temporal_cluster_1samp_test
from mne.minimum_norm import apply_inverse, read_inverse_operator
from mne.datasets import sample

###############################################################################
# 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'
subjects_dir = data_path + '/subjects'

tmin = -0.2
tmax = 0.3  # Use a lower tmax to reduce multiple comparisons

#   Setup for reading the raw data
raw = io.Raw(raw_fname)
events = mne.read_events(event_fname)

###############################################################################
# Read epochs for all channels, removing a bad one
raw.info['bads'] += ['MEG 2443']
picks = mne.pick_types(raw.info, meg=True, eog=True, exclude='bads')
event_id = 1  # L auditory
reject = dict(grad=1000e-13, mag=4000e-15, eog=150e-6)
epochs1 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                     baseline=(None, 0), reject=reject, preload=True)

event_id = 3  # L visual
epochs2 = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                     baseline=(None, 0), reject=reject, preload=True)

#    Equalize trial counts to eliminate bias (which would otherwise be
#    introduced by the abs() performed below)
equalize_epoch_counts([epochs1, epochs2])

###############################################################################
# Transform to source space

fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
snr = 3.0
lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)
inverse_operator = read_inverse_operator(fname_inv)
sample_vertices = [s['vertno'] for s in inverse_operator['src']]

fname_label = data_path + '/MEG/sample/labels/Aud-rh.label'

#    Prepare label
label = mne.read_label(fname_label)
label.values.fill(1.0)
label = label.morph('sample',  subject_to='fsaverage',
                    subjects_dir=subjects_dir)

#    Let's average and compute inverse to speed things up
condition1, condition2 = [(apply_inverse(e.average(), inverse_operator,
                                         lambda2, method)
                           .morph('fsaverage', subjects_dir=subjects_dir)
                           .in_label(label=label)
                           .crop(0, None))
                          for e in [epochs1, epochs2]]

###############################################################################
# Transform to common cortical space

n_vertices_sample, n_times = condition1.data.shape
n_subjects = 7
print('Simulating data for %d subjects.' % n_subjects)

#    Let's make sure our results replicate, so set the seed.
np.random.seed(0)

#    get connectivity on label
tris = mne.grade_to_tris(5)
tris_label = label.get_tris(tris)
connectivity = spatial_tris_connectivity(tris_label, remap_vertices=True)
label_vertices = label.get_vertices_used(np.arange(10242))

X = randn(len(label_vertices), n_times, n_subjects, 2) * 10
X[:, :, :, 0] += condition1.data[:, :, np.newaxis]
X[:, :, :, 1] += condition2.data[:, :, np.newaxis]
X = X.reshape(len(label_vertices), n_times, n_subjects, 2)
X = np.abs(X)  # only magnitude
X = X[:, :, :, 0] - X[:, :, :, 1]  # make paired contrast

###############################################################################
# Compute statistic

#    Note that X needs to be a multi-dimensional array of shape
#    samples (subjects) x time x space, so we permute dimensions
X = np.transpose(X, [2, 1, 0])

#    Now let's actually do the clustering. This can take a long time...
#    Here we set the threshold quite high to reduce computation.
p_threshold = 0.01
t_threshold = -stats.distributions.t.ppf(p_threshold / 2., n_subjects - 1)
print('Clustering.')
T_obs, clusters, cluster_p_values, H0 = clu = \
    spatio_temporal_cluster_1samp_test(X, connectivity=connectivity, n_jobs=2,
                                       threshold=t_threshold,
                                       n_permutations=1024)
#    Now select the clusters that are sig. at p < 0.05 (note that this value
#    is multiple-comparisons corrected).
good_cluster_inds = np.where(cluster_p_values < 0.05)[0]

###############################################################################
#    Visualize the clusters

from surfer import Brain

brain = Brain('fsaverage', 'rh', 'inflated')

colors = ['red', 'blue', 'purple']
times = condition1.times * 1e3

import matplotlib.pyplot as plt
fig = plt.figure()

for ii, (t_inds, s_inds) in enumerate(np.array(clusters)[good_cluster_inds]):
    color = colors[ii]
    s_inds = np.unique(s_inds)
    t_inds = np.unique(t_inds)
    brain.add_label(
        mne.Label(vertices=label_vertices[s_inds],
                  subject='fsaverage',
                  hemi='rh').smooth(subjects_dir=subjects_dir),
        color=color
    )

    plt.plot(times, X.mean(0)[:, s_inds].mean(-1), color=color,
             label='cluster-{}'.format(ii + 1))
    plt.fill_betweenx((-15, 30), times[t_inds][0], times[t_inds][-1],
                      color='gray', alpha=0.3)
plt.ylim([-15, 30])
plt.xlim(times[[0, -1]])
plt.legend(loc='upper right')
plt.xlabel('Time (ms)')
plt.ylabel('Source Activation (dSPM)')