Skip to content
Permalink
release
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
function [stat, cfg] = ft_statistics_montecarlo(cfg, dat, design, varargin)
% FT_STATISTICS_MONTECARLO performs a nonparametric statistical test by calculating
% Monte-Carlo estimates of the significance probabilities and/or critical values from
% the permutation distribution. This function should not be called directly, instead
% you should call the function that is associated with the type of data on which you
% want to perform the test.
%
% Use as
% stat = ft_timelockstatistics(cfg, data1, data2, data3, ...)
% stat = ft_freqstatistics (cfg, data1, data2, data3, ...)
% stat = ft_sourcestatistics (cfg, data1, data2, data3, ...)
%
% where the data is obtained from FT_TIMELOCKANALYSIS, FT_FREQANALYSIS or
% FT_SOURCEANALYSIS respectively, or from FT_TIMELOCKGRANDAVERAGE,
% FT_FREQGRANDAVERAGE or FT_SOURCEGRANDAVERAGE respectively
% and with cfg.method = 'montecarlo'
%
% The configuration options that can be specified are:
% cfg.numrandomization = number of randomizations, can be 'all'
% cfg.correctm = string, apply multiple-comparison correction, 'no', 'max', cluster', 'tfce', 'bonferroni', 'holm', 'hochberg', 'fdr' (default = 'no')
% cfg.alpha = number, critical value for rejecting the null-hypothesis per tail (default = 0.05)
% cfg.tail = number, -1, 1 or 0 (default = 0)
% cfg.correcttail = string, correct p-values or alpha-values when doing a two-sided test, 'alpha','prob' or 'no' (default = 'no')
% cfg.ivar = number or list with indices, independent variable(s)
% cfg.uvar = number or list with indices, unit variable(s)
% cfg.wvar = number or list with indices, within-cell variable(s)
% cfg.cvar = number or list with indices, control variable(s)
% cfg.feedback = string, 'gui', 'text', 'textbar' or 'no' (default = 'text')
% cfg.randomseed = string, 'yes', 'no' or a number (default = 'yes')
%
% If you use a cluster-based statistic, you can specify the following options that
% determine how the single-sample or single-voxel statistics will be thresholded and
% combined into one statistical value per cluster.
% cfg.clusterstatistic = how to combine the single samples that belong to a cluster, 'maxsum', 'maxsize', 'wcm' (default = 'maxsum')
% the option 'wcm' refers to 'weighted cluster mass', a statistic that combines cluster size and intensity;
% see Hayasaka & Nichols (2004) NeuroImage for details
% cfg.clusterthreshold = method for single-sample threshold, 'parametric', 'nonparametric_individual', 'nonparametric_common' (default = 'parametric')
% cfg.clusteralpha = for either parametric or nonparametric thresholding per tail (default = 0.05)
% cfg.clustercritval = for parametric thresholding (default is determined by the statfun)
% cfg.clustertail = -1, 1 or 0 (default = 0)
%
% To include the channel dimension for clustering of channel level data, you should specify
% cfg.neighbours = neighbourhood structure, see FT_PREPARE_NEIGHBOURS
% If you specify an empty neighbourhood structure, clustering will only be done
% over frequency and/or time and not over neighbouring channels.
%
% The statistic that is computed for each sample in each random reshuffling
% of the data is specified as
% cfg.statistic = 'indepsamplesT' independent samples T-statistic,
% 'indepsamplesF' independent samples F-statistic,
% 'indepsamplesregrT' independent samples regression coefficient T-statistic,
% 'indepsamplesZcoh' independent samples Z-statistic for coherence,
% 'depsamplesT' dependent samples T-statistic,
% 'depsamplesFmultivariate' dependent samples F-statistic MANOVA,
% 'depsamplesregrT' dependent samples regression coefficient T-statistic,
% 'actvsblT' activation versus baseline T-statistic.
% or you can specify your own low-level statistical function.
%
% You can also use a custom statistic of your choice that is sensitive to the
% expected effect in the data. You can implement the statistic in a "statfun" that
% will be called for each randomization. The requirements on a custom statistical
% function is that the function is called ft_statfun_xxx, and that the function returns
% a structure with a "stat" field containing the single sample statistical values.
% Have a look at the functions in the fieldtrip/statfun directory (e.g.
% FT_STATFUN_INDEPSAMPLEST) for the correct format of the input and output.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS, FT_SOURCESTATISTICS,
% FT_STATISTICS_ANALYTIC, FT_STATISTICS_STATS, FT_STATISTICS_MVPA,
% FT_STATISTICS_CROSSVALIDATE
% Undocumented local options:
% cfg.resampling permutation, bootstrap
% cfg.computecritval yes|no, for the statfun
% cfg.computestat yes|no, for the statfun
% cfg.computeprob yes|no, for the statfun
% cfg.voxelstatistic deprecated
% cfg.voxelthreshold deprecated
% cfg.precondition before|after|[], for the statfun
% Copyright (C) 2005-2015, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% FieldTrip is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$
% do a sanity check on the input data
assert(isnumeric(dat), 'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');
assert(isnumeric(design), 'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');
% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamed', {'factor', 'ivar'});
cfg = ft_checkconfig(cfg, 'renamed', {'unitfactor', 'uvar'});
cfg = ft_checkconfig(cfg, 'renamed', {'repeatedmeasures', 'uvar'});
cfg = ft_checkconfig(cfg, 'renamedval', {'clusterthreshold', 'nonparametric', 'nonparametric_individual'});
cfg = ft_checkconfig(cfg, 'renamedval', {'correctm', 'yes', 'max'});
cfg = ft_checkconfig(cfg, 'renamedval', {'correctm', 'none', 'no'});
cfg = ft_checkconfig(cfg, 'renamedval', {'correctm', 'bonferoni', 'bonferroni'});
cfg = ft_checkconfig(cfg, 'renamedval', {'correctm', 'holms', 'holm'});
cfg = ft_checkconfig(cfg, 'required', {'statistic'});
cfg = ft_checkconfig(cfg, 'forbidden', {'ztransform', 'removemarginalmeans', 'randomfactor', 'voxelthreshold', 'voxelstatistic'});
cfg = ft_checkconfig(cfg, 'renamedval', {'statfun', 'depsamplesF', 'ft_statfun_depsamplesFmultivariate'});
cfg = ft_checkconfig(cfg, 'renamedval', {'statfun', 'ft_statfun_depsamplesF', 'ft_statfun_depsamplesFmultivariate'});
% set the defaults for the main function
cfg.alpha = ft_getopt(cfg, 'alpha', 0.05);
cfg.tail = ft_getopt(cfg, 'tail', 0);
cfg.correctm = ft_getopt(cfg, 'correctm', 'no');
cfg.resampling = ft_getopt(cfg, 'resampling', 'permutation');
cfg.feedback = ft_getopt(cfg, 'feedback', 'text');
cfg.ivar = ft_getopt(cfg, 'ivar', 'all');
cfg.uvar = ft_getopt(cfg, 'uvar', []);
cfg.cvar = ft_getopt(cfg, 'cvar', []);
cfg.wvar = ft_getopt(cfg, 'wvar', []);
cfg.correcttail = ft_getopt(cfg, 'correcttail', 'no');
cfg.precondition = ft_getopt(cfg, 'precondition', []);
% explicit check for option 'yes' in cfg.correctail.
if strcmp(cfg.correcttail, 'yes')
ft_error('cfg.correcttail = ''yes'' is not allowed, use either ''prob'', ''alpha'' or ''no''')
end
if strcmp(cfg.correctm, 'tfce')
% TODO this could require some better defaults
cfg.connectivity = ft_getopt(cfg, 'connectivity', []);
cfg.tfce_h0 = ft_getopt(cfg, 'tfce_h0', 0);
cfg.tfce_H = ft_getopt(cfg, 'tfce_H', 2);
cfg.tfce_E = ft_getopt(cfg, 'tfce_E', 0.5);
cfg.tfce_nsteps = ft_getopt(cfg, 'tfce_nsteps', 100);
else
% these options only apply to tfce, to ensure appropriate configs they are forbidden when _not_ clustering
cfg = ft_checkconfig(cfg, 'unused', {'tfce_h0', 'tfce_H', 'tfce_E', 'tfce_nsteps'});
end
if strcmp(cfg.correctm, 'cluster')
% set the defaults for clustering
cfg.clusterstatistic = ft_getopt(cfg, 'clusterstatistic', 'maxsum');
cfg.clusterthreshold = ft_getopt(cfg, 'clusterthreshold', 'parametric');
cfg.clusteralpha = ft_getopt(cfg, 'clusteralpha', 0.05);
cfg.clustercritval = ft_getopt(cfg, 'clustercritval', []);
cfg.clustertail = ft_getopt(cfg, 'clustertail', cfg.tail);
cfg.connectivity = ft_getopt(cfg, 'connectivity', []); % the default is dealt with below
else
% these options only apply to clustering, to ensure appropriate configs they are forbidden when _not_ clustering
cfg = ft_checkconfig(cfg, 'unused', {'clusterstatistic', 'clusteralpha', 'clustercritval', 'clusterthreshold', 'clustertail'});
end
if any(strcmp(cfg.correctm, {'cluster' 'tfce'}))
% these options might require a spatial neighbourhood matrix
% deal with the neighbourhood of the channels/triangulation/voxels
if isempty(cfg.connectivity)
if isfield(cfg, 'dim') && ~isfield(cfg, 'channel') && ~isfield(cfg, 'tri')
% input data can be reshaped into a 3D volume, use bwlabeln/spm_bwlabel rather than clusterstat
ft_info('using connectivity of voxels in 3-D volume\n');
cfg.connectivity = nan;
elseif isfield(cfg, 'tri')
% input data describes a surface along which neighbours can be defined
ft_info('using connectivity of vertices along triangulated surface\n');
cfg.connectivity = triangle2connectivity(cfg.tri);
if isfield(cfg, 'insideorig')
cfg.connectivity = cfg.connectivity(cfg.insideorig, cfg.insideorig);
end
elseif isfield(cfg, 'avgoverchan') && istrue(cfg.avgoverchan)
% channel dimension has been averaged across, no sense in clustering across space
cfg.connectivity = true(1);
elseif isfield(cfg, 'channel')
cfg.neighbours = ft_getopt(cfg, 'neighbours', []);
cfg.connectivity = channelconnectivity(cfg);
else
% there is no connectivity in the spatial dimension
cfg.connectivity = false(size(dat,1));
end
else
% use the specified connectivity: this is not fully robust because
% there is no guarantee that the order of the spatial elements in the
% data is the same as the order of the spatial elements in the
% adjacency matrix
end
end
% for backward compatibility and other warnings relating correcttail
if isfield(cfg,'correctp') && strcmp(cfg.correctp,'yes')
ft_warning('cfg.correctp has been renamed to cfg.correcttail and the options have been changed')
disp('setting cfg.correcttail to ''prob''')
cfg.correcttail = 'prob';
cfg = rmfield(cfg,'correctp');
elseif isfield(cfg,'correctp') && strcmp(cfg.correctp,'no')
cfg = ft_checkconfig(cfg, 'renamed', {'correctp', 'correcttail'});
end
if strcmp(cfg.correcttail,'no') && cfg.tail==0 && cfg.alpha==0.05
ft_warning('Doing a two-sided test without correcting p-values or alpha-level, p-values and alpha-level will reflect one-sided tests per tail. See http://bit.ly/2YQ1Hm8')
end
% for backward compatibility
if size(design,2)~=size(dat,2)
design = transpose(design);
end
if ischar(cfg.ivar) && strcmp(cfg.ivar, 'all')
cfg.ivar = 1:size(design,1);
end
% fetch function handle to the low-level statistics function
statfun = ft_getuserfun(cfg.statistic, 'statfun');
if isempty(statfun)
ft_error('could not locate the appropriate statistics function');
else
ft_info('using "%s" for the single-sample statistics\n', func2str(statfun));
end
% construct the resampled design matrix or data-shuffling matrix
ft_info('constructing randomized design\n');
resample = resampledesign(cfg, design);
Nrand = size(resample,1);
% most of the statfuns result in this warning, which is not interesting
ws = ft_warning('off', 'MATLAB:warn_r14_stucture_assignment');
if strcmp(cfg.correctm, 'cluster')
% determine the critical value for cluster thresholding
if strcmp(cfg.clusterthreshold, 'nonparametric_individual') || strcmp(cfg.clusterthreshold, 'nonparametric_common')
ft_info('using a nonparametric threshold for clustering\n');
cfg.clustercritval = []; % this will be determined later
elseif strcmp(cfg.clusterthreshold, 'parametric') && isempty(cfg.clustercritval)
ft_info('computing a parametric threshold for clustering\n');
tmpcfg = cfg; % the next line does not pass on non-standard options that a statfun might use
% tmpcfg = keepfields(cfg, {'dim' 'dimord' 'clusteralpha' 'clustertail' 'ivar' 'uvar' 'cvar' 'wvar' 'contrastcoefs'});
tmpcfg.computecritval = 'yes'; % explicitly request the computation of the crtitical value
tmpcfg.computestat = 'no'; % skip the computation of the statistic
tmpcfg.alpha = cfg.clusteralpha; % the statfun uses cfg.alpha most likely
try
cfg.clustercritval = getfield(statfun(tmpcfg, dat, design), 'critval');
catch
disp(lasterr);
ft_error('could not determine the parametric critical value for clustering');
end
elseif strcmp(cfg.clusterthreshold, 'parametric') && ~isempty(cfg.clustercritval)
ft_info('using the specified parametric threshold for clustering\n');
cfg.clusteralpha = [];
end
end
% compute the statistic for the observed data
ft_progress('init', cfg.feedback, 'computing statistic');
% get an estimate of the time required per evaluation of the statfun
time_pre = cputime;
try
% the nargout function in MATLAB 6.5 and older does not work on function handles
num = nargout(statfun);
catch
num = 1;
end
if num==1
% only the statistic is returned
[statobs] = statfun(cfg, dat, design);
elseif num==2
% both the statistic and the (updated) configuration are returned
[statobs, cfg] = statfun(cfg, dat, design);
elseif num==3
% both the statistic and the (updated) configuration and the (updated) data are returned
tmpcfg = cfg;
if strcmp(cfg. precondition, 'before'), tmpcfg.preconditionflag = 1; end
[statobs, tmpcfg, dat] = statfun(tmpcfg, dat, design);
tmpcfg.preconditionflag = 0;
cfg = tmpcfg;
end
if isstruct(statobs)
% remember all details for later reference, continue to work with the statistic
statfull = statobs;
statobs = statobs.stat;
end
% remember the statistic for later reference, continue to work with the statistic
statfull.stat = statobs;
time_eval = cputime - time_pre;
ft_info('estimated time per randomization is %.2f seconds\n', time_eval);
% pre-allocate some memory
if strcmp(cfg.correctm, 'cluster')
statrand = zeros(size(statobs,1), size(resample,1), class(dat)); % this reduces the memory footprint, requires the user to use ft_struct2single on the input data
else
prb_pos = zeros(size(statobs));
prb_neg = zeros(size(statobs));
end
if strcmp(cfg.precondition, 'after')
tmpcfg = cfg;
tmpcfg.preconditionflag = 1;
[tmpstat, tmpcfg, dat] = statfun(tmpcfg, dat, design);
end
if any(strcmp(cfg.correctm, {'tfce' 'max'}))
% pre-allocate the memory to hold the distribution of most extreme positive (right) and negative (left) statistical values
posdistribution = nan(1,Nrand);
negdistribution = nan(1,Nrand);
end
% compute the statistic for the randomized data and count the outliers
for i = 1:Nrand
ft_progress(i/Nrand, 'computing statistic %d from %d\n', i, Nrand);
if strcmp(cfg.resampling, 'permutation')
tmpdesign = design(:,resample(i,:)); % the columns in the design matrix are reshufled by means of permutation
tmpdat = dat; % the data itself is not shuffled
if size(tmpdesign,1)==size(tmpdat,2)
tmpdesign = transpose(tmpdesign);
end
elseif strcmp(cfg.resampling, 'bootstrap')
tmpdesign = design; % the design matrix is not shuffled
tmpdat = dat(:,resample(i,:)); % the columns of the data are resampled by means of bootstrapping
end
if any(strcmp(cfg.correctm, {'cluster' 'tfce'}))
% keep each randomization in memory for cluster postprocessing
dum = statfun(cfg, tmpdat, tmpdesign);
if isstruct(dum)
statrand(:,i) = dum.stat;
else
statrand(:,i) = dum;
end
else
% do not keep each randomization in memory, but process them on the fly
statrand = statfun(cfg, tmpdat, tmpdesign);
if isstruct(statrand)
statrand = statrand.stat;
end
% the following line is for debugging
% stat.statkeep(:,i) = statrand;
if strcmp(cfg.correctm, 'max')
% compare each data element with the maximum statistic
prb_pos = prb_pos + (statobs<max(statrand(:)));
prb_neg = prb_neg + (statobs>min(statrand(:)));
posdistribution(i) = max(statrand(:));
negdistribution(i) = min(statrand(:));
else
% compare each data element with its own statistic
prb_pos = prb_pos + (statobs<statrand);
prb_neg = prb_neg + (statobs>statrand);
end
end
end
ft_progress('close');
if strcmp(cfg.correctm, 'cluster')
% do the cluster postprocessing
[stat, cfg] = clusterstat(cfg, statrand, statobs);
elseif strcmp(cfg.correctm, 'tfce')
[stat, cfg] = tfcestat(cfg, statrand, statobs);
else
if ~isequal(cfg.numrandomization, 'all')
% in case of random permutations (i.e., montecarlo sample, and NOT full
% permutation), the minimum p-value should not be 0, but 1/N
prb_pos = prb_pos + 1;
prb_neg = prb_neg + 1;
Nrand = Nrand + 1;
end
switch cfg.tail
case 1
clear prb_neg % not needed any more, free some memory
stat.prob = prb_pos./Nrand;
case -1
clear prb_pos % not needed any more, free some memory
stat.prob = prb_neg./Nrand;
case 0
% for each observation select the tail that corresponds with the lowest probability
prb_neg = prb_neg./Nrand;
prb_pos = prb_pos./Nrand;
stat.prob = min(prb_neg, prb_pos); % this is the probability for the most unlikely tail
end
end
% In case of a two tailed test, the type-I error rate (alpha) refers to
% both tails of the distribution, whereas the stat.prob value computed sofar
% corresponds with one tail, i.e. the probability, under the assumption of
% no effect or no difference (the null hypothesis), of obtaining a result
% equal to or more extreme than what was actually observed. The decision
% rule whether the null-hopothesis should be rejected given the observed
% probability therefore should consider alpha divided by two, to correspond
% with the probability in one of the tails (the most extreme tail). This
% is conceptually equivalent to performing a Bonferroni correction for the
% two tails.
%
% An alternative solution to distribute the alpha level over both tails is
% achieved by multiplying the probability with a factor of two, prior to
% thresholding it wich cfg.alpha. The advantage of this solution is that
% it results in a p-value that corresponds with a parametric probability.
% Below both options are realized
if strcmp(cfg.correcttail, 'prob') && cfg.tail==0
stat.prob = stat.prob .* 2;
stat.prob(stat.prob>1) = 1; % clip at p=1
% also correct the probabilities in the pos/negcluster fields
if isfield(stat, 'posclusters')
for i=1:length(stat.posclusters)
stat.posclusters(i).prob = stat.posclusters(i).prob*2;
if stat.posclusters(i).prob>1; stat.posclusters(i).prob = 1; end
end
end
if isfield(stat, 'negclusters')
for i=1:length(stat.negclusters)
stat.negclusters(i).prob = stat.negclusters(i).prob*2;
if stat.negclusters(i).prob>1; stat.negclusters(i).prob = 1; end
end
end
elseif strcmp(cfg.correcttail, 'alpha') && cfg.tail==0
cfg.alpha = cfg.alpha / 2;
end
% compute range of confidence interval p ? 1.96(sqrt(var(p))), with var(p) = var(x/n) = p*(1-p)/N
stddev = sqrt(stat.prob.*(1-stat.prob)/Nrand);
stat.cirange = 1.96*stddev;
if isfield(stat, 'posclusters')
for i=1:length(stat.posclusters)
stat.posclusters(i).stddev = sqrt(stat.posclusters(i).prob.*(1-stat.posclusters(i).prob)/Nrand);
stat.posclusters(i).cirange = 1.96*stat.posclusters(i).stddev;
if i==1 && stat.posclusters(i).prob<cfg.alpha && stat.posclusters(i).prob+stat.posclusters(i).cirange>=cfg.alpha
ft_warning('FieldTrip:posCluster_exceeds_alpha', sprintf('The p-value confidence interval of positive cluster #%i includes %.3f - consider increasing the number of permutations!', i, cfg.alpha));
end
end
end
if isfield(stat, 'negclusters')
for i=1:length(stat.negclusters)
stat.negclusters(i).stddev = sqrt(stat.negclusters(i).prob.*(1-stat.negclusters(i).prob)/Nrand);
stat.negclusters(i).cirange = 1.96*stat.negclusters(i).stddev;
if i==1 && stat.negclusters(i).prob<cfg.alpha && stat.negclusters(i).prob+stat.negclusters(i).cirange>=cfg.alpha
ft_warning('FieldTrip:negCluster_exceeds_alpha', sprintf('The p-value confidence interval of negative cluster #%i includes %.3f - consider increasing the number of permutations!', i, cfg.alpha));
end
end
end
if ~isfield(stat, 'prob')
ft_warning('probability was not computed');
else
switch lower(cfg.correctm)
case 'max'
% the correction is implicit in the method
ft_notice('using a maximum-statistic based method for multiple comparison correction\n');
ft_notice('the returned probabilities and the thresholded mask are corrected for multiple comparisons\n');
stat.mask = stat.prob<=cfg.alpha;
stat.posdistribution = posdistribution;
stat.negdistribution = negdistribution;
case 'tfce'
ft_notice('using a threshold free cluster enhancement based method for multiple comparison correction\n');
ft_notice('the returned probabilities and the thresholded mask are corrected for multiple comparisons\n');
stat.mask = stat.prob<=cfg.alpha;
case 'cluster'
% the correction is implicit in the method
ft_notice('using a cluster-based method for multiple comparison correction\n');
ft_notice('the returned probabilities and the thresholded mask are corrected for multiple comparisons\n');
stat.mask = stat.prob<=cfg.alpha;
case 'bonferroni'
ft_notice('performing Bonferroni correction for multiple comparisons\n');
ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
stat.mask = stat.prob<=(cfg.alpha ./ numel(stat.prob));
case 'holm'
% test the most significatt significance probability against alpha/N, the second largest against alpha/(N-1), etc.
ft_notice('performing Holm-Bonferroni correction for multiple comparisons\n');
ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
[pvals, indx] = sort(stat.prob(:)); % this sorts the significance probabilities from smallest to largest
k = find(pvals > (cfg.alpha ./ ((length(pvals):-1:1)')), 1, 'first'); % compare each significance probability against its individual threshold
mask = (1:length(pvals))'<k;
stat.mask = zeros(size(stat.prob));
stat.mask(indx) = mask;
case 'hochberg'
% test the most significant significance probability against alpha/N, the second largest against alpha/(N-1), etc.
ft_notice('performing Hochberg''s correction for multiple comparisons (this is *not* the Benjamini-Hochberg FDR procedure!)\n');
ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
[pvals, indx] = sort(stat.prob(:)); % this sorts the significance probabilities from smallest to largest
k = find(pvals <= (cfg.alpha ./ ((length(pvals):-1:1)')), 1, 'last'); % compare each significance probability against its individual threshold
mask = (1:length(pvals))'<=k;
stat.mask = zeros(size(stat.prob));
stat.mask(indx) = mask;
case 'fdr'
ft_notice('performing FDR correction for multiple comparisons\n');
ft_notice('the returned probabilities are uncorrected, the thresholded mask is corrected\n');
stat.mask = fdr(stat.prob, cfg.alpha);
otherwise
ft_notice('not performing a correction for multiple comparisons\n');
stat.mask = stat.prob<=cfg.alpha;
end
end
% return the observed statistic
if ~isfield(stat, 'stat')
stat.stat = statobs;
end
if exist('statrand', 'var')
stat.ref = mean(statrand,2);
end
% return optional other details that were returned by the statfun
stat = copyfields(statfull, stat, fieldnames(statfull));
ft_warning(ws); % revert to original state