eloreta does not behave as expected when using spectrally resolved data

[arnodelorme]

Environment and versions (please complete the following information):

• Operating System: Osx Mojave
• MATLAB version 2020a
• FieldTrip version 20200607

Considering the small snippet of code

load('dataPre_and_vol.mat');

% Compute spectrum
cfg                  = [];
cfg.output  = 'powandcsd';
cfg.channel = 'all';
cfg.method  = 'mtmfft';
cfg.taper   = 'boxcar';
dataFreq = ft_freqanalysis(cfg, dataPre);

% source reconstruction
cfg             = [];
cfg.method      = 'eloreta';
cfg.sourcemodel = sourcemodel;
cfg.headmodel   = vol.vol;
sourceFreq      = ft_sourceanalysis(cfg, dataFreq);  % compute the source model

% -----------------------------

% Compute an ERP
cfg                  = [];
cfg.covariance       = 'yes';
cfg.covariancewindow = [-1 0]; % calculate the average of the covariance matrices 
                                   % for each trial (but using the pre-event baseline  data only)
dataAvg = ft_timelockanalysis(cfg, dataPre);

% source reconstruction
cfg             = [];
cfg.method      = 'eloreta';
cfg.sourcemodel = sourcemodel;
cfg.headmodel   = vol.vol;
sourceTime      = ft_sourceanalysis(cfg, dataAvg);  % compute the source model

We have sourceFreq

>> sourceFreq = 

  struct with fields:

      freq: 32
       cfg: [1×1 struct]
       dim: [15 18 15]
    inside: [4050×1 logical]
       pos: [4050×3 double]
    method: 'average'
       avg: [1×1 struct]

>> sourceFreq.avg.mom(3938)

ans =

  1×1 cell array

    {3×31 double}

Now for time, we have

>> sourceTime = 

  struct with fields:

      time: [1×384 double]
       cfg: [1×1 struct]
       dim: [15 18 15]
    inside: [4050×1 logical]
       pos: [4050×3 double]
    method: 'average'
       avg: [1×1 struct]

>> sourceTime.avg.mom(3938)

ans =

  1×1 cell array

    {3×384 double}

One would expect sourceFreq.freq to be an array of frequencies of length 193 as it is the case for sourceTime.time which contains an array of latencies. For the moment for frequency source reconstruction, one would expect sourceFreq.avg.mom to be of size {3x193} since there are 193 frequencies. The size {3x31} above where 31 is the number of channels and makes little sense to me.

[arnodelorme]

dataPre_and_vol.mat.zip

[schoffelen]

I think that ft_sourceanalysis does not work with multiple frequency bins, and will probably (silently, and perhaps unexpectedly) average across frequencies in some way.

In addition, in order to be able to get some frequency resolved stuff in the 'mom' field, the input should at least contain fourierspctrm, rather than crsspctrm/powspctrm. This is admittedly not documented well, but I guess nobody actually ever really cared to try out spectrally resolved data in combination with eloreta.

Could you try out running ft_freqanalysis with cfg.output = 'fourier', and use this as input into ft_freqanalysis, and report back what came out of this ?

[arnodelorme]

I thought the cross-spectrum could be used by Loreta as the covariance matrix is used in the time domain, which is why I used 'powandcsd'. I have also tried using 'pow' only as output, and the ft_sourceanalysis function crashes (Conversion to double from cell is not possible) which is consistent with my intuition that the Loreta function uses cross spectral density.

WIth 'fourier', now I get

>> sourceFreq.avg.mom(3938)

ans =

  1×1 cell array

    {3×80 double}

For 80 trials (not frequencies). In my case, trials make little sense, because this is just segmented continuous data (to perform Welch method as indicated in this tutorial http://www.fieldtriptoolbox.org/workshop/madrid2019/tutorial_freq/). I have tried to set the option not to keep trials and got the error

Error using ft_freqanalysis (line 361)
cfg.output = 'fourier' requires cfg.keeptrials = 'yes' and cfg.keeptapers = 'yes'

(which makes sense I think; you cannot combine the complex Fourier transform from the multiple trials). BTW, the dimord for this is wrong though indicating "rpttap_chan_freq" (repeated tapers while these are trials).

I have tried transposing the frequency data, but it did not work either. Finally, I was able to get to a solution which is not ideal but seems to work. I used the Fourier method, then swapped the trial data and the frequency data (in the frequency domain, my guess is that it should be OK to average across trials instead of frequencies when performing source localization but I might be missing something).

cfg                  = [];
cfg.output  = 'fourier';
%cfg.taper   = 'hanning';
cfg.channel = 'all';
cfg.method  = 'mtmfft';
cfg.taper   = 'boxcar';
dataFreq = ft_freqanalysis(cfg, dataPre);

% swap trials and frequencies for source localization
dataFreq.freq = [1:80];
dataFreq.fourierspctrm = permute(dataFreq.fourierspctrm, [3 2 1]);
dataFreq.cumsumcnt = ones(1, size(dataFreq.fourierspctrm,1))*80;
dataFreq.cumtapcnt = ones(1, size(dataFreq.fourierspctrm,1));
dataFreq = rmfield(dataFreq, 'trialinfo');

% source reconstruction
cfg             = [];
cfg.method      = 'eloreta';
cfg.sourcemodel = sourcemodel;
cfg.headmodel   = vol.vol;
sourceFreq      = ft_sourceanalysis(cfg, dataFreq);  % compute the source model

Then I get moments of the desired size (I have 193 frequencies in my spectral decomposition).

>> sourceFreq.avg.mom(3938)

ans =

  1×1 cell array

    {3×193 double}

Comments appreciated.

[schoffelen]

The bottom line here seems to be that, despite the suggestions made by the description in the function's help, ft_sourceanalysis with cfg.method = 'eloreta' does not behave as per the user's expectation, when the input data contains spectrally resolved data.

• [ ] ft_sourceanalysis does not work with multiple frequency bins in general, since the original 'dics' and 'pcc' methods per design don't support this.
• [ ] personally I don't give a jot about eloreta, but it would be nice if indeed the code does what it kind of promises.
• [ ] the most easy fix would be to adjust the documentation.
• [ ] right now a frequency resolved eloreta analysis in principle would be possible running a for loop around ft_sourceanalysis, with cfg.frequency as the variable that changes in each iteration.
• [ ] there is a particular expectation that the output 'mom' field should contain a power spectrum, which is incorrect. If anything, the power should go in to the power field.
• [ ] the low-level eloreta code takes as an input a square matrix, no idea how this should work for multiple frequency bins, i.e. should a frequency-specific spatial filter be computed, or should a 'generic' spatial filter be applied to the individual frequency bins?

[arnodelorme]

right now a frequency resolved eloreta analysis in principle would be possible running a for loop around ft_sourceanalysis, with cfg.frequency as the variable that changes in each iteration.

I have tried that solution and it return garbage data. Also the moment array is of size 3 x chans which makes no sense. I would definitely not advise that solution.

eLoreta is what Stefan Haufe recommends for connectivity analysis. I do not think it is as outdated as you make it sound. Thanks for taking the time.

[schoffelen]

To reiterate: the mom is not a field in which power can be represented, since it is the output of the left multiplication of the dipole specific spatial filter with a sensor data matrix, so it should reflect either a dipole moment as in a time course of an event-related field/potential, or fourier coefficients (which need to be absoluted-squared-and-averaged-across-trials before it represents power). Apparently, the person who contributed the eloreta code, thought that in absence of a fourier-matrix, the csd matrix could be used as a surrogate, which results in a 3xnchan matrix in the output, which indeed is bogus. We look forward to a PR that fixes this.

So, in short, if the garbage you refer to in the case of a for-loop around ft_sourceanalysis across frequencies, I cannot agree more. Yet, I recommend looking at the 'pow' field instead.

[schoffelen]

BTW, I would recommend beamformers for connectivity analysis... Also, for the record, I have never said that erloreta is outdated, I just meant this as a nudge that people who actually care about this specific method are in a better position (and have more incentive possibly) to get it well-supported in software.

[schoffelen]

Note: on line 537 in ft_sourceanalysis the avg variable is specified as Cf, in the absence of a fourierspctrm. This is clearly nonsense: this explains the strangely shaped mom in the output

[schoffelen]

Ah, if you meant with 'garbage' data that the pow is not coming out as expected (in the sense that it is possible to reconstruct a nice spectrum per dipole, which has an interpretable shape) this would mean that a per-frequency spatial filter does not allow to obtain meaningful spectra without some post/ad hoc normalisation. That would mean that the best way would be to use the same spatial filter for each frequency bin. This could for instance be achieved with running ft_sourceanalysis once with cfg.keepfilter = 'yes' (and integrating across the whole frequency range), followed by a for loop across frequencies, while prespecifying cfg.sourcemodel.filter

[schoffelen]

Note that in case of frequency domain input data into sourceanalysis, on line 313 a subselection of the input data is made (as per ft_selectdata), where both the time and the frequency dimensions are force averaged.

[arnodelorme]

Yes I mean that the "pow" output is not coming out as expected when a single frequency is specified when calling ft_sourceanalysis. The output varies based on the frequency parameter but the resulting solution look like a checkerboard when plotted - while it should obviously be smooth. This might be related to the bogus code using "csd" you are mentioning.

Yes I understand that the moments need to be absolute-squared-and-averaged-across-trials before it represents power. One quick fix right now would be to disallow using eLoreta when only power and csd matrices are provided and force the use of the Fourier coefficient. Thanks you Jan for taking the time.

[schoffelen]

No worries Arno, I think your (and my) confusion are reasonable, so it makes sense to sort this out. The only slightly annoying thing is that it's most likely that other people (than you and me) created the mess, and that's up to us to fix it. I'll keep you posted, keep an eye on #1456