Add trial-masked detrending

[romquentin]

Function that create a weight matrix (n_channels * n_times) with masked periods (value of zero) in order to mask the trials of interest during detrending. Related to #35

[romquentin]

Is this ok with you @nbara ?

[romquentin]

And now @nbara ? Not sure to get what is the issue.

[codecov]

Codecov Report

Merging #41 (f6ae62f) into master (852391b) will increase coverage by 0.02%.
The diff coverage is 85.71%.

@@            Coverage Diff             @@
##           master      #41      +/-   ##
==========================================
+ Coverage   78.20%   78.22%   +0.02%     
==========================================
  Files          20       20              
  Lines        2202     2209       +7     
==========================================
+ Hits         1722     1728       +6     
- Misses        480      481       +1     

Impacted Files	Coverage Δ
meegkit/detrend.py	90.14% <85.71%> (-0.23%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 852391b...f6ae62f. Read the comment docs.

[nbara]

hi @romquentin yes that's good for me but I'd like the new function to be used at least once in a unit test (tests/test_detrend.py) before we merge. That way we can make sure that it won't break silently in the future

[romquentin]

Hi @nbara, I made a test on tests/test_detrend.py. It does not crash but the detrending does not seem to work well. It seems to have the same behavior than when you don't use weights few lines before:

# Without weights – detrending fails
y, _, _ = detrend(x, 3, None, show=show)

Do you have an idea why?

Thanks,

Romain

[nbara]

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.

[nbara]

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

[nbara]

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 :

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

[nbara]

@romquentin are you still here ?

[romquentin]

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?

[nbara]

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

[romquentin]

Sounds good ! Thank you