[FieldTrip] Non-parametric granger causality stripe artifact

[Schoffelen, J.M. (Jan Mathijs)]

Hi Walter,

Here are some thoughts:


The zigzag is a well-known convergence issue that may be caused by too large collinearity between your signals. This may occur if:
- the multivariate signal is too multi, (i.e. there are many channels that go into the non-parametric spectral factorization at once) -> possible solution use cfg.sfmethod = ‘bivariate’;
- there is a large common (instantaneous) component underlying your signals. -> not much to do about this, apart perhaps in your case use local bipolar rereference, rather than a common reference

Also, zigzags are more likely to occur if the factorization is applied to time-frequency data, as computed with mtmconvol, which probably has to do with spectral leakage -> if you want time resolved GC, run the shifting time window in an outer loop, and use mtmfft for the spectral transformation. One possible reason is that users might be a bit too creative in their specification of the options for mtmconvol which are suboptimal for the spectral factorization algorithm.

Moreover, the spectra should ideally contain all frequencies from DC up to (close to Nyquist) with a sufficiently high sampling -> use cfg.pad that is relatively large, e.g. 4
Moreover, the spectra should be relatively smooth, without ’spikes’ -> ensure that your spectra don’t contain spikes, and use a sufficiently large number of trials + ideally some multitapering.

Finally, there's an option cfg.granger.stabilityfix (in ft_connectivity_csd2transfer) which might alleviate the zigzags a bit, although I would only use this if the zigzag is minor.

I quickly browsed the archive of this discussion list, and I noticed that optimal analysis strategies for spectrally resolved GC have been discussed before, however they don’t seem to have made it into documentation  on the fieldtrip website. If you feel like it, it would be great if you could make a suggestion for a FAQ that collects the above mentioned (and other tips and tricks), for future reference.

Best wishes,
Jan-Mathijs


On 18 Jul 2023, at 11:16, Canedo Riedel, Walter via fieldtrip <fieldtrip at science.ru.nl<mailto:fieldtrip at science.ru.nl>> wrote:


Dear community,

My name is Walter Cañedo Riedel and I am working at Kelschlab (Mainz University). I am currently analyzing LFP data; there are several channels and I am trying to obtain their Granger causality spectra (by pairs of channels) with the non-parametric approach starting with a wavelet transform. My LFP data is at a resolution of 1000Hz.


I would appreciate if you help me understand how to run the non-parametric Granger causality correctly.  The GC spectra that I obtain have all the following artifact: every second frequency is exalted with respect to its immediate neighbor frequencies, forming a zebra-stripe pattern (image: GC from a channel to another).

I made a simple simulation to look into it. Channel 1 and channel 2 are modulated by different frequencies plus noise. Ch1 from time t-3 to t-1 influences Ch2 linearly at time t.
If Ch1 is modulated by a frequency below 10 Hz, then I get the stripe artifact for its Granger causality to Ch2. (Image: GC has an artifact...)  The script for the simulation is found below after my signature.


Things that I've tried to no avail:


  *   Doing the analysis up to the Nyquist limit (by down-sampling the data to 200Hz and analyzing frequencies 1:100).
  *   High-pass filtering my data (I still get stripes even after filtering out frequencies below 15Hz).
  *   Trying out 'mtmconvol' with hanning or with dpss instead of 'wavelet' for ft_freqanalysis cfg.method.
  *   Trying cfg.output = 'fourier' or 'powandcsd'.
  *   Varying the size of the width for the wavelet or of the respective window for mtmconvol.


<pastedImage.png>

<pastedImage.png>


All the best,
Walter Cañedo Riedel

Central Institute of Mental Health (ZI)
AG Kelsch
+49 (0)621 1703 6203



%% Simulation
clear
lfpchannels = NaN(1200000, 2); % Two channels (CH1 & CH2) with million two hundred thousand datapoints
frequency = 1000; % I have a 1000Hz resolution

SaverGC = [];
SaverPW = [];
for fr = 1:10
    lfpchannels(:, :) = rand(size(lfpchannels)); % Initial noise in the channels

    % Modulate CH1 at (fr)Hz
    dt = 1/frequency;
    StopTime = size(lfpchannels,1)/frequency;
    t = (0:dt:StopTime-dt)';
    F = fr; % Sine wave frequency (Hz)
    Modulator1 = sin(2*pi*F*t) / 10;
    lfpchannels(1:numel(Modulator1), 1) = lfpchannels(1:numel(Modulator1), 1) + Modulator1;
    % Modulate CH2 at 38Hz
    F = 38; % Sine wave frequency (Hz)
    Modulator = sin(2*pi*F*t) / 10;
    lfpchannels(1:numel(Modulator), 2) = lfpchannels(1:numel(Modulator), 2) + Modulator;

    % Random weights for CH1(t-3:t-1) influencing CH2(t).
    FRA_s = 3;
    FRA = rand(FRA_s, 1);
    % Make CH1 influence CH2
    for p = FRA_s+1:size(lfpchannels,1)
        lfpchannels(p, 2) = lfpchannels(p, 2) + mean(FRA.*lfpchannels(p-FRA_s : p-1, 1));
    end
    % z score channels to have same variance
    lfpchannels = (lfpchannels - mean(lfpchannels,1)) ./ std(lfpchannels, 0, 1);

    %%% Make data structure with trials
    TrialMarker = 1:numel(-1.999:0.001:14)/2:(numel(-1.999:0.001:14)/2)*150 + 1;
    data            = struct;
    data.label      = {'1';'2'};
    data.fsample    = 1000;
    data.trial      = {};
    data.time       = {};
    for tr = 1:50
        data.trial{tr} = lfpchannels(TrialMarker(tr):TrialMarker(tr+2)-1, :)';
        data.time{tr} = -1.999:0.001:14;
    end
    data.sampleinfo = [1:50;101:150];
    data.trialinfo = ones(1, 50);

    %%% Run wavelet transform
    cfg        = [];
    cfg.method = 'wavelet';
    cfg.output = 'powandcsd';
    cfg.foi    = 1:100;
    cfg.width  = 7;
    cfg.toi    = (3);
    cfg.pad    = 20;
    freq = ft_freqanalysis(cfg, data);

    cfg           = [];
    cfg.method    = 'granger';
    cfg.granger.sfmethod = 'bivariate';
    Granger = ft_connectivityanalysis(cfg, freq);
    SaverGC(fr, :) = Granger.grangerspctrm(1, :);
    SaverPW(fr, :) = freq.powspctrm(1, :);
end

figure
tiledlayout(1, 2)
nexttile
imagesc(SaverPW)
xlabel("Power spectrum (Hz)")
ylabel("Modulating signal (Hz)")
nexttile
imagesc(SaverGC)
xlabel("GC spectrum (Hz)")
sgtitle("GC has an artifact for the channel modulated by a signal of less than 10Hz")



_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://urldefense.com/v3/__https://doi.org/10.1371/journal.pcbi.1002202__;!!HJOPV4FYYWzcc1jazlU!5M0VHKF356oxF_mZnzCwH0xUkXKCJJl_TBxDQx8L83BkI8GevSMrDsHx0pYTZrZeWb-WkDYPWyup4cbGf3_X2_GuxEMFAJDeWT6ZNQ$

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20230808/975b1fa0/attachment.htm>