Today I'm interested in seeing how the source localization on the simulated data performs.
Below, I try to solve the inverse problems using three different methods employed by MNE,
MNE's namesake, dSPM, and sLORETA. I've plotted the results below.
It looks like the left hemisphere performs better than the right hemisphere for all methods,
whereas the right hemisphere picks up a lot of signal from the left. Also, in all cases,
the right amygdala looks a lot like the left cortex (in this data, there should be no
amygdala activity at all.
One way to 'score' each method is to perform a pearson's correlation of the source time
with the input time series. We would expect each hemisphere to be correlated with its input,
incorrelated with the opposite hemisphere's input, and the amygdala to be uncorrelated with
each. Based on these criteria, MNE performs the best out of the 3, but the right amygdala
still has a very high correlation with the right input data.
Some things to try out tomorrow: 1.) adding both the left and right amygdalas as source spaces
and seeing if that has an effect, 2.) trying beamformer methods.
%matplotlib inline
%load /home/alan/PythonEEG/gsoc-subcortical/spm_test_subcortical_on_simulated_data.py
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt
import mne
from mne.datasets import spm_face
from mne.time_frequency import morlet
mne.set_log_level(False)
# read in the data files
data_path = spm_face.data_path()
fname_string = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D'
sim = mne.read_evokeds(fname_string + '-sim-ave.fif')[0]
inv = mne.minimum_norm.read_inverse_operator(fname_string + '-inv.fif')
# set some parameters for the inverse operator
sfreq = sim.info['sfreq']
snr = 5.0
lambda2 = 1.0 / snr ** 2
# try 3 different methods
for method in ['MNE', 'dSPM', 'sLORETA']:
# apply inverse operator to simulated data
stc = mne.minimum_norm.apply_inverse(sim, inv, lambda2, method)
# plot the results
fig = plt.figure()
# right amygdala
ax1 = fig.add_subplot(221)
t = stc.times
x = stc.data[len(stc.vertno[0])+len(stc.vertno[1]):].mean(0)
ax1.plot(t, x)
ax1.set_title('%s Right Amygdala' % method)
# left cortex
ax2 = fig.add_subplot(222)
x = stc.lh_data.mean(0)
ax2.plot(t, x)
ax2.set_title('%s Left Hemisphere' % method)
# right cortex
ax3 = fig.add_subplot(223)
x = stc.rh_data.mean(0)
ax3.plot(t, x)
ax3.set_title('%s Right Hemisphere' % method)
# input data
ax4 = fig.add_subplot(224)
x = np.zeros((2, 600))
Ws = morlet(sfreq, [3, 10], n_cycles=[1, 1.5])
x[0][:len(Ws[0])] = np.real(Ws[0])
x[1][:len(Ws[1])] = np.real(Ws[1])
x[1] = np.roll(x[1], 80)
l1, = ax4.plot(t, x[0])
l2, = ax4.plot(t, x[1])
ax4.legend([l1, l2], ['Left', 'Right'])
ax4.set_title('%s Input Data' % method)
print method
r, p = stats.pearsonr(stc.lh_data.mean(0), x[0])
print "Left Hemi - Left Input: r=%.2f, p=%.2f" % (r, p)
r, p = stats.pearsonr(stc.lh_data.mean(0), x[1])
print "Left Hemi - Right Input: r=%.2f, p=%.2f" % (r, p)
r, p = stats.pearsonr(stc.rh_data.mean(0), x[0])
print "Right Hemi - Left Input: r=%.2f, p=%.2f" % (r, p)
r, p = stats.pearsonr(stc.rh_data.mean(0), x[1])
print "Right Hemi - Right Input: r=%.2f, p=%.2f" % (r, p)
r, p = stats.pearsonr(stc.data[len(stc.vertno[0])+len(stc.vertno[1]):].mean(0), x[0])
print "Right Amyg - Left Input: r=%.2f, p=%.2f" % (r, p)
r, p = stats.pearsonr(stc.data[len(stc.vertno[0])+len(stc.vertno[1]):].mean(0), x[1])
print "Right Amyg - Right Input: r=%.2f, p=%.2f" % (r, p)
print ''
plt.show()
No comments:
Post a Comment