<jupyter_start><jupyter_text>EEG-03 Solutions<jupyter_code>import os
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')
sns.set_context('notebook', font_scale=1.25)
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Today's demonstration will introduce time-frequency analysis of EEG data using `mne-python`. We will start by loading the epoched data from last session.<jupyter_code>from mne import read_epochs

## Load data.
f = os.path.join('..','eeg-data','sub-01_task-audvis-epo.fif')
epochs = read_epochs(f, preload=True, verbose=False)<jupyter_output><empty_output><jupyter_text>Section 1: Frequency AnalysisAs a rough first-pass, we can look at the frequency composition, collapsing across time, over the scalp topography.<jupyter_code>## Define frequency bands.
bands = [(8, 12, 'Alpha'), (12, 30, 'Beta')]

## Extract epochs.
LV = epochs['LV'].copy()
RV = epochs['RV'].copy()

## Iteratively plot.
print('Left Visual Stimuli')
fig, axes = plt.subplots(1,2,figsize=(8,4))
LV.plot_psd_topomap(bands, axes=axes, verbose=False);

print('Right Visual Stimuli')
fig, axes = plt.subplots(1,2,figsize=(8,4))
RV.plot_psd_topomap(bands, axes=axes, verbose=False);<jupyter_output>Left Visual Stimuli<jupyter_text>Section 2: Time-Frequency DecompositionIn this section, we will apply spectral decomposition to the visual evoked potentials from last class. Specifically, we will use Morlet wavelet decomposition to study the time-frequency composition of the evoked potentials.First we need to define the parameters of the Morlet wavelet decomposition. This includes:- *freqs*: frequencies of interest- *n_cycles*: the number of cycles used for estimation, either globally or for each frequency<jupyter_code>## Define parameters.
freqs = np.logspace(np.log(8), np.log(30), num=10, base=np.e) # Logarithmic sampling
n_cycles = 3<jupyter_output><empty_output><jupyter_text>Next we will perform the decomposition for each condition separately.<jupyter_code>from mne.time_frequency import tfr_morlet

## Extract epochs.
LV = epochs['LV'].copy()
RV = epochs['RV'].copy()

## Perform time-frequency decomposition.
LV_power = tfr_morlet(LV, freqs, n_cycles, average=True, return_itc=False, verbose=False)
RV_power = tfr_morlet(RV, freqs, n_cycles, average=True, return_itc=False, verbose=False)<jupyter_output><empty_output><jupyter_text>It is important to perform baseline correlations for time-frequency data. There are many methods for performing baseline correction. Here we will use the log-ratio technique:$$ \text{decibel power} = 10 \times \log_{10}{\frac{\text{power}}{\text{baseline}}} $$<jupyter_code>## Define parameters.
baseline = (-0.150, -0.050)
method = 'logratio'

## Apply baseline correction.
LV_power = LV_power.apply_baseline(baseline, 'logratio')
RV_power = RV_power.apply_baseline(baseline, 'logratio')

## Convert to decibels.
LV_power_arr = 10 * LV_power.data
RV_power_arr = 10 * RV_power.data<jupyter_output>Applying baseline correction (mode: logratio)
Applying baseline correction (mode: logratio)<jupyter_text>Similar to ERP analysis, we can plot out the difference between conditions across the entire scalp.<jupyter_code>## Compute difference between conditions.
averageTFR = RV_power - LV_power

## Plot topographic difference.
averageTFR.plot_topo(baseline=None);<jupyter_output>No baseline correction applied<jupyter_text>Now let's visualize the difference between the right and left visual stimulus presentation for a left posterior channel, EEG 057.<jupyter_code>## initialize canvas.
fig, ax = plt.subplots(1,1,figsize=(16,4))

## Compute difference.
diff = RV_power_arr[epochs.ch_names.index('EEG 057')] - LV_power_arr[epochs.ch_names.index('EEG 057')]

## Plotting.
sns.heatmap(diff, center=0, vmin=-5, vmax=5, yticklabels=freqs.round(1), ax=ax)
ax.set(xticks=np.arange(epochs.times.size)[::10], xticklabels=epochs.times[::10].round(2),
       xlabel='Time (s)', ylabel='Freq (Hz)', title='RV - LV (Left Visual Cortex)')
ax.invert_yaxis()<jupyter_output><empty_output><jupyter_text>Section 3: Phase-Locked vs. Induced PowerTo compute nonphase-locked power (also sometimes called “induced” activity), we must subtract the event-related potential (ERP) from single trials.<jupyter_code>from mne.time_frequency import tfr_morlet

## Extract epochs.
LV = epochs['LV'].copy()
RV = epochs['RV'].copy()

## Subtracted evoked potential
LV = LV.subtract_evoked()
RV = RV.subtract_evoked()

## Perform time-frequency decomposition.
LV_power = tfr_morlet(LV, freqs, n_cycles, average=True, return_itc=False, verbose=False)
RV_power = tfr_morlet(RV, freqs, n_cycles, average=True, return_itc=False, verbose=False)<jupyter_output>Subtracting Evoked from Epochs
[done]
Subtracting Evoked from Epochs
[done]<jupyter_text>It is important to perform baseline correlations for time-frequency data. There are many methods for performing baseline correction. Here we will use the log-ratio technique:$$ \text{decibel power} = 10 \times \log_{10}{\frac{\text{power}}{\text{baseline}}} $$<jupyter_code>## Define parameters.
baseline = (-0.150, -0.050)
method = 'logratio'

## Apply baseline correction.
LV_power = LV_power.apply_baseline(baseline, 'logratio')
RV_power = RV_power.apply_baseline(baseline, 'logratio')

## Convert to decibels.
LV_power_arr = 10 * LV_power.data
RV_power_arr = 10 * RV_power.data<jupyter_output>Applying baseline correction (mode: logratio)
Applying baseline correction (mode: logratio)<jupyter_text>Now let's visualize the difference between the right and left visual stimulus presentation for a left posterior channel, EEG 057.<jupyter_code>## initialize canvas.
fig, ax = plt.subplots(1,1,figsize=(16,4))

## Compute difference.
diff = RV_power_arr[epochs.ch_names.index('EEG 057')] - LV_power_arr[epochs.ch_names.index('EEG 057')]

## Plotting.
sns.heatmap(diff, center=0, vmin=-5, vmax=5, yticklabels=freqs.round(1), ax=ax)
ax.set(xticks=np.arange(epochs.times.size)[::10], xticklabels=epochs.times[::10].round(2),
       xlabel='Time (s)', ylabel='Freq (Hz)', title='RV - LV (Left Visual Cortex)')
ax.invert_yaxis()<jupyter_output><empty_output><jupyter_text>Section 4: Permutation TestingFinally, we demonstrate 2d time-frequency permutation testing. Of course for permutation testing, we compute our statistics at the trial level. As such, we re-decompose our data this time at the trial level.<jupyter_code>from mne.time_frequency import tfr_morlet

## Extract epochs.
LV = epochs['LV'].copy()
RV = epochs['RV'].copy()

## Subtracted evoked potential
LV = LV.subtract_evoked()
RV = RV.subtract_evoked()

## Perform time-frequency decomposition.
LV_power = tfr_morlet(LV, freqs, n_cycles, average=False, return_itc=False, verbose=False)
RV_power = tfr_morlet(RV, freqs, n_cycles, average=False, return_itc=False, verbose=False)

## Apply baseline correction.
LV_power = LV_power.apply_baseline((-0.150, -0.050), 'logratio')
RV_power = RV_power.apply_baseline((-0.150, -0.050), 'logratio')

## Convert to decibels.
LV_power_arr = 10 * LV_power.data
RV_power_arr = 10 * RV_power.data<jupyter_output>Subtracting Evoked from Epochs
[done]
Subtracting Evoked from Epochs
[done]
Applying baseline correction (mode: logratio)
Applying baseline correction (mode: logratio)<jupyter_text>Now let's perform permutation testing. We above, we restrict our analysis to our channel of interest. We then pass both arrays, as a list, to `mne.stats.permutation_cluster_test`.<jupyter_code>from mne.stats import permutation_cluster_test

## Limit to channel of interest (EEG 057).
ix = epochs.ch_names.index('EEG 057')
LV_power_057 = LV_power_arr[:,ix]
RV_power_057 = RV_power_arr[:,ix]

## Perform permutation testing.
F_obs, clusters, cluster_pv, H0 = permutation_cluster_test([RV_power_057, LV_power_057], 
                                                           n_permutations=1024)<jupyter_output>Using a threshold of 3.908741
stat_fun(H1): min=0.000000 max=15.403558
Running initial clustering
Found 16 clusters
Permuting 1023 times...
[                                        ]   0.00%  |<jupyter_text>First let's visualize the clusters the permutation test detected.<jupyter_code>## initialize canvas.
fig, ax = plt.subplots(1,1,figsize=(16,4))

## Compute difference.
diff = diff = RV_power_057.mean(axis=0) - LV_power_057.mean(axis=0)
diff *= np.sum(clusters, axis=0)

## Plotting.
sns.heatmap(diff, center=0, vmin=-5, vmax=5, yticklabels=freqs.round(1), ax=ax)
ax.set(xticks=np.arange(epochs.times.size)[::10], xticklabels=epochs.times[::10].round(2),
       xlabel='Time (s)', ylabel='Freq (Hz)', title='RV - LV (Left Visual Cortex)')
ax.invert_yaxis()<jupyter_output><empty_output><jupyter_text>Now let's threshold to only the significant clusters.<jupyter_code>## initialize canvas.
fig, ax = plt.subplots(1,1,figsize=(16,4))

## Compute difference.
diff = diff = RV_power_057.mean(axis=0) - LV_power_057.mean(axis=0)
diff *= np.sum([cluster if p < 0.05 else np.zeros_like(cluster) 
                for p, cluster in zip(cluster_pv, clusters)], axis=0)

## Plotting.
sns.heatmap(diff, center=0, vmin=-5, vmax=5, yticklabels=freqs.round(1), ax=ax)
ax.set(xticks=np.arange(epochs.times.size)[::10], xticklabels=epochs.times[::10].round(2),
       xlabel='Time (s)', ylabel='Freq (Hz)', title='RV - LV (Left Visual Cortex)')
ax.invert_yaxis()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1259607173356106, 'fraction_numerical': 0.025085397096498718, 'mean_word_length': 6.039068369646882, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '535274', 'license': 'MIT', 'n_tokens_mistral': 3436, 'n_tokens_neox': 3128, 'n_words': 811, 'path': 'eeg-03/eeg-03-solutions.ipynb', 'repo_name': '2020-NEU502b/demos', 'size': '9368'}