[FieldTrip] directionality for phase-slope index (PSI)

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

Dear Jed,

Sorry about the confusion, and well-spotted. After thinking about it for a while, I wasn’t sure myself anymore (although I’d have assumed that the documentation in the function is incorrect). Therefore, I decided to simulated some data (with a known delay), in order to evaluate the results from the connectivity computation.
The code is pasted below, but I would appreciate it if this were made into a piece of documentation on the website (hint). Also, based on the outcome, we should update the help documentation in the relevant function (in this case: ft_connectivity_psi).

If you run this piece of code, you’ll see that the (2,1) combination (both when represented as a ‘full’ 2x2xnfreq, and as a ‘sparse’ 1xnfreq matrix, as indexed with {‘chan02’ ‘chan01’}) will have a positive psi-peak. Since by construction the second channel is leading in time, I’d say that a positive psi (as outputted by fieldtrip) means ‘leading’ in time. This is confirmed by the granger causality analysis.

Could I ask you (and anyone else interested) to verify this, and if you agree, to make a suggestion for an update of the help documentation in ft_connectivity_psi, as well as a draft for a FAQ on the website? Both should be submitted as pull requests on https://github.com/fieldtrip/fieldtrip (for code changes) and https://github.com/fieldtrip/website (for website changes), respectively.

Thanks and with best wishes,

Jan-Mathijs



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE BELOW CODE SHOULD IDEALLY GO INTO A FAQ OR EXAMPLE SCRIPT ON THE WEBSITE

% simulate some data
dat = ft_preproc_bandpassfilter(randn(1,100050),1000,[15 25],[],'firws');

dat1 = reshape(dat(1:100000),[],100)';
dat2 = reshape(dat(51:end),[],100)'; % 2 is leading 1

tlck.trial = zeros(100,2,1000);
tlck.trial(:,1,:) = dat1+randn(100,1000)./4;
tlck.trial(:,2,:) = dat2+randn(100,1000)./4;
tlck.label  = {'chan01';'chan02'};
tlck.dimord = 'rpt_chan_time';
tlck.time   = (0:999)./1000;

data = ft_checkdata(tlck,'datatype', 'raw');

% spectral decomposition
cfg        = [];
cfg.method = 'mtmfft';
%cfg.foilim = [0 100];
cfg.tapsmofrq = 2;
cfg.output = 'fourier';
cfg.pad    = 4;
freq       = ft_freqanalysis(cfg, data);


cfg.output = 'powandcsd';
cfg.channelcmb = {'all' 'all'};
freq_csd   = ft_freqanalysis(cfg, data);


% connectivity estimation
cfg           = [];
cfg.method    = 'psi';
cfg.bandwidth = 5;
psi           = ft_connectivityanalysis(cfg, freq);
psi2          = ft_connectivityanalysis(cfg, freq_csd);

% bonus: estimate delay from the phase difference spectrum
cfg           = [];
cfg.method    = 'coh';
cfg.complex   = 'complex';
coh           = ft_connectivityanalysis(cfg, freq);

findx = nearest(freq.freq,18):nearest(freq.freq,22);
X     = [ones(1,numel(findx));freq.freq(findx)-mean(freq.freq(findx))];
b     = unwrap(squeeze(angle(coh.cohspctrm(2,1,findx))))'/X;
delay = 1000.*b(2)./(2.*pi);

% another bonus: estimate granger causality
cfg        = [];
cfg.method = 'granger';
granger    = ft_connectivityanalysis(cfg, freq);
granger2   = ft_connectivityanalysis(cfg, freq_csd);

cfg = [];
cfg.xlim = [0 50];
cfg.parameter = 'psispctrm';
ft_connectivityplot(cfg, psi);


On 29 Apr 2019, at 23:27, Jed Meltzer <jedmeltzer at yahoo.com<mailto:jedmeltzer at yahoo.com>> wrote:

Hello,
I'm Jed Meltzer, a researcher at the Rotman Research Institute of Baycrest Hospital in Toronto. A student in my lab is working on an MEG project involving connectivity between the left and right motor cortex during hand movement, and we are interested in quantifying causal influences between the two regions using methods such as Granger Causality and phase-slope index (PSI). When calculating PSI in field-trip, we are a bit confused about which direction is which, and are seeking clarification. Here is our cfg setup:

data.label   = {'LM1' 'RM1' 'LV1' 'RV1'};
data.trial   = tdata;
data.time    = tdata_time;
data.cfg     = cfg;

%% time-frequency analysis

% tf
cfg = [];
cfg.output = 'fourier';
cfg.method = 'mtmconvol';
cfg.keeptrials = 'yes';
cfg.taper = 'hanning'; %
cfg.toi = timesec;
cfg.channel = 'all'
cfg.foi = 1:0.5:80;
cfg.t_ftimwin = ones(length(cfg.foi),1).*0.5;
cfg.pad = 10 % this is the parameter to output evenly spaced freq bins

freq = ft_freqanalysis(cfg, data)

freqdim = freq.freq;
timedim = freq.time;
                                        % psi
cfg = [];
cfg.method = 'psi'; %
cfg.bandwidth = 5; %
stat = ft_connectivityanalysis(cfg, freq);
RM1toLM1_psi = squeeze(stat.psispctrm(1,2,:,:));
LM1toRM1_psi = squeeze(stat.psispctrm(2,1,:,:));


We are not sure about those last two lines of the code. We would like to have the causal influence of channel 1 to channel 2 be reflected as a positive number. We have read through the documentation and we find two parts of it to seem contradictory. First, there is this page:

http://www.fieldtriptoolbox.org/faq/in_what_way_can_frequency_domain_data_be_represented_in_fieldtrip/

It states:
------

Note that this representation lacks a ‘labelcmb’ field, and that the ‘dimord’ is ‘chan_chan_freq’. This means that the numeric data now implicitly contains both the combinations {‘a’ ‘b’} (in coh.cohspctrm(1,2,:) ) , and the combination {‘b’ ‘a’} (in coh.cohspctrm(2,1,:) ). For a quantity like the coherence spectrum the values across the diagonal are symmetric, but for complex-valued quantities, as well as for directional measures of interaction, the values at the entries across the diagonal are typically different. The convention used by FieldTrip is that the row-channel ‘causes’ the column-channel.

--------
So, we interpret this to mean that if channel 1 is LM1, and channel 2 is RM1, then we would probably want:
LM1toRM1_psi = squeeze(stat.psispctrm(1,2,:,:));
because the first dimension is the "row" channel and the second dimension is the "column" channel.

However, in the help for the psi function:
http://www.fieldtriptoolbox.org/reference/ft_connectivity_psi/
it says:
-----

If the phase slope
 index is positive, then the first chan (1st dim) becomes more lagged (or
 less leading) with higher frequency, indicating that it is causally
 driven by the second channel (2nd dim)

------
So following this logic, if we want a positive number to reflect a causal influence from LM1 (channel 1) to RM1 (channel 2), then we would instead want to write:
LM1toRM1_psi = squeeze(stat.psispctrm(2,1,:,:));

Can you help me figure out which one is correct? Our data have periods of both positive and negative PSI (quite interestingly) so it's impossible for us to "guess" which one is correct just from the pattern.

Thanks,
Jed
_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://doi.org/10.1371/journal.pcbi.1002202

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20190501/ed2c642d/attachment-0001.html>
-------------- next part --------------
_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://doi.org/10.1371/journal.pcbi.1002202