Add standardized measurement error (SME)

[cbrnr]

Fixes #10647. This PR adds a new function compute_sme() (in mne.stats.erp), which (as the name suggests) computes the analytic standardized mean error (SME). Analytic in this context means that the target measure is mean amplitude within a defined time window. Other target measures like peak amplitude or peak latency will require a different approach using bootstrapping, which is not implemented in this PR.

Here's how to use it:

import mne

from mne.stats.erp import compute_sme

fpath = mne.datasets.sample.data_path()
raw = (
    mne.io.read_raw_fif(fpath / "MEG" / "sample" / "sample_audvis_raw.fif")
    .crop(tmax=60)
)
events = mne.find_events(raw, stim_channel="STI 014")

epochs = mne.Epochs(raw.pick("eeg"), events, tmin=-0.3, tmax=0.7, preload=True)

# condition "1", all channels
compute_sme(epochs["1"], 0.3, 0.5)

# condition "2", single channel
compute_sme(epochs["2"].pick("EEG 031"), 0.3, 0.5)

# multiple conditions and channels
compute_sme(epochs[["1", "2", "3", "4"]].pick(["EEG 001", "EEG 031"]), -0.2, 0)

The SME is a data quality measure for ERP analysis. I think the current interface quite nicely integrates into the MNE API, but there is one (important) use case that I'd like to discuss. By default (i.e., without any arguments), ERPLAB divides the entire epoch duration into 100-millisecond segments and computes the SME within each segment (see here for how this looks like). This is a quick way to get some data quality metrics very quickly. How could we achieve this in MNE? Of course, it is possible to do it manually like this:

import numpy as np

tmin, tmax = epochs.tmin, epochs.tmax

for t1, t2 in zip(np.arange(tmin, tmax, 0.1), np.arange(tmin + 0.1, tmax + 0.1, 0.1)):
    sme = compute_sme(epochs["1"].pick("EEG 001"), t1, t2).item()
    print(f"Time window: {t1:4.1f} - {t2:4.1f}, SME: {sme * 1e6:.3f} µV")

This throws an error, because the last time window is out of bounds. Sure this can be fixed, but I can't think of a simple solution right now. Also, I'm not sure if it would be preferable to provide a built-in way to compute this. If so, what API would you suggest?

[hoechenberger]

is it MEAN or MEASUREMENT error?

[cbrnr]

Ha! Of course it is measurement! Thanks! It is referenced correctly everywhere in the code though I think.

[larsoner]

Would it make more sense to add it as part of #12434 ?

[cbrnr]

@larsoner not sure what you mean, I wouldn't want to wait until everything else in that PR is ready. If you meant moving the function into mne.stats.erp, I wouldn't really mind, although I think as a method it would be more readily available.

[mscheltienne]

Maybe a single method to compute different kind of metrics for epochs, instead of one method per metric to avoid clutering the object API? Just a thought, no strong feeling.

[larsoner]

If you meant moving the function into mne.stats.erp, I wouldn't really mind, although I think as a method it would be more readily available.

I'd prefer to move for now. If we start adding methods for every type of ERP measure to epochs and evoked then they will become too cluttered. I think this is roughly the direction we decided to go after a lot of discussions in #12434 . So it would be nice if this PR followed that design.

[cbrnr]

Done, the function is now available as mne.stats.erp.compute_sme().

Do you have any input/thoughts on the default data quality use case I mentioned?

[hoechenberger]

@cbrnr Would you consider adding a tiny example perhaps? just to get an idea of what the output would look like and how to interpret it?

[cbrnr]

Yes, sure, good idea!

[cbrnr]

We usually make functions available in the parent package, i.e. mne.stats.compute_sme() in this case. Do we even want that here?

[cbrnr]

And BTW, I'm totally confused by the doc error. I'm assuming it has to do with :func:`~mne.stats.erp.compute_sme` in the changelog, but I have no idea why.

[mscheltienne]

And BTW, I'm totally confused by the doc error. I'm assuming it has to do with :func:`~mne.stats.erp.compute_sme` in the changelog, but I have no idea why.

I haven't looked at the error, but it might be:
:func:`~mne.stats.compute_sme`
since sphinx requires always the top level import, and it looks like you import it in the parent module?

[larsoner]

Suggested change

	ERP-related statistics:
	ERP-related statistics
	----------------------
	:py:mod:`mne.stats.erp`:
	.. automodule:: mne.stats.erp
	:no-members:
	:no-inherited-members:
	.. currentmodule:: mne.stats.erp

[larsoner]

... or something like this would work. Locally you can iteratively do make html-noplot (will take ~10 min the first time then should be quick on follow-up runs). I can push a commit if you want

[larsoner]

I think this should probably live in mne/stats/erp/tests/test_<something>.py

[hoechenberger]

if you don't want to expose the "child" module, we might as well consider making it private

[larsoner]

We should expose it, that was the plan in #12434 that this PR is meant to follow

[hoechenberger]

but then why clutter the stats namespace and not leave it in stats.erp 🧐

[larsoner]

but then why clutter the stats namespace and not leave it in stats.erp 🧐

I agree it shouldn't live in mne.stats, the changes I suggested should move it to only mne.stats.erp

[hoechenberger]

ok cool. so then we agree :)

[cbrnr]

Before I implement the separate mne.stats.erp namespace (as a dedicated directory as opposed to just a module erp.py), I think it would be useful to clarify/discuss why the other stats functions do not get to live in their own namespace. What's the rationale behind it? Is it the number of (expected) functions?

[larsoner]

Looks good, marking for merge when green, thanks @cbrnr

[withmywoessner]

Sorry, I have not worked on #12434 in a while! My role at my lab changed quite a lot when we started a new study. I think I can have everything finished next week though!

[cbrnr]

I also wanted to add that we didn't resolve the question how to compute the SME for regularly spaced time segments (e.g. every 100ms). I think this would be important, and even though I can tackle this in a new PR, I'd really appreciate any input on how to best implement it.

[larsoner]

You mean like if you have 1 sec of 1000 Hz data and you want the SME for each 100 ms interval, i.e., 10 values? I think some suitable reshape and axis arguments should get you there

[cbrnr]

@larsoner see my example at the bottom of the initial comment. I think this is too onerous for users.

[larsoner]

Maybe allow tmin and tmax to be lists of time interval pairs?

[cbrnr]

Maybe allow tmin and tmax to be lists of time interval pairs?

Yes, that would work. I guess my main issue is that I don't know how to create 100ms segments reliably though. Even epochs defined as -0.2 to 0.8 seconds relative to an event tend to have time ranges that cannot be divided by an integer number of 100ms segments. I guess it is probably not worth implementing a logic to do that automatically, so I'd rather leave it to the users. But is there a simpler way to create such segments other than in my example (which also doesn't even check bounds, especially in the last segment).

[drammock]

is there a simpler way to create such segments other than in my example (which also doesn't even check bounds, especially in the last segment).

Seems like we have some code that is close enough to doing this already that it shouldn't be too tricky to refactor it to be useful here too. E.g. for maxwell_filter if you do tSSS it creates equal-length chunks and lumps any remainder into the last chunk. Also make_fixed_length_epochs.

So I think it wouldn't be too hard to add this capability... but I'm unsure what the API should be:

• new flexible public helper func get_time_windows or so, that returns sample numbers?
• new option to the SME function that calls private helper under the hood?
• new convenience function that wraps the SME function and adds the sliding window capability?

[cbrnr]

Great, make_fixed_length_epochs sounds very promising, I'll take a look!

Regarding the API, I don't know either. I kind of want to suggest that simply calling mne.stats.erp.compute_sme(epochs) without any start or stop arguments should probably return the SME in 100ms segments over the entire epochs duration. At least from a user perspective, this would be extremely easy to do, and I'd say given how important this use case is (and this is also how it works in ERPLAB), it's maybe not such a bad idea. But I'm open for other opinions of course!

[drammock]

I kind of want to suggest that simply calling mne.stats.erp.compute_sme(epochs) without any start or stop arguments should probably return the SME in 100ms segments over the entire epochs duration.

hmm... to me tmin=None, tmax=None feels like it should do the whole file as one window. Maybe we could take inspiration from epoching baseline param? There, None means no baselining, and (None, None) means whole-file baselining. Maybe here we could replace tmin, tmax with a single param (None or "auto" -> 100ms windows, tuple -> single window, array-like -> multiple windows)?