From c52dcaef3d8113a620437eb994c399dd90f04cb8 Mon Sep 17 00:00:00 2001 From: Beliz Gokmen <48097107+BelizSertcan@users.noreply.github.com> Date: Tue, 25 Aug 2026 12:13:50 +0200 Subject: [PATCH] Add more ERP measures Co-authored-by: Jacob Woessner --- doc/api/statistics.rst | 4 + doc/changes/dev/14221.newfeature.rst | 1 + doc/changes/names.inc | 1 + mne/stats/erp.py | 379 ++++++++++++++++++++++++++- mne/stats/tests/test_erp.py | 74 +++++- mne/utils/docs.py | 13 + 6 files changed, 469 insertions(+), 3 deletions(-) create mode 100644 doc/changes/dev/14221.newfeature.rst diff --git a/doc/api/statistics.rst b/doc/api/statistics.rst index f098b7206db..c6f174152b3 100644 --- a/doc/api/statistics.rst +++ b/doc/api/statistics.rst @@ -53,6 +53,10 @@ ERP-related statistics: .. autosummary:: :toctree: ../generated/ + erp.compute_area + erp.compute_frac_area_latency + erp.compute_frac_peak_latency + erp.compute_peak erp.compute_sme Compute ``adjacency`` matrices for cluster-level statistics: diff --git a/doc/changes/dev/14221.newfeature.rst b/doc/changes/dev/14221.newfeature.rst new file mode 100644 index 00000000000..ee0ddd6f678 --- /dev/null +++ b/doc/changes/dev/14221.newfeature.rst @@ -0,0 +1 @@ +Added 4 new ERP measures: peak latency, peak amplitude, fractional peak latency, and fractional area latency. By `Jacob Woessner`_ and :newcontrib:`Beliz Sertcan Gökmen`. \ No newline at end of file diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 1fac2d421d2..414d728c9b2 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -49,6 +49,7 @@ .. _Baris Talar: https://github.com/baris-talar .. _Basile Pinsard: https://github.com/bpinsard .. _Beige Jerry Jin: https://github.com/BeiGeJin +.. _Beliz Sertcan Gökmen: https://github.com/BelizSertcan .. _Ben Beasley: https://github.com/musicinmybrain .. _Ben Tang: https://github.com/bentang18 .. _Benedikt Ehinger: https://www.benediktehinger.de diff --git a/mne/stats/erp.py b/mne/stats/erp.py index 9fc7c3bb4cd..b33f9ac81a0 100644 --- a/mne/stats/erp.py +++ b/mne/stats/erp.py @@ -5,8 +5,17 @@ # Copyright the MNE-Python contributors. import numpy as np +from scipy import integrate -from mne.utils import _validate_type +from mne._fiff.pick import _picks_to_idx +from mne.utils import ( + _check_option, + _check_pandas_installed, + _time_mask, + _validate_type, + fill_doc, + warn, +) def compute_sme(epochs, start=None, stop=None): @@ -84,3 +93,371 @@ def compute_sme(epochs, start=None, stop=None): data = epochs.get_data(tmin=start, tmax=stop) return data.mean(axis=2).std(axis=0) / np.sqrt(data.shape[0]) + + +def _compute_peak( + evoked, start=None, stop=None, picks="all", mode="abs", average=False, strict=True +): + """Locate the peak shared by compute_peak and compute_frac_peak_latency.""" + data = evoked.get_data(picks=picks) + picked_idx = _picks_to_idx(evoked.info, picks, "all", exclude=()) + ch_names = [evoked.ch_names[i] for i in picked_idx] + times = evoked.times + mask = _time_mask(times, start, stop, evoked.info["sfreq"]) + data_masked = data[:, mask] + + if average: + data = np.mean(data, axis=0, keepdims=True) + data_masked = np.mean(data_masked, axis=0, keepdims=True) + ch_names = ["Average"] + + if mode == "abs": + data_masked = np.abs(data_masked) + elif mode == "neg": + if strict and not np.any(data_masked < 0): + raise ValueError( + "No negative values encountered. Cannot operate in neg mode." + ) + data_masked = -data_masked + elif mode == "pos": + if strict and not np.any(data_masked > 0): + raise ValueError( + "No positive values encountered. Cannot operate in pos mode." + ) + + max_indices = np.argmax(data_masked, axis=1) + peak_amplitudes = data[np.arange(data.shape[0]), max_indices + np.where(mask)[0][0]] + peak_latencies = times[max_indices + np.where(mask)[0][0]] + + return peak_latencies, peak_amplitudes, data_masked, mask, times, ch_names + + +@fill_doc +def compute_peak( + evoked, + start=None, + stop=None, + picks="all", + mode="abs", + average=False, + strict=True, +): + """Compute the peak amplitude and latency of an evoked response. + + Parameters + ---------- + evoked : instance of Evoked + The evoked response object. + %(erp_evoked_start_stop)s + %(picks_all)s + mode : str + Specifies how the peak amplitude should be determined. Can be one of: + + ``'abs'`` + The peak amplitude is the maximum absolute value. + ``'neg'`` + The peak amplitude is the maximum negative value. If there are + no negative values and ``strict`` is True, a ValueError is raised. + ``'pos'`` + The peak amplitude is the maximum positive value. If there are + no positive values and ``strict`` is True, a ValueError is raised. + + Defaults to ``'abs'``. + average : bool + If True, the peak amplitude is computed by averaging the data across + channels before finding the peak. Defaults to False. + %(erp_strict)s + + Returns + ------- + peak_df : pandas.DataFrame + A DataFrame with columns 'channel', 'latency', and 'amplitude' + containing the peak amplitude and latency for each channel. If + ``average=True``, contains a single row whose 'channel' value is + ``'Average'``. + """ + pd = _check_pandas_installed(strict=True) + _check_option("mode", mode, ["abs", "neg", "pos"]) + peak_latencies, peak_amplitudes, _, _, _, channel = _compute_peak( + evoked, start, stop, picks, mode, average, strict + ) + + peak_df = pd.DataFrame( + { + "channel": channel, + "latency": peak_latencies, + "amplitude": peak_amplitudes, + } + ) + + return peak_df + + +@fill_doc +def compute_area( + evoked, + start=None, + stop=None, + picks="all", + mode="abs", + average=False, +): + """ + Compute the area under the curve of an evoked response within a given time window. + + Parameters + ---------- + evoked : instance of Evoked + The evoked response object. + %(erp_evoked_start_stop)s + %(picks_all)s + mode : str + Specifies how the area should be computed. Can be one of: + + ``'abs'`` + The absolute value of the data is used. + ``'neg'`` + Only negative values are considered. + ``'pos'`` + Only positive values are considered. + ``'intg'`` + The integral of the data is computed without rectification. + + Defaults to ``'abs'``. + average : bool + If True, the area is computed by averaging the data across channels + before integration. Defaults to False. + + Returns + ------- + area_df : pandas.DataFrame + A DataFrame with columns 'channel' and 'area' containing the area + under the curve for each channel. If ``average=True``, contains a + single row whose 'channel' value is ``'Average'``. + """ + pd = _check_pandas_installed(strict=True) + _check_option("mode", mode, ["abs", "neg", "pos", "intg"]) + data = evoked.get_data(picks=picks) + picked_idx = _picks_to_idx(evoked.info, picks, "all", exclude=()) + channel = [evoked.ch_names[i] for i in picked_idx] + times = evoked.times + mask = _time_mask(times, start, stop, evoked.info["sfreq"]) + data_masked = data[:, mask] + + if average: + data_masked = np.mean(data_masked, axis=0, keepdims=True) + channel = ["Average"] + if mode == "abs": + data_masked = np.abs(data_masked) + elif mode == "neg": + data_masked = np.clip(data_masked, None, 0) + elif mode == "pos": + data_masked = np.clip(data_masked, 0, None) + + area = integrate.trapezoid(data_masked, times[mask], axis=1) + area_df = pd.DataFrame({"channel": channel, "area": area}) + + return area_df + + +@fill_doc +def compute_frac_peak_latency( + evoked, + frac=0.5, + start=None, + stop=None, + picks="all", + mode="abs", + average=False, + strict=True, +): + """Compute the latency at which a fraction of the peak amplitude is reached. + + Parameters + ---------- + evoked : instance of Evoked + The evoked response object. + frac : float + The fraction of the peak amplitude at which to compute the latency. + Defaults to 0.5. + %(erp_evoked_start_stop)s + %(picks_all)s + mode : str + Specifies how the peak amplitude should be determined. Can be one of: + + ``'abs'`` + The peak amplitude is the maximum absolute value. + ``'neg'`` + The peak amplitude is the maximum negative value. If there are + no negative values and ``strict`` is True, a ValueError is raised. + ``'pos'`` + The peak amplitude is the maximum positive value. If there are + no positive values and ``strict`` is True, a ValueError is raised. + + Defaults to ``'abs'``. + average : bool + If True, the fractional peak latency is computed by averaging the data + across channels before finding the latency. Defaults to False. + %(erp_strict)s + + Returns + ------- + frac_peak_df : pandas.DataFrame + A DataFrame with columns 'channel', 'fractional_peak_onset', + 'fractional_peak_offset', and 'amplitude' containing the latency at + which the peak amplitude reaches the fractional threshold. If + ``average=True``, contains a single row whose 'channel' value is + ``'Average'``. + """ + pd = _check_pandas_installed(strict=True) + _check_option("mode", mode, ["abs", "neg", "pos"]) + + _, peak_amplitudes, data_masked, mask, times, channel = _compute_peak( + evoked, start, stop, picks, mode, average, strict + ) + + peak_idx = np.argmax(data_masked, axis=1) + transformed_peak = data_masked[np.arange(data_masked.shape[0]), peak_idx] + frac_amplitudes = frac * transformed_peak[:, np.newaxis] + + # Find the first time point before the peak where the signal reaches the + # fractional threshold + frac_peak_onset = np.argmax(data_masked >= frac_amplitudes, axis=1) + frac_peak_onset_latency = times[mask][frac_peak_onset] + + # Find the first time point after the peak where the signal falls back to + # the fractional threshold; NaN if it never does before the window ends + frac_peak_offset_latency = np.full(data_masked.shape[0], np.nan) + nan_channels = [] + for i in range(data_masked.shape[0]): + below_threshold = np.where(data_masked[i, peak_idx[i] :] <= frac_amplitudes[i])[ + 0 + ] + if len(below_threshold) > 0: + frac_peak_offset_latency[i] = times[mask][peak_idx[i] + below_threshold[0]] + else: + nan_channels.append(channel[i]) + if nan_channels: + warn( + f"The signal never fell back below the fractional threshold before " + f"the end of the window for {len(nan_channels)} channel(s) " + f"({', '.join(nan_channels)}); fractional_peak_offset is NaN for " + "these channels." + ) + + frac_peak_df = pd.DataFrame( + { + "channel": channel, + "fractional_peak_onset": frac_peak_onset_latency, + "fractional_peak_offset": frac_peak_offset_latency, + "amplitude": peak_amplitudes, + } + ) + + return frac_peak_df + + +@fill_doc +def compute_frac_area_latency( + evoked, + frac=0.5, + start=None, + stop=None, + picks="all", + mode="abs", + average=False, +): + """Compute the latency at which a fraction of the total area is reached. + + Parameters + ---------- + evoked : instance of Evoked + The evoked response object. + frac : float + The fraction of the area at which to compute the latency. Defaults to 0.5. + %(erp_evoked_start_stop)s + %(picks_all)s + mode : str + Specifies how the area should be computed. Can be one of: + + ``'abs'`` + The absolute value of the data is used. + ``'neg'`` + Only negative values are considered. + ``'pos'`` + Only positive values are considered. + ``'intg'`` + The integral of the data is computed without rectification. + + Defaults to ``'abs'``. + average : bool + If True, the fractional area latency is computed by averaging the data + across channels before finding the latency. Defaults to False. + + Returns + ------- + frac_area_df : pandas.DataFrame + A DataFrame with columns 'channel', 'fractional_area_latency', + and 'area' containing the latency at which the area under the curve + reaches the fractional threshold. If ``average=True``, contains a + single row whose 'channel' value is ``'Average'``. + + Notes + ----- + With ``mode='intg'`` the running signed area is not guaranteed to + increase monotonically, so for a channel whose positive and negative + portions nearly cancel, the reported latency may not correspond to any + visually meaningful point in the waveform. Only a channel whose total + area is *exactly* zero is guarded against (yielding ``NaN``); a total + area that is merely small relative to the channel's overall activity is + not. The earliest sample satisfying the fractional threshold is + returned. + """ + pd = _check_pandas_installed(strict=True) + _check_option("mode", mode, ["abs", "neg", "pos", "intg"]) + data = evoked.get_data(picks=picks) + picked_idx = _picks_to_idx(evoked.info, picks, "all", exclude=()) + channel = [evoked.ch_names[i] for i in picked_idx] + times = evoked.times + mask = _time_mask(times, start, stop, evoked.info["sfreq"]) + data_masked = data[:, mask] + times = times[mask] + if average: + data_masked = np.mean(data_masked, axis=0, keepdims=True) + channel = ["Average"] + if mode == "abs": + data_masked = np.abs(data_masked) + elif mode == "neg": + data_masked = np.clip(data_masked, None, 0) + elif mode == "pos": + data_masked = np.clip(data_masked, 0, None) + + cum_area = integrate.cumulative_trapezoid(data_masked, times, axis=1, initial=0) + area = cum_area[:, -1] + + frac_area_latency = np.full(data_masked.shape[0], np.nan) + nan_channels = [] + for ch in range(data_masked.shape[0]): + if area[ch] == 0: + # Nothing accumulated; no latency can be defined + nan_channels.append(channel[ch]) + continue + # Normalize + idx = np.where(cum_area[ch] / area[ch] >= frac)[0] + if len(idx) > 0: + frac_area_latency[ch] = times[idx[0]] + if nan_channels: + warn( + f"No area was accumulated for {len(nan_channels)} channel(s) " + f"({', '.join(nan_channels)}); fractional_area_latency is NaN " + "for these channels." + ) + + frac_area_df = pd.DataFrame( + { + "channel": channel, + "fractional_area_latency": frac_area_latency, + "area": area, + } + ) + return frac_area_df diff --git a/mne/stats/tests/test_erp.py b/mne/stats/tests/test_erp.py index d0dea27f43c..39e66adc5d8 100644 --- a/mne/stats/tests/test_erp.py +++ b/mne/stats/tests/test_erp.py @@ -4,11 +4,21 @@ from pathlib import Path +import numpy as np import pytest +from numpy.testing import assert_allclose -from mne import Epochs, read_events +from mne import Epochs, EvokedArray, create_info, read_events from mne.io import read_raw_fif -from mne.stats.erp import compute_sme +from mne.stats.erp import ( + compute_area, + compute_frac_area_latency, + compute_frac_peak_latency, + compute_peak, + compute_sme, +) + +pytest.importorskip("pandas") base_dir = Path(__file__).parents[2] / "io" / "tests" / "data" raw = read_raw_fif(base_dir / "test_raw.fif") @@ -29,3 +39,63 @@ def test_compute_sme(): compute_sme(epochs, -1.2, 0.3) with pytest.raises(ValueError, match="out of bounds"): compute_sme(epochs, -0.1, 0.8) + + +def _triangle_evoked(sfreq=1000.0): + """Return an Evoked with triangular peaks of exactly known area and latency. + + Each triangle has its apex at 0.3 s and a base spanning 0.2-0.4 s, so the + area is 0.5 * 0.2 * height, and the half-amplitude points are at 0.25 and + 0.35 s. + """ + times = np.arange(-0.1, 0.6, 1 / sfreq) + triangle = np.clip(1 - np.abs(times - 0.3) / 0.1, 0, None) + data = np.array([triangle * 1e-6, triangle * 2e-6, triangle * -1e-6]) + info = create_info(["ch0", "ch1", "ch2"], sfreq, "eeg") + return EvokedArray(data, info, tmin=times[0]) + + +def test_compute_peak(): + """Test peak amplitude and latency against a known triangular peak.""" + evoked = _triangle_evoked() + peaks = compute_peak(evoked, start=0.2, stop=0.4, mode="pos") + assert list(peaks["channel"]) == ["ch0", "ch1", "ch2"] + assert_allclose(peaks["latency"][:2], [0.3, 0.3], atol=2e-3) + assert_allclose(peaks["amplitude"][:2], [1e-6, 2e-6], rtol=1e-3) + + +def test_compute_area(): + """Test area computation against a triangle of known area.""" + evoked = _triangle_evoked() + areas = compute_area(evoked, start=0.2, stop=0.4, mode="pos") + assert_allclose(areas["area"][:2], [1e-7, 2e-7], rtol=1e-3) + + signed = compute_area(evoked, start=0.2, stop=0.4, mode="intg") + assert_allclose(signed["area"][2], -1e-7, rtol=1e-3) + + +def test_compute_frac_peak_latency(): + """Test fractional peak latency against known half-amplitude crossings.""" + evoked = _triangle_evoked() + latencies = compute_frac_peak_latency( + evoked, frac=0.5, start=0.2, stop=0.4, mode="pos" + ) + assert_allclose(latencies["fractional_peak_onset"][0], 0.25, atol=2e-3) + assert_allclose(latencies["fractional_peak_offset"][0], 0.35, atol=2e-3) + + +def test_compute_frac_area_latency(): + """Test fractional area latency on a symmetric triangle.""" + evoked = _triangle_evoked() + latencies = compute_frac_area_latency( + evoked, frac=0.5, start=0.2, stop=0.4, mode="pos", picks=["ch0"] + ) + assert_allclose(latencies["fractional_area_latency"][0], 0.3, atol=2e-3) + + # ch2 is negative-going, so mode="pos" accumulates no area for it and the + # result is NaN with a warning + with pytest.warns(RuntimeWarning, match="No area was accumulated"): + nan_latency = compute_frac_area_latency( + evoked, frac=0.5, start=0.2, stop=0.4, mode="pos", picks=["ch2"] + ) + assert np.isnan(nan_latency["fractional_area_latency"][0]) diff --git a/mne/utils/docs.py b/mne/utils/docs.py index 868ee4fe0e4..4f4aba465d7 100644 --- a/mne/utils/docs.py +++ b/mne/utils/docs.py @@ -1380,6 +1380,19 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): .. versionadded:: 1.8 """ +docdict["erp_evoked_start_stop"] = """ +start, stop : float + Start and end time of the ERP computation window in seconds. Defaults to + ``None`` and ``None``, which corresponds to the entire Evoked object. +""" + +docdict["erp_strict"] = """ +strict : bool + If True, raise an error if values are all positive when detecting + a minimum (mode='neg'), or all negative when detecting a maximum + (mode='pos'). Defaults to True. +""" + docdict["estimate_plot_psd"] = """\ estimate : str, {'power', 'amplitude'} Can be "power" for power spectral density (PSD; default), "amplitude" for