第7章: ERP解析¶

これまでのすべての目標が積み重なってきた、まさにこの章です。まず(EEGデータを使わずに)、エポックを平均することがなぜ隠れた信号を明らかにするのか、その直感を養います。その後、実際に適用します。Epochs を Evoked オブジェクト(ERP)に変え、それを可視化し、条件同士を比較します。

パート1: なぜ平均化がうまくいくのか(シミュレーション)¶

核心にある問題はこうです。単一試行の脳反応は、たいてい信号中の継続的なノイズ(他の脳活動、筋緊張、わずかな体動)よりずっと小さいのです。単一試行だけを見ても、反応がまったく目に見えないことがよくあります。

しかし、その反応はイベントに対して時間的に固定されている(毎回同じ潜時)のに対し、ノイズはランダム(イベントとは無相関)です。多数の試行を平均すると、ランダムなノイズは部分的に打ち消し合う一方、一貫した反応は生き残り、積み重なっていきます。その平均こそが ERP です。

import numpy as np
import matplotlib.pyplot as plt

# 日本語の文字がグラフ内で正しく表示されるよう、フォントを設定する
plt.rcParams["font.sans-serif"] = ["Yu Gothic", "Meiryo", "MS Gothic", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False

rng = np.random.default_rng(42)

sfreq = 250
tmin, tmax = -0.2, 0.8
times = np.arange(tmin, tmax, 1 / sfreq)

# 300msでピークを迎える架空の「真の」脳反応(P300のようなもの)。
# 実際には、これを直接目にすることは決してない -- ノイズに埋もれているからだ。
true_response = 3 * np.exp(-0.5 * ((times - 0.3) / 0.08) ** 2)

n_trials = 40
noise_amplitude = 8  # 意図的に信号よりずっと大きくしている -- 単一試行のEEGとしては現実的

single_trials = np.array(
    [true_response + rng.normal(0, noise_amplitude, size=times.shape) for _ in range(n_trials)]
)

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

axes[0].plot(times, single_trials.T, color="gray", alpha=0.3)
axes[0].plot(times, true_response, color="black", linewidth=2, label="真の(隠れた)反応")
axes[0].set_title(f"{n_trials}回の単一試行 (それぞれほぼノイズ)")
axes[0].legend()

axes[1].plot(times, single_trials.mean(axis=0), color="crimson", linewidth=2, label="試行間の平均")
axes[1].plot(times, true_response, color="black", linewidth=1, linestyle="--", label="真の(隠れた)反応")
axes[1].set_title("試行の平均 = ERP")
axes[1].legend()

for ax in axes:
    ax.set_xlabel("イベントからの経過時間 (秒)")
    ax.axvline(0, color="k", linewidth=0.5)
axes[0].set_ylabel("振幅 (a.u.)")
plt.tight_layout()
plt.show()

png

2つのパネルの間で、フィルタリングもクリーニングも一切行っていません。行った操作は平均化だけです。これをMNEのオブジェクトに対応づけると、上の single_trials は Epochs が保持しているもの(1行が1試行)、single_trials.mean(axis=0) はまさに epochs.average() が計算するもので、Evoked オブジェクトを生成します。第4章から第6章で行ったこと(フィルタリング、ICA、除外)はすべて、個々のエポックのノイズを減らし、この平均がより速く、よりきれいに収束するようにするためのものにすぎません。しかし、ERPを実際に生み出しているのは、この平均化のステップそのものです。

パート2: 実際のERP¶

第4章と第6章のクリーンな、エポック化されたデータを再構築します(このノートブック単体で動作するように、自己完結させています)。

import mne
from pathlib import Path

mne.set_log_level("WARNING")

sample_folder = Path(mne.datasets.sample.data_path()) / "MEG" / "sample"
raw = mne.io.read_raw_fif(sample_folder / "sample_audvis_raw.fif", preload=True)
events = mne.find_events(raw, stim_channel="STI 014")
event_id = {
    "auditory/left": 1,
    "auditory/right": 2,
    "visual/left": 3,
    "visual/right": 4,
}

raw.pick(["eeg", "eog"])
raw.filter(l_freq=1.0, h_freq=40.0)

epochs = mne.Epochs(
    raw,
    events,
    event_id=event_id,
    tmin=-0.2,
    tmax=0.5,
    baseline=(-0.2, 0.0),
    preload=True,
    reject=dict(eeg=150e-6, eog=250e-6),
)
epochs

	General
	MNE object type	Epochs
	Measurement date	2002-12-03 at 19:01:10 UTC
	Participant	Unknown
	Experimenter	MEG
	Acquisition
	Total number of events	269
	Events counts	auditory/left: 63 auditory/right: 69 visual/left: 72 visual/right: 65
	Time range	-0.200 – 0.499 s
	Baseline	-0.200 – 0.000 s
	Sampling frequency	600.61 Hz
	Time points	421
	Metadata	No metadata set
	Channels
	EEG	and
	EOG
	Head & sensor digitization	146 points
	Filters
	Highpass	1.00 Hz
	Lowpass	40.00 Hz

`epochs.average()` → `Evoked`¶

これは、上のシミュレーションの1行版を、架空の試行ではなく実際のエポックに対して実行したものです。

evoked_aud_left = epochs["auditory/left"].average()
print(evoked_aud_left)
evoked_aud_left.plot(picks="eeg")
plt.show()

<Evoked | 'auditory/left' (average, N=63), -0.1998 – 0.49949 s, baseline -0.2 – 0 s, 60 ch, ~3.1 MiB>

png

プロットの読み方: 何が「本物の」ERP成分なのか?¶

音が鳴ってから約100ms後に、明確な負の落ち込みが見えるはずです。これは古典的なN100聴覚反応です。plot_joint は、鍵となる潜時における頭皮上のトポグラフィ(電極間の電圧分布)を波形と並べて重ね合わせて表示します。ある変化が、ノイズではなく空間的に一貫した脳反応であることを確かめる、良い健全性チェックです。本物の成分は、頭皮上のマップで散らばってランダムに見えるのではなく、なめらかで焦点の定まった見え方をするはずです。

evoked_aud_left.plot_joint(picks="eeg")
plt.show()

png

成分の定量化: ピーク振幅と潜時¶

目視での確認だけでなく、get_peak() は指定した時間窓内での最大変位の、正確なチャンネル・時刻・振幅を見つけます。これは実際に報告したり、統計にかけたりする類の数値です。

ch_name, latency, amplitude = evoked_aud_left.get_peak(
    tmin=0.05, tmax=0.15, mode="neg", return_amplitude=True
)
print(f"50-150ms間で最も負に振れた変位: {amplitude * 1e6:.2f} uV ({ch_name}、{latency * 1000:.0f} ms)")

50-150ms間で最も負に振れた変位: -12.00 uV (EEG 014、90 ms)

条件の比較¶

ほとんどのERP研究の要点は、1つを見るだけでなく条件同士を比較することにあります。ここでは、左耳と右耳の聴覚クリック音を比較します。

evoked_aud_right = epochs["auditory/right"].average()

mne.viz.plot_compare_evokeds(
    dict(left=evoked_aud_left, right=evoked_aud_right),
    picks="eeg",
    combine="mean",
)
plt.show()

png

統計に関する補足¶

2つの条件の曲線を目で見比べるのは出発点であって、結論ではありません。実際のERP研究では、条件間の差が統計的に信頼できるかどうかを検証します。典型的にはクラスターベースの並べ替え検定(mne.stats.permutation_cluster_test および関連の関数)を使います。これは、相関し合う多数の時点/チャンネルを同時に検定しているという事実を、単一の比較としてではなく適切に扱う手法です。これはこのチュートリアルの範囲を超えていますが、上のような比較から実際に結論を導く準備ができたときのために、名前だけでも知っておく価値があります。

まとめ¶

epochs.average() は Evoked、つまりあなたのERPを生成します。最も重要な考え方は、平均化がランダムノイズを打ち消す一方で、一貫した時間固定の反応は生き残る、ということです。
可視化には evoked.plot() / evoked.plot_joint()、成分の定量化には evoked.get_peak() を使います。
条件の比較には mne.viz.plot_compare_evokeds() を使い、正式な比較にはクラスターベースの並べ替え統計が必要です。

次へ: 第8章 — 自分のデータに適用する