Statistical Analysis of Event-Related Potentials

Background

The most straightforward approach to assess differences between conditions is to compare the average event-related fields (ERFs) either visually or by computing a difference ERF, which is the channel-by-channel subtraction of two conditions. The difference ERF highlights where and when two conditions diverge.

For MEG magnetometers, topomaps of individual ERFs and their difference show a characteristic dipolar pattern: a neural current dipole generates a magnetic field that exits the skull on one side of the source and re-enters on the other, producing adjacent positive and negative poles. This contrasts with EEG, where the resistive skull smears the electric potential into a broader, less focal distribution.

Exercises

In the following exercises, you are going to compute ERFs for individual conditions, compute and visualize their difference, and compare them at specific channels over auditory and visual cortex. Here are the code snippets you will need:

Example: Compute and plot the ERF for all "auditory/left" epochs.

erf = epochs["auditory/left"].average()
erf.plot();

Exercise: Compute and plot the ERF for all "auditory/right" epochs.

Solution

erf = epochs["auditory/right"].average()
erf.plot();

There clearly is a difference between the responses to auditory stimuli to the left and right ears. To compare those responses we first have to compute the ERF for every condition. Run the cell below to compute the ERFs with by_event_type=True and convert the result to a dictionary. We can index the resulting dictionary of erfs like we did with the epochs.

erfs = epochs.average(by_event_type=True)
erfs = {e.comment: e for e in erfs}
erfs

{'auditory/left': <Evoked | 'auditory/left' (average, N=55), -0.1998 – 0.49949 s, baseline off, 102 ch, ~3.1 MiB>,
 'auditory/right': <Evoked | 'auditory/right' (average, N=55), -0.1998 – 0.49949 s, baseline off, 102 ch, ~3.1 MiB>,
 'visual/left': <Evoked | 'visual/left' (average, N=55), -0.1998 – 0.49949 s, baseline off, 102 ch, ~3.1 MiB>,
 'visual/right': <Evoked | 'visual/right' (average, N=55), -0.1998 – 0.49949 s, baseline off, 102 ch, ~3.1 MiB>}

Example: Use the combine_evoked() function with weights=[1, -1] to subtract the ERFs for the "auditory/left" and "auditory/right" conditions and use the .plot_joint() method to visualize the difference.

erf_diff = mne.combine_evoked( 
    [erfs["auditory/left"], erfs["auditory/right"]],
    weights=[1, -1]
    )
erf_diff.plot_joint();

Projections have already been applied. Setting proj attribute to True.

Exercise: Use the combine_evoked() function with weights=[1, -1] to subtract the ERFs for the "visual/left" and "visual/right" conditions and use the .plot_joint() method to visualize the difference.

Which region shows the largest difference ERF?

Solution

erf_diff = mne.combine_evoked( 
    [erfs["visual/left"], erfs["visual/right"]],
    weights=[1, -1]
    )
erf_diff.plot_joint(times=[0.09, 0.125]);

Projections have already been applied. Setting proj attribute to True.

Example: Compare the ERFs of the "auditory/left" and "auditory/right" at the channel ["MEG 1411"]. This channel is located over the right temporal region and picks up activity generated in the right auditory cortex. Which stimulus produces a larger response at this channel?

mne.viz.plot_compare_evokeds(
    [erfs["auditory/left"], erfs["auditory/right"]],
    picks=["MEG 1411"]
    );

Exercise: Compare the ERFs of the "visual/left" and "visual/right" at the channel ["MEG 1731"]. This channel is located over the left occipital region and picks up activity generated in the left visual cortex. How does the left/right asymmetry of the visual response compare to the auditory response from the previous example?

Solution

mne.viz.plot_compare_evokeds(
    [erfs["visual/left"], erfs["visual/right"]],
    picks=["MEG 1731"]
    );

Background

The one-sample t-test tests whether the mean of a sample is significantly different from a reference value — in our case, whether the mean difference between two conditions is different from zero. The t-statistic is computed as:

where $\bar{x}$ is the sample mean, $\mu_0$ is the reference value (here 0), $s$ is the sample standard deviation, and $n$ is the number of observations (trials). The t-statistic therefore measures how many standard errors the mean is away from zero. The p-value is the probability of observing a t-statistic at least as extreme as the measured one, under the null hypothesis. By convention, effects with $p < 0.05$ are considered statistically significant.

By computing a t-test at every time point, we obtain a t-value time course that shows when the two conditions diverge. Extending this to all channels gives a t-value and p-value for every channel and time point, which can be visualized as a topographic map to show where the effect is significant.

Exercises

In the following exercises, you will extract epoch data as NumPy arrays, compute t-tests at single channels and across all channels, and visualize the results as time courses and topographic maps. Here are the code snippets you will need:

Example: Get the data for all epochs in the "auditory/left" and "auditory/right" conditions at the central channel "MEG 0611", compute their difference and plot the average difference.

Hint: diff.mean(axis=0).T averages the difference across epochs and transposes the result so it can be passed to plt.plot().

data_aud_l = epochs["auditory/left"].get_data(picks=["MEG 0611"])
data_aud_r = epochs["auditory/right"].get_data(picks=["MEG 0611"])
diff = data_aud_l - data_aud_r

plt.plot(epochs.times, diff.mean(axis=0).T)
plt.ylabel("Amplitude [T]")
plt.xlabel("Time [s]")

Exercise: Get the data for all epochs in the "auditory/left" and "auditory/right" conditions at the right temporal channel ["MEG 1411"], compute their difference and plot the average difference.

Solution

data_aud_l = epochs["auditory/left"].get_data(picks=["MEG 1411"])
data_aud_r = epochs["auditory/right"].get_data(picks=["MEG 1411"])
diff = data_aud_l - data_aud_r

plt.plot(epochs.times, diff.mean(axis=0).T)
plt.ylabel("Amplitude [T]")
plt.xlabel("Time [s]")

Example: Perform a t-test to test whether the difference is significantly different from 0 at every point in time and plot the resulting t-values across time.

How does this curve compare to the average difference between the conditions from the previous exercise?

t, p = ttest_1samp(diff, popmean=0)

plt.plot(epochs.times, t.T)
plt.ylabel("T-Value [a.u.]")
plt.xlabel("Time [s]")

Exercise: The cell below computes the difference between epochs in the "visual/left" and "visual/right" conditions. Perform a t-test to test whether the difference is significantly different from 0 and plot the t-values across time.

data_vis_l = epochs["visual/left"].get_data(picks=["MEG 1731"])
data_vis_r = epochs["visual/right"].get_data(picks=["MEG 1731"])
diff = data_vis_l - data_vis_r

Solution

t, p = ttest_1samp(diff, popmean=0)

plt.plot(epochs.times, t.T)
plt.ylabel("T-Value [a.u.]")
plt.xlabel("Time [s]")

Example: Plot the t-values and highlight times where the p-value is lower than p_thresh=0.05 in red.

Hint: (p < p_thresh).any(axis=0) checks if the p-value at ANY channel is smaller than the threshold.

p_thresh= 0.05
plt.plot(epochs.times, t.T)
plt.fill_between(epochs.times, t.min(), t.max(), where=(p < p_thresh).any(axis=0), alpha=0.2, color="r")

Exercise: Plot the t-values and highlight times where the p-value is lower than 0.001 in red.

Solution

p_thresh= 0.001
plt.plot(epochs.times, t.T)
plt.fill_between(epochs.times, t.min(), t.max(), where=(p < p_thresh).any(axis=0), alpha=0.2, color="r")

The cell below computes the difference between the "visual/left" and "visual/right" conditions for ALL channels and runs a t-test. This gives us a t-value and p-value for every channel and time point.

data_vis_l = epochs["visual/left"].get_data()
data_vis_r = epochs["visual/right"].get_data()
diff = data_vis_l - data_vis_r
t, p = ttest_1samp(diff, popmean=0)
t.shape, p.shape

((102, 106), (102, 106))

Example: Plot the topographical distribution of t-values at time=0.05 seconds and mask all channels where the p-value is smaller than 0.05.

time = 0.05
p_thresh = 0.05
idx = np.argmin(np.abs(epochs.times-time))
mne.viz.plot_topomap(t[:, idx], epochs.info, size=4, mask=p[:,idx]<p_thresh);

Exercise: Plot the topographical distribution of t-values at time=0.1 seconds and mask all channels where the p-value is smaller than 0.001.

Solution

time = 0.1
p_thresh = 0.001
idx = np.argmin(np.abs(epochs.times-time))
mne.viz.plot_topomap(t[:, idx], epochs.info, size=4, mask=p[:,idx]<p_thresh);

Exercise: How low would you have to set p_thresh to only mask the channels that cover the prominent dipole over the occipital region?

Note: This is just for demonstration purposes, p-values should always be selected independently of the data!

Solution

time = 0.1
p_thresh = 0.00001
idx = np.argmin(np.abs(epochs.times-time))
mne.viz.plot_topomap(t[:, idx], epochs.info, size=4, mask=p[:,idx]<p_thresh);

Background

When performing a statistical test at every time point and channel, the probability of obtaining at least one false positive increases with the number of tests. If we conduct $m$ independent tests each at significance level $\alpha$, the probability of at least one false positive across all tests is $1 - (1 - \alpha)^m$, which approaches 1 as $m$ grows large.

The Bonferroni correction controls the family-wise error rate (FWER) — the probability of making even a single false positive across all tests — by dividing the significance threshold by the number of tests:

This guarantees that the probability of any false positive is at most $\alpha$, regardless of how many tests are performed. The correction becomes increasingly strict as $m$ grows: with 102 channels and 106 time points, the corrected threshold is $\alpha / 10812 \approx 4.6 \times 10^{-6}$, meaning only very strong effects will survive.

Exercises

In the following exercises, you will apply the Bonferroni correction to the p-values from the previous section and visualize how the corrected threshold changes with the number of tests. Here are the code snippets you will need:

The code below tests the difference between the "auditory/left" and "auditory/right" conditions at the right temporal channel "MEG 1441" and plots the t-values while highlighting all samples where the uncorrected p-value is smaller than 0.05.

data_aud_l = epochs["auditory/left"].get_data(picks=["MEG 1441"])
data_aud_r = epochs["auditory/right"].get_data(picks=["MEG 1441"])
diff = data_aud_l - data_aud_r
t, p = ttest_1samp(diff, popmean=0)

p_thresh = 0.05
plt.plot(epochs.times, t.T)
plt.fill_between(epochs.times, t.min(), t.max(), where=(p < p_thresh).any(axis=0), alpha=0.2, color="r")

Exercise: Apply the bonferroni_correction to the p-values then plot the t-values and highlight all samples where the Bonferroni-corrected p-value is smaller than 0.05.

Solution

_, p_corr = mne.stats.bonferroni_correction(p)
plt.plot(epochs.times, t.T)
plt.fill_between(epochs.times, t.min(), t.max(), where=(p_corr < p_thresh).any(axis=0), alpha=0.2, color="r")

Exercise: The code below tests the difference between "auditory/left" and "auditory/right" for all channels. Apply the bonferroni_correction and plot the topographic map of t-values at time=0.08 while masking all channels where the Bonferroni-corrected p-value is lower than 0.05.

data_aud_l = epochs["auditory/left"].get_data()
data_aud_r = epochs["auditory/right"].get_data()
diff = data_aud_l - data_aud_r
t, p = ttest_1samp(diff, popmean=0)

Solution

_, p_corr = mne.stats.bonferroni_correction(p)
time = 0.08
p_thresh = 0.05
idx = np.argmin(np.abs(epochs.times-time))
mne.viz.plot_topomap(t[:, idx], epochs.info, size=4, mask=p_corr[:,idx]<p_thresh);

Exercise: The code below plots the distribution of p-values and marks both the uncorrected and Bonferroni-corrected thresholds. Increase n_channels and n_timepoints to reflect more realistic numbers of tests and observe how the corrected threshold changes and how many results survive.

Note: The hemispheric asymmetry for unilateral auditory stimuli is a very strong effect and may survive even very strict corrections.

n_channels = 1
n_timepoints = 50
p_thresh = 0.05
p_thresh_adjusted = p_thresh / (n_channels * n_timepoints)

bins = np.logspace(np.log10(p.min()), 0, 50)
plt.hist(p.flatten(), bins=bins)
plt.axvline(p_thresh, color="c", label="Original Threshold")
plt.axvline(p_thresh_adjusted, color="r", label="Bonferroni-Corrected")
plt.xscale("log")
plt.xlabel("p-value")
plt.ylabel("Count")
plt.yscale("log")
plt.title(f"Number of Significant Results: {(p.flatten()<p_thresh_adjusted).sum()}")
plt.legend()

Background

The permutation cluster test is a non-parametric approach that exploits the spatial and temporal structure of the data. Rather than testing each time point and channel independently, it groups adjacent significant time points and channels into clusters and assesses their significance as a whole. By assessing the cluster as a whole, the permutation cluster test avoids the multiple comparisons problem entirely (see Maris & Oostenveld, 2007). The procedure has three steps:

1. Form clusters: compute t-statistics at each time point and channel and group adjacent elements that exceed a threshold into clusters. The sum of the t-values within a cluster is the cluster statistic.
2. Estimate the null distribution: randomly shuffle the trials between conditions and find the largest cluster in the shuffled data. Repeat this across many permutations to get the distribution of maximum cluster statistics under the null hypothesis.
3. Compute cluster p-values: the p-value is the proportion of permutations in which the maximum null cluster statistic was at least as large as the observed cluster statistic. A small p-value means that a cluster of the given size is unlikely under randomly shuffled data.

Once we have established the significance of a cluster, we can visualize its spatial and temporal extent to link it to the underlying neural phenomena.

Exercises

In the following exercises, you will run the spatio-temporal permutation cluster test to compare two conditions across all channels and time points, and visualize the significant clusters. Here are the code snippets you will need:

To run the cluster test we have to prepare a list X that contains the epochs for the conditions we want to compare. Because the spatio_temporal_cluster_test() function expects the data to be shaped (n_epochs, n_times, n_channels) we have to .transpose() the array.

data_aud_l = epochs["auditory/left"].get_data()
data_aud_r = epochs["auditory/right"].get_data()
X = [data_aud_l.transpose(0, 2, 1), data_aud_r.transpose(0, 2, 1)]

The cluster test needs to know which channels count as adjacent. We can use the find_ch_adjacency() function to get this information and plot the adjacency matrix.

adjacency, ch_names = mne.channels.find_ch_adjacency(epochs.info, ch_type="mag")
mne.viz.plot_ch_adjacency(epochs.info, adjacency, ch_names);

Reading adjacency matrix for neuromag306mag.

Finally, we need to set a t-value threshold: times and channels where the t-value exceeds the threshold will be grouped into clusters. Here, we choose the value from the t-distribution t_dist that corresponds to a p-value of p=0.01.

p_thresh = 0.01
n = data_aud_l.shape[0] + data_aud_r.shape[0]
t_thresh = t_dist.ppf(1 - p_thresh / 2, df=n - 2)

Example: Run the spatio_temporal_cluster_test() on the list of epochs X with 100 permutations using the ttest_ind_no_p() function as stat_fun. Then print the smallest p-value.

t, clusters, p, _= mne.stats.spatio_temporal_cluster_test(
    X,
    n_permutations=100,
    threshold=t_thresh,
    adjacency=adjacency,
    stat_fun=mne.stats.ttest_ind_no_p,
    seed=42
)
p.min()

stat_fun(H1): min=-7.8487696973852 max=5.745549017832592
Running initial clustering …
Found 37 clusters

np.float64(0.01)

Exercise: Run the spatio_temporal_cluster_test() again but increase the number of permutations to 1000 and print the smallest p-value. Why does increasing the number of permutations change the p-value?

t, clusters, p, _= mne.stats.spatio_temporal_cluster_test(
    X,
    n_permutations=1000,
    threshold=t_thresh,
    adjacency=adjacency,
    stat_fun=mne.stats.ttest_ind_no_p,
    seed=42
)
p.min()

stat_fun(H1): min=-7.8487696973852 max=5.745549017832592
Running initial clustering …
Found 37 clusters

np.float64(0.007)

Example: Get the index of the cluster with the smallest p-value and visualize the cluster using the utils.plot_cluster_topo() function. The topographical plot shows the t-value of each channel, averaged across the cluster duration.

idx = np.argsort(p)[0]
utils.plot_cluster_topo(clusters[idx], t, p[idx], epochs)

Exercise: Get the index of the cluster with the second smallest p-value and visualize the cluster using the utils.plot_cluster_topo() function.

Solution

idx = np.argsort(p)[1]
utils.plot_cluster_topo(clusters[idx], t, p[idx], epochs)

Example: Get the index of the cluster with the smallest p-value, compute the difference ERF between the "auditory/left" and "auditory/right" conditions and visualize it using utils.plot_cluster_channels() which marks the temporal extent of the cluster.

idx = np.argsort(p)[0]
diff = mne.combine_evoked([erfs["auditory/left"], erfs["auditory/right"]], weights=(-1, 1))
utils.plot_cluster_channels(clusters[idx], p[idx], diff)

Exercise: Get the index of the cluster with the second smallest p-value and visualize the difference ERF using utils.plot_cluster_channels().

Solution

idx = np.argsort(p)[1]
utils.plot_cluster_channels(clusters[idx], p[idx], diff)

The t-value threshold determines the spatial and temporal sensitivity of the cluster test. With a lower threshold, the sensitivity is higher but clusters are more spread out. Let’s increase the threshold to the t-value that corresponds to p_thresh=0.0001.

Note: The cluster test establishes only the significance of the cluster as a whole, not the precise latency or location of the effect, see Sassenhagen & Draschkow, 2019.

p_thresh = 0.0001
n = data_aud_l.shape[0] + data_aud_r.shape[0]
t_thresh = t_dist.ppf(1 - p_thresh / 2, df=n - 2)
t_thresh

np.float64(4.040874280449325)

Exercise: Run the spatio_temporal_cluster_test again with 1000 permutations using the new t_thresh.

t, clusters, p, _= mne.stats.spatio_temporal_cluster_test(
    X,
    n_permutations=1000,
    threshold=t_thresh,
    adjacency=adjacency,
    stat_fun=mne.stats.ttest_ind_no_p,
    seed=42
)

stat_fun(H1): min=-7.8487696973852 max=5.745549017832592
Running initial clustering …
Found 8 clusters

Exercise: Use utils.plot_cluster_topo and utils.plot_cluster_channels again to visualize the spatial and temporal extent of the cluster with the second lowest p-value. How does it compare to the previous cluster?

Solution

idx = np.argsort(p)[1]
utils.plot_cluster_topo(clusters[idx], t, p[idx], epochs)

idx = np.argsort(p)[1]
utils.plot_cluster_channels(clusters[idx], p[idx], diff)

Bonus: Put everything together: run the spatio_temporal_cluster_test to compare the responses to the "visual/left" and "visual/right" conditions and plot the results.

Solution

data_vis_l = epochs["visual/left"].get_data()
data_vis_r = epochs["visual/right"].get_data()
X = [data_vis_l.transpose(0, 2, 1), data_vis_r.transpose(0, 2, 1)]

t, clusters, p, _= mne.stats.spatio_temporal_cluster_test(
    X,
    n_permutations=100,
    tail=1,
    threshold=t_thresh,
    adjacency=adjacency,
    stat_fun=mne.stats.ttest_ind_no_p,
    seed=42
)

stat_fun(H1): min=-19.19529989359747 max=19.313787724127828
Running initial clustering …
Found 8 clusters

idx = np.argsort(p)[0]
utils.plot_cluster_topo(clusters[idx], t, p[idx], epochs)

idx = np.argsort(p)[0]
diff = mne.combine_evoked([erfs["visual/left"], erfs["visual/right"]], weights=[-1, 1])
utils.plot_cluster_channels(clusters[idx], p[idx], diff)