[FieldTrip] Questions about using ft_freqstatistics with 'indepsamplesZcoh'
Francesco Torricelli
francesco.torricelli at unife.it
Wed Aug 6 16:08:21 CEST 2025
Dear Jan Mathijs,
Thank you very much for your detailed reply and for providing a useful code
excerpt as a good starting point.
We'll most likely try to pursue this: in the hopeful case of a positive
result, I'll submit a version of the code as a PR on the FT github repo, as
you suggested.
Best,
Francesco Torricelli
Il giorno mer 6 ago 2025 alle ore 12:27 Schoffelen, J.M. (Jan Mathijs) via
fieldtrip <fieldtrip at science.ru.nl> ha scritto:
> 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://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%7Cbd10aea227f74f35bdac08ddd4f2be54%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900861216110175%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=SFP3aEssXIgXE2l2EC9TZ%2Bq6Ue0IqQBPfuobKnLbdTM%3D&reserved=0
> <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%7Cbd10aea227f74f35bdac08ddd4f2be54%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900861216130963%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=%2BS2IMRQQXf7ii7x94Iic5%2BClk03%2FYyXapWlphYY%2FlD4%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://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmailman.science.ru.nl%2Fmailman%2Flistinfo%2Ffieldtrip&data=05%7C02%7Cfieldtrip%40science.ru.nl%7Cbd10aea227f74f35bdac08ddd4f2be54%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900861216149947%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=M1OhgBJmEGa8A34zxm0nxwEQd7Jc78SzDgR6QHYaij0%3D&reserved=0
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1371%2Fjournal.pcbi.1002202&data=05%7C02%7Cfieldtrip%40science.ru.nl%7Cbd10aea227f74f35bdac08ddd4f2be54%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900861216169879%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=VC3Dvs0p02sqyhmmqs7%2FOQw7aKkiZFiZmfygEAGD3zI%3D&reserved=0
>
>
> _______________________________________________
> fieldtrip mailing list
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmailman.science.ru.nl%2Fmailman%2Flistinfo%2Ffieldtrip&data=05%7C02%7Cfieldtrip%40science.ru.nl%7Cbd10aea227f74f35bdac08ddd4f2be54%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900861216190655%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=zZKRsY21ZIK%2BEW0yEN5Co22al8s1Vv8zSWM7t3splx4%3D&reserved=0
> https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1371%2Fjournal.pcbi.1002202&data=05%7C02%7Cfieldtrip%40science.ru.nl%7Cbd10aea227f74f35bdac08ddd4f2be54%7C084578d9400d4a5aa7c7e76ca47af400%7C1%7C0%7C638900861216205881%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=WOHJbAEuimYXaQYhhKsHq8ECussYwgMqA7pdrwtKZnQ%3D&reserved=0
>
--
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20250806/591ceebf/attachment.htm>
More information about the fieldtrip
mailing list