mne-tools/mne-python — pull request #1407 by rgoj (2014-07-10 13:33:59+00:00)
Title: WIP/ENH: Time-frequency beamformer Wilcoxon across epochs

This is a WIP implementation of Wilcoxon tests across trials comparing each time window with a baseline period (see, e.g. [1]), which was a feature suggested by @adykstra.

I'm posting this as a WIP, because I'd like to discuss the API and the way tf_lcmv now treats baselines and noise, and also to ease contributions (I might be slow in finishing this).

`tf_lcmv` currently (before this PR) returns these quantities:
* `tf_bin_power / noise_power` or,
* `tf_bin_power / baseline_power`,
depending on whether noise_cov is calculated from empty room noise or a baseline period. (The division is not actual division of estimated source power, instead covariance is whitened and the result is equivalent to van Veen's neural activity index).

We'd like to be able to return also:
* z score from Wilcoxon test across trials (neural activity index in active time window vs. neural activity index in baseline)
* logratio of `tf_bin_power / noise_power` or logratio of `tf_bin_power / baseline_power`

We could also compute what the original 2008 Dalal 5d t-f beamforming publication uses:
* `Fdb = 10 log10 ((tf_bin_power - noise_power) / (baseline_power - noise_power)`

...and many more pseudo t, z, F options...

What keeps confusing me is how to cleanly make the user choose between these options and make clear when the baseline period is used and when only noise covariance is used. At the moment there are two new parameters for `tf_lcmv`: `baseline` and `mode`. There's a slight ambiguity, because a user could still pass baseline-calculated covariance through the `noise_cov` parameters. And if the user wants to get `tf_bin_power / baseline_power` then `noise_cov` should also be explicitly allowed to be None. Any thoughts welcome, or I can just go on and implement something completely and then we can discuss that.

Also, I had to copy the code of the wilcoxon function from scipy, because it didn't return the z scores we'd like to have. There's even an issue about it in scipy, so I'll submit a PR there, but in the meantime we'd probably have to keep the modified code locally, no? What's the best place to put it? Separate file in mne/stats or keep here in _lcmv.py where it's used?

Here's a list of TODOs before I'd see this PR as complete:
* [ ] Implement choosing baseline period.
* [ ] The wilcoxon test should be run on time windows rather than (as is the case now) on tf bins calculated by smoothing across multiple time windows. For this, the baseline period's covariance needs a separate calculation outside of the main loop across all tf bins.
* [ ] Using subtract_evoked doesn't seem to change much. Why? Look into this.
* [ ] The wilcoxon code returns only negative z-scores irrespective of how it is invoked. I would expect both negative and positive z-scores. Compare the scipy implementation with NUTMEG code.
* [ ] Should the modified wilcoxon code be placed somewhere else? Any procedures for approaching switching back to using scipy code once the changes are (hopefully) implemented upstream?
* [ ] Submit PR to scipy implementing returning z-scores in scipy.stats.wilcoxon.

[1]: Dalal et al. Simultaneous MEG and intracranial EEG recordings during attentive reading. NeuroImage (2009) vol. 45 (4) pp. 1289-1304