Use SCoT to estimate effective connectivity on MNE labels

mnescot.py

"""
=============================================================
Compute PDC spectrum source space connectivity between labels
=============================================================

The connectivity is computed between 4 labels across the spectrum
between 5 and 40 Hz.
"""

# Authors: Alexandre Gramfort <gramfort@nmr.mgh.harvard.edu>
#          Martin Billinger <martin.billinger@tugraz.at>
#
# License: BSD (3-clause)

import sys
import os

try:
    import mne
except ImportError:
    sys.path.insert(1, os.path.expanduser('~/git-working/mne-python'))
    import mne

try:
    import scot
except ImportError:
    sys.path.insert(1, os.path.expanduser('~/git-working/SCoT'))
    import scot    

from mne.datasets import sample
from mne.fiff import Raw, pick_types
from mne.minimum_norm import apply_inverse_epochs, read_inverse_operator
from mne.connectivity import spectral_connectivity

from scot.backend.sklearn import VAR
from scot.connectivity import connectivity
import scot.connectivity_statistics as scs
import scot.plotting as splt
import numpy as np

data_path = sample.data_path()
subjects_dir = data_path + '/subjects'
fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
fname_raw = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_event = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'

# Load data
inverse_operator = read_inverse_operator(fname_inv)
raw = Raw(fname_raw)
events = mne.read_events(fname_event)

# Add a bad channel
raw.info['bads'] += ['MEG 2443']

# Pick MEG channels
picks = pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,
                   exclude='bads')

# Define epochs for left-auditory condition
event_id, tmin, tmax = 1, -0.2, 0.5
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,
                                                    eog=150e-6))

# Compute inverse solution and for each epoch. By using "return_generator=True"
# stcs will be a generator object instead of a list.
snr = 1.0  # use lower SNR for single epochs
lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)
stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, method,
                            pick_ori="normal", return_generator=True)

# Read some labels
names = ['Aud-lh', 'Aud-rh', 'Vis-lh', 'Vis-rh']
labels = [mne.read_label(data_path + '/MEG/sample/labels/%s.label' % name)
          for name in names]

# Average the source estimates within each label using sign-flips to reduce
# signal cancellations
src = inverse_operator['src']
label_ts = mne.extract_label_time_course(stcs, labels, src, mode='mean_flip',
                                         return_generator=False)

fmin, fmax = 5, 40.
sfreq = raw.info['sfreq']  # the sampling frequency

# rearrange data to fit scot's format
label_ts = np.asarray(label_ts).transpose(2, 1, 0)

# remove mean over epochs (evoked response) to improve stationarity
label_ts -= label_ts.mean(axis=2, keepdims=True)

# create VAR model of order 9
mvar = VAR(9)

# generate connectivity surrogates under the null-hypothesis of no connectivity
c_surrogates = scs.surrogate_connectivity('PDC', label_ts, mvar)
c0 = np.percentile(c_surrogates, 95, axis=0)
fig = splt.plot_connectivity_spectrum([c0], fs=sfreq, freq_range=[fmin, fmax], diagonal=-1)

# estimate and plot connectivity
mvar.fit(label_ts)
con = connectivity('PDC', mvar.coef, mvar.rescov)
splt.plot_connectivity_spectrum(con, fs=sfreq, freq_range=[fmin, fmax], diagonal=-1, fig=fig)

splt.show_plots()