[FieldTrip] Questions about using ft_freqstatistics with 'indepsamplesZcoh'

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

Hi Francesco,

In summary, you would like to do statistical inference using a cluster-based permutation test, using (delta-Z)coherence as a test statistic, computed between a fixed reference channel (hopefully not an EEG channel) and all other channels.

With the current state of the fieldtrip code, this will not work out-of-the-box.
The functionality that you require has not been well developed in FieldTrip, see for instance this thread: https://mailman.science.ru.nl/pipermail/fieldtrip/2014-September/021278.html<https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmailman.science.ru.nl%2Fpipermail%2Ffieldtrip%2F2014-September%2F021278.html&data=05%7C02%7Cfieldtrip%40science.ru.nl%7C655b34ee04f344990d3b08ddd4cea6bf%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900706234576769%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=7Qiqa1X8k4Kjs3ciq7jAZ4fqE0s%2FRrEWq086JBhrME0%3D&reserved=0>

Given the lack of interest in general for this specific way of doing statistical inference, there has been little incentive to ensure that the data bookkeeping is consistent (i.e. reshaping vectors back into appropriately sized matrices etc). An added complication in your case, is the fact that the black-box operation that currently converts the input data into a mass-bivariate test-statistic (ft_statfun_indepsamplesZcoh) blows up the dimensionality (going from Nchan to Nchan x Nchan), which indeed makes it challenging to come up with a good way to interpret (and create) clusters.

Part of the challenges would go away, if you were to create a version of ft_statfun_indepsampesZcoh that allows for an additional input argument that lists the index of the channel that is to be used as a reference channel, so that the output of the adjusted ’statfun’ will be reshapeable into a Nchan x Nfreq x Ntime matrix, thus subsequently allowing for ’traditional’ spatial clustering. What would be needed in addition, is some ‘glue code’ in the higher level functions: ft_statistics_montecarlo/analytic/etc + possibly ft_freqstatistics to translate a meaningfully named cfg option (e.g.: cfg.refchannel, which as a default value could be defined as ‘all’) into a vector/scalar of indices that is then passed as input into ft_statfun_indepsamplesZcoh)

I think I tried something like this once in the past, but I don’t remember why this did not make it into FT.

See below (with no guarantee that it works, nor that it - if it works - the results are correct) for some inspiration. Note that the example code probably will not treat data with more than a singleton time point in a correct way.

If you decide to pursue this, please feel free to submit an improved version of the code as a PR on our github repo.

Best wishes,
Jan-Mathijs



===============================

function [s,cfg] = statfun_indepsamplesZcoh(cfg, dat, design);

% STATFUN_indepsamplesZcoh calculates the independent samples coherence Z-statistic
% on the biological data in dat (the dependent variable), using the information on
% the independent variable (iv) in design.
%
% The external interface of this function has to be
%   [s,cfg] = statfun_indepsamplesT(cfg, dat, design);
% where
%   dat    contains the biological data, Nsamples x Nreplications
%          dat must contain fourier representations.
%   design contains the independent variable (iv), Nfac x Nreplications
%
% Configuration options:
%   cfg.computestat    = 'yes' or 'no', calculate the statistic (default='yes')
%   cfg.computecritval = 'yes' or 'no', calculate the critical values of the test statistics (default='no')
%   cfg.computepval    = 'yes' or 'no', calculate the p-values (default='no')
%
% The following options are relevant if cfg.computecritval='yes' and/or
% cfg.computepval='yes'.
%   cfg.alpha = critical alpha-level of the statistical test (default=0.05)
%   cfg.tail = -1, 0, or 1, left, two-sided, or right (default=1)
%              cfg.tail in combination with cfg.computecritval='yes'
%              determines whether the critical value is computed at
%              quantile cfg.alpha (with cfg.tail=-1), at quantiles
%              cfg.alpha/2 and (1-cfg.alpha/2) (with cfg.tail=0), or at
%              quantile (1-cfg.alpha) (with cfg.tail=1).
%

% Copyright (C) 2006, Eric Maris
%
% $Log: statfun_indepsamplesZcoh.m,v $
% Revision 1.3  2006/05/17 11:59:55  erimar
% Corrected bugs after extensive checking of the properties of this
% statfun.
%
% Revision 1.2  2006/05/12 15:32:40  erimar
% Added functionality to calculate one- and two-sided critical values and
% p-values.
%
% Revision 1.1  2006/05/05 13:06:47  erimar
% First version of statfun_indepsamplesZcoh.
%

% set defaults
if ~isfield(cfg, 'computestat'),       cfg.computestat    ='yes';  end;
if ~isfield(cfg, 'computecritval'),    cfg.computecritval ='no';   end;
if ~isfield(cfg, 'computepval'),       cfg.computepval    ='no';   end;
if ~isfield(cfg, 'alpha'),             cfg.alpha          =0.05;   end;
if ~isfield(cfg, 'tail'),              cfg.tail           =1;      end;
if ~isfield(cfg, 'ivar'),              cfg.ivar           =1;      end;

% perform some checks on the configuration
if strcmp(cfg.computepval,'yes') & strcmp(cfg.computestat,'no')
    error('P-values can only be calculated if the test statistics are calculated.');
end;

% original data
nsgn = cfg.dim(1);
if length(cfg.dim)>1,
  nfrq = cfg.dim(2);
else
  nfrq = 1;
end
if length(cfg.dim)>2,
  ntim = cfg.dim(3);
else
  ntim = 1;
end

if ~isfield(cfg, 'cmbindx'),
  cmb1 = 1:nsgn;
  cmb2 = 1:nsgn;
  cmbindx = reshape(1:(length(cmb1)*length(cmb2)),[length(cmb1) length(cmb2)]);
  cmbindx = tril(cmbindx, -1);
  cmbindx = cmbindx(:);
  cmbindx(find(cmbindx==0)) = [];
  ncmb    = length(cmbindx);
else
  cmb1 = unique(cfg.cmbindx(:,1));
  cmb2 = unique(cfg.cmbindx(:,2));
  ncmb = size(cfg.cmbindx,1);
  tmpcmbindx = reshape(1:(length(cmb1)*length(cmb2)),[length(cmb1) length(cmb2)]);
  for kk = 1:ncmb
    cmbindx(kk, 1) = tmpcmbindx(find(cmb1==cfg.cmbindx(kk,1)), ...
                                find(cmb2==cfg.cmbindx(kk,2)));
  end
end

% perform some checks on the design
selc1   = find(design(cfg.ivar,:)==1);
selc2   = find(design(cfg.ivar,:)==2);
nreplc1 = length(selc1);
nreplc2 = length(selc2);
nrepl   = nreplc1 + nreplc2;
if nrepl<size(design,2)
  error('Invalid specification of the independent variable in the design array.');
end;
if nreplc1<2 | nreplc2<2
    error('Every condition must contain at least two trials/tapers.');
end;
dfc1 = nreplc1*2;
dfc2 = nreplc2*2;

if strcmp(cfg.computestat, 'yes')
  tmpstat = zeros(ncmb, nfrq);
  for j = 1:nfrq
    tmpdat = dat([j-1]*nsgn+1:j*nsgn, :);
    % compute the statistic

    %---condition 1
    csdc1    = tmpdat(cmb1,selc1)*tmpdat(cmb2,selc1)'/nreplc1;
    powerc11 = mean(abs(tmpdat(cmb1,selc1)).^2,2);
    powerc12 = mean(abs(tmpdat(cmb2,selc1)).^2,2)';
    csdc1    = csdc1./sqrt(powerc11*powerc12);

    %---condition 2
    csdc2    = tmpdat(cmb1,selc2)*tmpdat(cmb2,selc2)'/nreplc2;
    powerc21 = mean(abs(tmpdat(cmb1,selc2)).^2,2);
    powerc22 = mean(abs(tmpdat(cmb2,selc2)).^2,2)';
    csdc2    = csdc2./sqrt(powerc21*powerc22);

    %---bias-correction
    biasc1   = 1./(dfc1-2);
    biasc2   = 1./(dfc2-2);
    denomZ   = sqrt(1./(dfc1-2) + 1./(dfc2-2));
    tmp      = (atanh(abs(csdc1))-biasc1-atanh(abs(csdc2))+biasc2)./denomZ;
    tmp      = tmp(:);
    tmpstat(:,j) = tmp(cmbindx);
  end
  s.stat = tmpstat;
end

if strcmp(cfg.computecritval,'yes')
  % also compute the critical values
  if cfg.tail==-1
    s.critval = norminv(cfg.alpha);
  elseif  cfg.tail==0
    s.critval = [norminv(cfg.alpha/2),norminv(1-cfg.alpha/2)];
  elseif cfg.tail==1
    s.critval = norminv(1-cfg.alpha);
  end;
end

if strcmp(cfg.computepval,'yes')
  % also compute the p-values
  if cfg.tail==-1
    s.pval = normcdf(s.stat);
  elseif  cfg.tail==0
    s.pval = 2*normcdf(-abs(s.stat));
  elseif cfg.tail==1
    s.pval = 1-normcdf(s.stat);
  end;
end


===============================



On 18 Jul 2025, at 15:23, Francesco Torricelli via fieldtrip <fieldtrip at science.ru.nl> wrote:

Dear FieldTrip Community,

We would like to test differences in coherence at the single-subject level (i.e., with trial as unit-of-observation) with cluster-based permutation test. To do this, we thought we could use the function ft_freqstatistics combined with statfun_indepsamplesZcoh. However, when running the code, we get an error, which is probably due to the fact that statfun_indepsamplesZcoh works on channels combinations. It is therefore non-obvious how to form clusters in the spatial dimension. Since we are interested in computing coherence between just one channel j and all the other channels N-j, we would like to form clusters in the spatial dimension across the N-j channels (plus the frequency and time dimensions). Has anyone in this community ever dealt with this kind of issue before? And, if so, how did you manage to solve the problem?
Here is a more detailed description of how the input data looks like (data1 and data2, one for each condition):
<structure.png>
Both data1 and data2 structures are the output of ft_freqanalysis.
What follows is the code we used:
cfg                  = [];
cfg.method           = 'montecarlo';
cfg.parameter        = 'fourierspctrm';
cfg.statistic        = 'indepsamplesZcoh';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan        = 2;
cfg.neighbours       = neighbours;
cfg.tail             = 0;
cfg.clustertail      = 0;
cfg.alpha            = 0.025;
cfg.numrandomization = 100;
n_M1_categ1          = size(data1.fourierspctrm, 1);
n_M1_categ2          = size(data2.fourierspctrm, 1);
cfg.design           = [ones(1,n_M1_categ1), ones(1,n_M1_categ2)*2]';
cfg.ivar             = 1;

[stat] = ft_freqstatistics(cfg, data1, data2);
It may be helpful to report also that, when running the same code without performing cluster-based correction (i.e., cfg.correctm = 'no'), the stat output that we obtain is a row vector with length equal to the number of channels pairs (reported in stat.cfg.chancmbindx) multiplied by the frequency dimension and the time dimension. This output shape does not correspond to the dimord field reported in the stat structure (i.e., 'chancmb_freq_time').
Thank you very much for your attention.
Best,

Francesco Torricelli


--
Francesco Torricelli, MSc, PhD
Section of Physiology,
Department of Neuroscience and Rehabilitation,
University of Ferrara
Via Fossato di Mortara 19, 44121 Ferrara (Italy)
Tel: +39 347 980 9976
_______________________________________________
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/20250806/66f58492/attachment.htm>