Skip to content

Add trial-masked detrending #41

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nbara merged 9 commits into from
Jun 2, 2021
Merged

Add trial-masked detrending #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
8e07355
Merge pull request #2 from nbara/master
romquentin Apr 30, 2021
6605b63
add create_maksed_weight
romquentin May 3, 2021
54c644a
Merge branch 'masked_detrend' of https://github.com/romquentin/python…
romquentin May 3, 2021
c19dfc5
Formatting docstring
romquentin May 5, 2021
ef5c5da
Merge branch 'master' into masked_detrend
nbara May 11, 2021
29fc48c
format public function
romquentin May 11, 2021
1eb44f3
Merge branch 'masked_detrend' of https://github.com/romquentin/python…
romquentin May 11, 2021
d9d371c
add test + correct dimension of weights
romquentin May 11, 2021
f6ae62f
docstring
nbara May 12, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 42 additions & 0 deletions meegkit/detrend.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -297,3 +297,45 @@ def _plot_detrend(x, y, w):
ax2.set_ylabel('ch. weights')
ax2.set_xlabel('samples')
plt.show()


def create_masked_weight(x, events, tmin, tmax, sfreq):
"""Output a weight matrix for trial-masked detrending.

Creates a (n_times, n_channels) weight matrix with masked
periods (value of zero) in order to mask the trials of interest during
detrending [1]_.

Parameters
----------
x : ndarray, shape=(n_times, n_channels)
Raw data matrix.
events : ndarray, shape=(n_events)
Time samples of the events.
tmin : float
Start time before event (in seconds).
tmax : float
End time after event (in seconds).
sfreq : float
The sampling frequency of the data.

Returns
-------
weights : ndarray, shape=(n_times, n_channels)
Weight for each channel and each time sample (zero is masked).

References
----------
.. [1] van Driel, J., Olivers, C. N., & Fahrenfort, J. J. (2021). High-pass
filtering artifacts in multivariate classification of neural time series
data. Journal of Neuroscience Methods, 352, 109080.

"""
if x.ndim != 2:
raise ValueError('The shape of x must be (n_times, n_channels)')

weights = np.ones(x.shape)
for e in events:
weights[int(e + tmin * sfreq): int(e + tmax * sfreq + 1), :] = 0

return weights
12 changes: 11 additions & 1 deletion tests/test_detrend.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""Test robust detrending."""
import numpy as np

from meegkit.detrend import regress, detrend, reduce_ringing
from meegkit.detrend import regress, detrend, reduce_ringing, create_masked_weight

from scipy.signal import butter, lfilter

Expand Down Expand Up @@ -89,6 +89,16 @@ def test_detrend(show=False):
x += 2 * np.sin(2 * np.pi * np.arange(1000) / 200)[:, None]
y, _, _ = detrend(x, 5, basis='sinusoids', show=True)

# trial-masked detrending
trend = np.linspace(0, 100, 1000)[:, None]
data = 3 * np.random.randn(*trend.shape)
data[:100, :] = 100
x = trend + data
events = np.arange(30, 970, 40)
Copy link
Owner

@nbara nbara May 12, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't looked into it in detail, but since your events are entirely equidistant from one another, the "masked weight" is not going to have much of an effect.

In other words, the trend inferred from the masked data is pretty much exactly the same. Especially since you're only trying to fit a linear trend in this example.

If we want a more striking difference between the trial mask and the vanilla implementation we would probably need a slightly more sophisticated simulated dataset.

Copy link
Owner

@nbara nbara May 12, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For instance, maybe we could add some pseudo-ERPs during the events, and see how it might mess up the detrending if we do not use a trial-mask

Copy link
Owner

@nbara nbara May 12, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider the slightly more complex example, with slightly longer 'events' to be masked out :

# trial-masked detrending
trend = np.cumsum(np.random.randn(1000, 1) + 0.1, axis=0)
data = 3 * np.random.randn(*trend.shape)
x = trend + data
events = np.arange(30, 970, 200)
tmin, tmax, sfreq = -1., 2., 20  # 3s 'events' 
for e in events:
    x[e - 10:e + 30] += 20  # add some significant deviations from baseline during events

# without trial-mask
w = np.ones(x.shape)
y1, _, _ = detrend(x, 10, w, basis='polynomials', show=show)

# with trial-mask
w = create_masked_weight(x, events, tmin, tmax, sfreq)
y2, _, _ = detrend(x, 10, w, basis='polynomials', show=show)

f, ax = plt.subplots(2, sharex=True,
                        gridspec_kw={'height_ratios': [.9, .1]})
ax[0].plot(y1, label='no trial-mask')
ax[0].plot(y2, label='trial-mask')
ax[0].legend()
ax[1].pcolormesh(w.T)
plt.show()

This produces very different detrending results :

Screenshot 2021-05-12 at 15 25 07

I'm not sure which is necessarily better, but the two lines are different, meaning the trial-mask is clearly doing something

Copy link
Owner

@nbara nbara May 28, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@romquentin are you still here ?

Copy link
Contributor Author

@romquentin romquentin May 28, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the delay. It is hard to find time working on it right now. Will try in the next coming weeks.
I will try on a real MEG dataset to see if it's worth or not. What do you think is the next step?

Copy link
Owner

@nbara nbara May 31, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To me this is fine, so I would merge. We can always revisit this issue in a future PR.

Copy link
Contributor Author

@romquentin romquentin Jun 1, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good ! Thank you

tmin, tmax, sfreq = -0.2, 0.3, 20
w = create_masked_weight(x, events, tmin, tmax, sfreq)
y, _, _ = detrend(x, 1, w, basis='polynomials', show=show)


def test_ringing():
"""Test reduce_ringing function."""
Expand Down