Add standardized measurement error (SME)

doc/api/statistics.rst

@@ -48,6 +48,13 @@ Non-parametric (clustering) resampling methods:
    summarize_clusters_stc
    bootstrap_confidence_interval
 
+ERP-related statistics:
+
+.. autosummary::
+   :toctree: ../generated/
+
+   erp.compute_sme
+
 Compute ``adjacency`` matrices for cluster-level statistics:
 
 .. currentmodule:: mne

doc/changes/devel/12707.newfeature.rst

@@ -0,0 +1 @@
+Add :func:`~mne.stats.erp.compute_sme` to compute the analytical standardized measurement error (SME) as a data quality measure for ERP studies, by `Clemens Brunner`_.

doc/references.bib

@@ -1158,6 +1158,17 @@ @misc{Lowry2014
   year = {2014}
 }
 
+@article{LuckEtAl2021,
+  author = {Luck, Steven J. and Stewart, Andrew X. and Simmons, Aaron M. and Rhemtulla, Mijke},
+  journal = {Psychophysiology},
+  title = {Standardized Measurement Error: A Universal Metric of Data Quality for Averaged Event-Related Potentials},
+  volume = {58},
+  number = {6},
+  pages = {e13793},
+  year = {2021},
+  doi = {10.1111/psyp.13793}
+}
+
 @article{MaessEtAl2016,
   author = {Maess, Burkhard and Schröger, Erich and Widmann, Andreas},
   doi = {10.1016/j.jneumeth.2015.12.003},

ignore_words.txt

@@ -38,3 +38,4 @@ pres
 aas
 vor
 connec
+sme

mne/stats/__init__.pyi

@@ -5,6 +5,7 @@ __all__ = [
     "bonferroni_correction",
     "bootstrap_confidence_interval",
     "combine_adjacency",
+    "erp",
     "f_mway_rm",
     "f_oneway",
     "f_threshold_mway_rm",
@@ -20,6 +21,7 @@ __all__ = [
     "ttest_1samp_no_p",
     "ttest_ind_no_p",
 ]
+from . import erp
 from ._adjacency import combine_adjacency
 from .cluster_level import (
     _st_mask_from_s_inds,

mne/stats/erp.py

@@ -0,0 +1,80 @@
+import numpy as np
+
+from mne.utils import _validate_type
+
+
+def compute_sme(epochs, start=None, stop=None):
+    """Compute standardized measurement error (SME).
+
+    The standardized measurement error :footcite:`LuckEtAl2021` can be used as a
+    universal measure of data quality in ERP studies.
+
+    Parameters
+    ----------
+    epochs : mne.Epochs
+        The epochs containing the data for which to compute the SME.
+    start : int | float | None
+        Start time (in s) of the time window used for SME computation. If ``None``, use
+        the start of the epoch.
+    stop : int | float | None
+        Stop time (in s) of the time window used for SME computation. If ``None``, use
+        the end of the epoch.
+
+    Returns
+    -------
+    sme : array, shape (n_channels,)
+        SME in given time window for each channel.
+
+    Notes
+    -----
+    Currently, only the mean value in the given time window is supported, meaning that
+    the resulting SME is only valid in studies which quantify the amplitude of an ERP
+    component as the mean within the time window (as opposed to e.g. the peak, which
+    would require bootstrapping).
+
+    References
+    ----------
+    .. footbibliography::
+
+    Examples
+    --------
+    Given an :class:`~mne.Epochs` object, the SME for the entire epoch duration can be
+    computed as follows:
+
+        >>> compute_sme(epochs)  # doctest: +SKIP
+
+    However, the SME is best used to estimate the precision of a specific ERP measure,
+    specifically the mean amplitude of an ERP component in a time window of interest.
+    For example, the SME for the mean amplitude of the P3 component in the 300-500 ms
+    time window could be computed as follows:
+
+        >>> compute_sme(epochs, start=0.3, stop=0.5)  # doctest: +SKIP
+
+    Usually, it will be more informative to compute the SME for specific conditions
+    separately. This can be done by selecting the epochs of interest as follows:
+
+        >>> compute_sme(epochs["oddball"], 0.3, 0.5)  # doctest: +SKIP
+
+    Note that the SME will be reported for each channel separately. If you are only
+    interested in a single channel (or a subset of channels), select the channels
+    before computing the SME:
+
+        >>> compute_sme(epochs.pick("Pz"), 0.3, 0.5)  # doctest: +SKIP
+
+    Selecting both conditions and channels is also possible:
+
+        >>> compute_sme(epochs["oddball"].pick("Pz"), 0.3, 0.5)  # doctest: +SKIP
+
+    In any case, the output will be a NumPy array with the SME value for each channel.
+    """
+    _validate_type(start, ("numeric", None), "start", "int or float")
+    _validate_type(stop, ("numeric", None), "stop", "int or float")
+    start = epochs.tmin if start is None else start
+    stop = epochs.tmax if stop is None else stop
+    if start < epochs.tmin:
+        raise ValueError("start is out of bounds.")
+    if stop > epochs.tmax:
+        raise ValueError("stop is out of bounds.")
+
+    data = epochs.get_data(tmin=start, tmax=stop)
+    return data.mean(axis=2).std(axis=0) / np.sqrt(data.shape[0])

mne/stats/tests/test_erp.py

@@ -0,0 +1,27 @@
+from pathlib import Path
+
+import pytest
+
+from mne import Epochs, read_events
+from mne.io import read_raw_fif
+from mne.stats.erp import compute_sme
+
+base_dir = Path(__file__).parents[2] / "io" / "tests" / "data"
+raw = read_raw_fif(base_dir / "test_raw.fif")
+events = read_events(base_dir / "test-eve.fif")
+
+
+def test_compute_sme():
+    """Test SME computation."""
+    epochs = Epochs(raw, events)
+    sme = compute_sme(epochs, start=0, stop=0.1)
+    assert sme.shape == (376,)
+
+    with pytest.raises(TypeError, match="int or float"):
+        compute_sme(epochs, "0", 0.1)
+    with pytest.raises(TypeError, match="int or float"):
+        compute_sme(epochs, 0, "0.1")
+    with pytest.raises(ValueError, match="out of bounds"):
+        compute_sme(epochs, -1.2, 0.3)
+    with pytest.raises(ValueError, match="out of bounds"):
+        compute_sme(epochs, -0.1, 0.8)