WIP: Inverse-transform ReceptiveField #4546

nbara · 2017-09-11T11:26:11Z

This is an attempt to implement inverse coefficient computation according to [1], which has been discussed in #3884 . This is mostly useful to visualise weights for backward models (aka stimulus reconstruction).

ATM the computation is done in ReceptiveField.fit(), so I had to expose the data covariance x_xt_ in TimeDelayingRidge(). I can compute the inverse coefficients directly inside TDR if you prefer.

These are specific points where your input would be appreciated:

proper unit tests. I was thinking to compute a forward model, and compare the coefficients to an inverse-transformed backward model, but do we have any idea how close these two sets of coefficient should be numerically?
API compatibility with the decoding framework in MNE (cf ENH: More features for ReceptiveField #3884 (comment)). Since the inverse coefs are computed by default here, is it just a matter of renaming inverse_coefs_ to patterns_ ?
Do I need to handle intercept explicitly before the final dot product?
Update example

Thanks in advance.

[1] | Haufe, S., Meinecke, F., Görgen, K., Dähne, S., Haynes, J.-D., Blankertz, B., & Bießmann, F. (2014). On the interpretation of weight vectors of linear models in multivariate neuroimaging. NeuroImage, 87, 96–110. doi:10.1016/j.neuroimage.2013.10.067

codecov-io · 2017-09-11T12:42:13Z

Codecov Report

Merging #4546 into master will decrease coverage by 1.19%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master    #4546     +/-   ##
=========================================
- Coverage   87.05%   85.85%   -1.2%     
=========================================
  Files         349      349             
  Lines       65780    65830     +50     
  Branches    10203    10208      +5     
=========================================
- Hits        57263    56518    -745     
- Misses       5634     6456    +822     
+ Partials     2883     2856     -27

kingjr · 2017-09-11T13:34:31Z

I would change x_xt_ -for cov_

I can compute the inverse coefficients directly inside TDR if you prefer.

+1 for computing them directly. How much longer does it take?

I was thinking to compute a forward model, and compare the coefficients to an inverse-transformed backward model, but do we have any idea how close these two sets of coefficient should be numerically?

Should be identical if no regularization is applied

API compatibility with the decoding framework in MNE (cf #3884 (comment)). Since the inverse coefs are computed by default here, is it just a matter of renaming inverse_coefs_ to patterns_ ?

yes

Do I need to handle intercept explicitly before the final dot product?

I would just center X and Y a priori, same as in LinearModel

proper unit tests. What do you think they should be?

Check X and Y data type and dims
Check error handling when not enough data (pinv may not be solvable is p >> n and no regularization is applied)

kingjr · 2017-09-11T13:31:12Z

mne/decoding/receptive_field.py

+                                            np.linalg.inv(np.dot(y.T, y))])
+        else:
+            inv_coef = np.linalg.multi_dot([self.estimator_.x_xt_, coef,
+                                            np.linalg.inv(np.dot(y.T, y))])


You need to ensure that y is centered to 0.

Prefer using np.cov when possible to simplify equation

larsoner · 2017-09-11T13:37:37Z

+1 for computing them directly. How much longer does it take?

The autocorrelation calculation can be very slow, especially for long recordings.

nbara · 2017-09-11T13:39:27Z

Thanks for the comments @kingjr.

+1 for computing them directly. How much longer does it take?

It's exactly the same computation, the only difference is that I wouldn't expose x_xt outside of TimeDelayingRidge() (but it would be computed anyway)

kingjr · 2017-09-11T13:40:13Z

@larsoner Ok. @nbara, could you list the exact use-cases in the PR description, to identify whether we should have an option to automatically compute the inverse in the object, notably to ensure CV compatibility.

agramfort · 2017-09-11T18:10:19Z

mne/decoding/receptive_field.py

+        coef = np.reshape(self.coef_, (n_feats * n_delays, n_outputs))
+        if not isinstance(self.estimator_, TimeDelayingRidge):
+            inv_coef = np.linalg.multi_dot([np.dot(X.T, X), coef,
+                                            np.linalg.inv(np.dot(y.T, y))])


use linalg.inv from scipy.linalg not numpy

nbara · 2017-09-12T13:50:28Z

OK, I rebased after #4543. I fixed your first few comments, but I need more tests still.

I'm still unclear as to what the consensus is regarding API : are you OK with the way things are now, or should patterns_ only be computed upon request?

If it's the latter, should it be an extra parameter to ReceptiveField(patterns=True) or a separate class function altogether, such as rf.inverse_transform(X, y)?

The current code doesn't add more than a couple seconds tops, but I guess this can add up when doing CVs and such.

larsoner · 2017-09-12T15:32:30Z

I'd rather compute them because I worry about the speed of recalculating. But @kingjr understands the sklearn ramifications better, will it make things worse there?

larsoner · 2017-09-12T15:32:51Z

Sorry, I mean I'd rather store the XtX and Xy than recompute them.

larsoner · 2017-09-12T16:18:12Z

Can some example be updated to show the utility of this?

nbara · 2017-09-12T16:32:59Z

Can some example be updated to show the utility of this?

Yep, sure. Can update the existing example to show Backward modelling+inverse patterns on the same data.

kingjr · 2017-09-12T17:19:56Z

I'm in favour of ReceptiveField(patterns=True)

larsoner · 2017-09-12T17:21:51Z

Default False to save time or True for convenience?

kingjr · 2017-09-12T17:22:12Z

Default False for time

larsoner · 2017-09-14T15:47:56Z

@nbara the positive-lags PR was merged so feel free to rebase

larsoner · 2017-09-14T15:48:07Z

(and update the examples/tutorials as necessary)

nbara · 2017-09-15T17:41:09Z

@larsoner I've rebased again and updated the example. Let me know what you think.

FWIW I can't reproduce FIG5 in Crosse et al. I've tried using sklearn's Ridge regression, as well as different hyperparameter values to match the paper, but it doesn't look similar. Maybe I've missed something.

agramfort · 2017-09-16T06:29:56Z

examples/decoding/plot_receptive_field.py

-Y, _ = raw[:]  # Outputs for the model
-Y = Y.T
+EEG, _ = raw[:]  # Outputs for the model
+EEG = EEG.T


please don't rename variables as PR can become controversial for little gain.

we usually use data (not EEG) and here it is written Y as it is the standard name for predicted target variable in ML

OK. The EEG is not the predicted variable in the backward case so I thought it'd make more sense.

nbara · 2017-09-18T07:41:52Z

If you have an idea of what a good test would be to check that the inverse coefficients have the right values besides checking that np.dot(coef.T, patterns) = np.eye(...) please let me know.

PS: squeezed in PR #4205 in whats_new as I'd forgotten back then. let me know if that's not kosher.

larsoner

Looks like good progress, some comments for you.

In case others want to look (@choldgraf?):

https://5738-1301584-gh.circle-artifacts.com/0/tmp/circle-artifacts.YbzgFFu/html/auto_examples/decoding/plot_receptive_field.html

larsoner · 2017-09-18T13:41:54Z

examples/decoding/plot_receptive_field.py

-linear encoding model using the continuously-varying speech envelope to
-predict activity of a 128 channel EEG system.
+can perform a similar function along with scikit-learn. We will first fit a
+linear encoding (or forward) model using the continuously-varying speech


You are going to get an angry GIF from @choldgraf, I had a similar "encoding (or forward)" in my PR and it did not go over well

Actually I am fine with it if they're both specified, as long as encoding/decoding come first ;-)

also now I feel like a curmudgeonly old academic

larsoner · 2017-09-18T13:45:12Z

examples/decoding/plot_receptive_field.py

+
+# Initialize the model. Here the features are the EEG data. We also specify
+# ``patterns=True`` to compute forward-transformed coefficients during model
+# fitting (cf. next section). We'll use a Ridge regression estimator with a


ridge regression (no caps)

larsoner · 2017-09-18T13:46:33Z

examples/decoding/plot_receptive_field.py

+# recreate `figure 5`_ from [1]_. The backward model weights reflect the
+# channels that contribute most toward reconstructing the stimulus signal, but
+# are not directly interpretable in a neurophysiological sense. Here we also
+# look at the inversed model weights, which are easier to understand [2]_.


This last sentence should probably start another section (so the plot is separate) and also give some explanation about what you mean about being easier to understand -- one or two sentences is fine

agree - we want one figure per section

larsoner · 2017-09-18T13:48:24Z

mne/decoding/receptive_field.py

+                y = y.reshape(-1, y.shape[-1], order='F')
+            else:
+                X = X - X.mean(0, keepdims=True)
+                cov_ = np.cov(X.T)


Can you change the TimeDelayingRidge class to have a patterns option? Then this can be stored instead of being recomputed.

could create a _compute_patterns(model, X, y) that would be shared between them.

@larsoner Not sure I understand. What is being recomputed here ? I'm only computing the pattern once either way...

Do you want the pattern/inverse coef to be computed inside TDR if pattern=True ? That means there's going to be some code redundancy between receptive_field.py and TDR.py, is that OK?

@larsoner Sorry if I'm being slow but can you elaborate on what it is you want to see?

Or are you happy with the current state of the code ?

It looks like you've already changed the code to be fast for TDR so I think we're good

larsoner · 2017-09-18T13:49:21Z

mne/decoding/receptive_field.py

+            del X
+
+            # Inverse output covariance
+            inv_Y = 1. / float(n_times * n_epochs - 1)


this should go in an else clause because it's not used in the y.ndim == 2 and y.shape[1] != 1 case

larsoner · 2017-09-18T13:50:17Z

mne/decoding/tests/test_receptive_field.py

+        # assert_allclose(inv_rf.coef_, np.transpose(rf.patterns_, (1, 0, 2)),
+        #                 rtol=0.1)
+        # assert_allclose(rf.coef_, np.transpose(inv_rf.patterns_, (1, 0, 2)),
+        #                 atol=0.1)


This might only be true if regularization is 0., can you check?

yeah my sense is that the regularization will do some strange things. Plus we're including extra information in there when we do the inverse transform. If we just do encoding model ridge, we're getting weights by normalizing w/ the stimulus covariance matrix, if we do an inverted decoding model, we're getting weights normalized w/ the channels covariance matrix, then normalized again by the stimulus covariance matrix, no?

@larsoner Yeah I checked and this also fails without regularization. Like @choldgraf I tend to think that they're something going on with the regularisation, although I can't put my finger on what exactly.

larsoner · 2017-09-18T13:52:56Z

examples/decoding/plot_receptive_field.py

+mne.viz.plot_topomap(np.mean(mean_patterns[:, ix_plot], axis=1),
+                     pos=info, axes=ax[1],
+                     show=False, vmin=-max_patterns, vmax=max_patterns)
+ax[1].set(title="Forward-transformed coefficients")


Can you add the "for delay ..." so there is a better correspondence with the earlier section?

choldgraf · 2017-09-20T15:46:16Z

@kingjr left another block on this PR...I'm hesitant to merge until he unblocks it...

choldgraf · 2017-09-20T15:48:21Z

Here's a link to the rendered page so that it's a top-level link:

https://5786-1301584-gh.circle-artifacts.com/0/tmp/circle-artifacts.Z9QLSD8/html/auto_examples/decoding/plot_receptive_field.html

choldgraf · 2017-09-20T15:52:48Z

LGTM - my only question is the delays in the final plot (-.14 and -.125). I think they're negative because stimulus activity lags in front of the brain activity, which in the model will show up as positive delays (aka, from the perspective of the stimulus, brain activity in the future is predicting stimulus activity now). I think we should mention this so that it doesn't confuse people.

nbara · 2017-09-20T16:08:07Z

@choldgraf How about this as a comment just prior to specifying the lags :

We use the same lags as in [1]. Negative lags now index the relationship between the neural response and the speech envelope earlier in time, whereas positive lags would index how a unit change in the amplitude of the EEG would aﬀect later stimulus activity (obviously this should have an amplitude of zero).

kingjr · 2017-09-20T17:40:29Z

LGTM - my only question is the delays in the final plot (-.14 and -.125) ...

? IMO - should always be the past.

kingjr · 2017-09-20T17:46:15Z

@kingjr left another block on this PR...I'm hesitant to merge until he unblocks it...

I can't (because it's been rebranched?), but my comment has been addressed

He said it's okay

larsoner · 2017-09-20T17:51:08Z

@nbara I like the proposed change, let's see if @choldgraf agrees then merge

choldgraf · 2017-09-20T18:01:45Z

https://youtu.be/2fwy7v7BK2Q?t=7s

only works with estimator=Ridge() at the moment

- demean prior to computing coefficients - rename inverse_coef_ to patterns_ - scipy.linalg - some tests (incomplete)

Can't reproduce Lalor's figures.

+ squeezed in GDF support in whatsnew

- replace `backward` by `decoding` - rewrite example - fix section headers in example TODO: add patterns option to TimeDelayingRidge

nbara · 2017-09-21T07:36:40Z

@larsoner OK it's pushed, but CircleCI fails. I'm not sure why. Let me know if there's anything to do on my end.

larsoner · 2017-09-21T15:43:14Z

Okay LGTM https://5815-1301584-gh.circle-artifacts.com/0/tmp/circle-artifacts.CszqOkn/html/auto_examples/decoding/plot_receptive_field.html

larsoner · 2017-09-21T15:43:32Z

Thanks @nbara !

choldgraf · 2017-09-21T15:47:22Z

woooo way to go! nice addition!

kingjr previously requested changes Sep 11, 2017

View reviewed changes

agramfort reviewed Sep 11, 2017

View reviewed changes

nbara force-pushed the inverse-transform-rf branch from 0f5f9f5 to 72df78f Compare September 12, 2017 07:34

nbara force-pushed the inverse-transform-rf branch from 8ba21f8 to e6ced48 Compare September 12, 2017 14:24

larsoner mentioned this pull request Sep 12, 2017

ENH: More features for ReceptiveField #3884

Open

10 tasks

nbara force-pushed the inverse-transform-rf branch from a776fcd to a402d19 Compare September 14, 2017 13:41

nbara force-pushed the inverse-transform-rf branch from a402d19 to 9405aea Compare September 15, 2017 17:36

agramfort reviewed Sep 16, 2017

View reviewed changes

nbara force-pushed the inverse-transform-rf branch 4 times, most recently from 87eb7a7 to cac4223 Compare September 18, 2017 13:38

larsoner reviewed Sep 18, 2017

View reviewed changes

nbara force-pushed the inverse-transform-rf branch from 786983e to 6d1b6e4 Compare September 20, 2017 07:25

larsoner added this to the 0.15 milestone Sep 20, 2017

nbara and others added 8 commits September 21, 2017 09:03

first commit

c708452

only works with estimator=Ridge() at the moment

support TimeDelayingRidge()

58ccbf4

address comments

aef5d58

- demean prior to computing coefficients - rename inverse_coef_ to patterns_ - scipy.linalg - some tests (incomplete)

made pattern computation optional + doc

73a5844

Add example

8cf483d

Can't reproduce Lalor's figures.

rename variable in example + fix travis unicode error

33df9bd

+ squeezed in GDF support in whatsnew

address comments

0dcadce

- replace `backward` by `decoding` - rewrite example - fix section headers in example TODO: add patterns option to TimeDelayingRidge

comments

3766dbe

nbara force-pushed the inverse-transform-rf branch 2 times, most recently from 3a2310d to 9d39873 Compare September 21, 2017 07:08

doc for negative lags

4bb7955

nbara force-pushed the inverse-transform-rf branch from 9d39873 to 4bb7955 Compare September 21, 2017 07:11

larsoner merged commit ce040f0 into mne-tools:master Sep 21, 2017

nbara deleted the inverse-transform-rf branch September 22, 2017 07:13

nbara restored the inverse-transform-rf branch November 22, 2018 14:30

Uh oh!

WIP: Inverse-transform ReceptiveField #4546

WIP: Inverse-transform ReceptiveField #4546

Uh oh!

Conversation

nbara commented Sep 11, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov-io commented Sep 11, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

kingjr commented Sep 11, 2017

Uh oh!

Choose a reason for hiding this comment

Uh oh!

larsoner commented Sep 11, 2017

Uh oh!

nbara commented Sep 11, 2017

Uh oh!

kingjr commented Sep 11, 2017

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nbara commented Sep 12, 2017

Uh oh!

larsoner commented Sep 12, 2017

Uh oh!

larsoner commented Sep 12, 2017

Uh oh!

larsoner commented Sep 12, 2017

Uh oh!

nbara commented Sep 12, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kingjr commented Sep 12, 2017

Uh oh!

larsoner commented Sep 12, 2017

Uh oh!

kingjr commented Sep 12, 2017

Uh oh!

larsoner commented Sep 14, 2017

Uh oh!

larsoner commented Sep 14, 2017

Uh oh!

nbara commented Sep 15, 2017

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nbara commented Sep 18, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

larsoner left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

choldgraf Sep 18, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nbara Sep 19, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nbara Sep 20, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

nbara commented Sep 11, 2017 •

edited

Loading

codecov-io commented Sep 11, 2017 •

edited

Loading

nbara commented Sep 12, 2017 •

edited

Loading

nbara commented Sep 18, 2017 •

edited

Loading

choldgraf Sep 18, 2017 •

edited

Loading

nbara Sep 19, 2017 •

edited

Loading

nbara Sep 20, 2017 •

edited

Loading