% STD_STAT - compute statistics for ERP/spectral traces or ERSP/ITC images
% of a component or channel cluster in a STUDY.
% Usage:
% >> [pcond, pgroup, pinter, statscond, statsgroup, statsinter] = std_stat( data, 'key', 'val', ...)
% Inputs:
% data - [cell array] mean data for each subject group and/or data
% condition. For example, to compute mean ERPs statistics from a
% STUDY for epochs of 800 frames in two conditions from three
% groups of 12 subjects:
%
% >> data = { [800x12] [800x12] [800x12];... % 3 groups, cond 1
% [800x12] [800x12] [800x12] }; % 3 groups, cond 2
% >> pcond = std_stat(data, 'condstats', 'on');
%
% By default, parametric statistics are computed across subjects
% in the three groups. See below and >> help statcond
% for more information about the statistical computations.
%
% Statistics options (EEGLAB):
% 'groupstats' - ['on'|'off'] Compute (or not) statistics across groups.
% {default: 'off'}
% 'condstats' - ['on'|'off'] Compute (or not) statistics across groups.
% {default: 'off'}
% 'method' - ['parametric'|'permutation'] Type of statistics to use
% default is 'parametric'. 'perm' and 'param' legacy
% abbreviations are still functional.
% 'naccu' - [integer] Number of surrogate averages to accumulate when
% computing permutation-based statistics. For example, to
% test p<0.01 use naccu>=200; for p<0.001, use naccu>=2000.
% If a non-NaN 'threshold' is set (see below) and 'naccu'
% is too low, it will be automatically increased. This
% keyword available only from the command line {default:500}
% 'alpha' - [NaN|p-value] threshold for computing p-value. In this
% function, it is only used to compute naccu above. NaN
% means that no threshold has been set.
% 'mcorrect' - ['none'|'fdr'] apply correcting for multiple comparisons.
% 'mode' - ['eeglab'|'fieldtrip'] statistical framework to use.
% 'eeglab' uses EEGLAB statistical functions and 'fieldtrip'
% uses Fieldtrip statistical functions. Default is 'eeglab'.
%
% Fieldtrip statistics options:
% 'fieldtripnaccu' - 'numrandomization' Fieldtrip parameter
% 'fieldtripalpha' - 'alpha' Fieldtrip parameter. Default is 0.05.
% 'fieldtripmethod' - 'method' Fieldtrip parameter. Default is 'analytic'
% 'fieldtripmcorrect' - ['cluster'|'max'|'fdr'|'holms'|'bonferoni'|'none'].
% Default is 'none'.
% 'fieldtripclusterparam' - string or cell array for optional parameters
% for cluster correction method, see function
% ft_statistics_montecarlo for more information.
% 'fieldtripchannelneighbor' - Fieldtrip channel neighbour structure for
% cluster correction method, see function
% std_prepare_neighbors for more information.
% Legacy parameters:
% 'threshold' - now 'alpha'
% 'statistics' - now 'method'
%
% Outputs:
% pcond - [cell] condition pvalues or mask (0 or 1) if an alpha value
% is selected. One element per group.
% pgroup - [cell] group pvalues or mask (0 or 1). One element per
% condition.
% pinter - [cell] three elements, group pvalues (condition pooled),
% condition pvalues (group pooled) and interaction pvalues.
% statcond - [cell] condition statistic values (F or T).
% statgroup - [cell] group pvalues or mask (0 or 1). One element per
% condition.
% statinter - [cell] three elements, condition statistics (group pooled),
% group statistics (condition pooled) and interaction F statistics.
%
% Author: Arnaud Delorme, CERCO, CNRS, 2006-
%
% See also: STATCOND
% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
function [pcond, pgroup, pinter, statscond, statsgroup, statsinter] = std_stat(data, varargin)
pgroup = {};
pcond = {};
pinter = {};
if nargin < 1
help std_stat;
return;
end
% decode inputs
% -------------
if ~isempty(varargin) && isstruct(varargin{1})
opt = varargin{1};
varargin(1) = [];
else
opt = [];
end
if ~isempty(varargin) ||isempty(opt)
opt = pop_statparams(opt, varargin{:});
end
if ~isfield(opt, 'paired'), opt.paired = { 'off' 'off' }; end
if isfield(opt, 'cluster'), error('Unknown parameter ''cluster'''); end
if ~isnan(opt.eeglab.alpha(1)) && isempty(opt.eeglab.naccu), opt.eeglab.naccu = 1/opt.eeglab.alpha(end)*2; end
if any(any(cellfun('size', data, 2)==1)), opt.groupstats = 'off'; opt.condstats = 'off'; end
if strcmpi(opt.eeglab.mcorrect, 'fdr'), opt.eeglab.naccu = opt.eeglab.naccu*20; end
if isempty(opt.eeglab.naccu), opt.eeglab.naccu = 2000; end
if ~isreal(data{1})
fprintf('*** ITC significance - converting complex values to absolute amplitude ***\n');
for ind = 1:length(data(:))
data{ind} = abs(data{ind});
end
end
nc = size(data,1);
ng = size(data,2);
% compute significance mask
% -------------------------
pcond = {};
pgroup = {};
pinter = {};
statscond = {};
statsgroup = {};
statsinter = {};
if strcmpi(opt.mode, 'eeglab')
% EEGLAB statistics
% -----------------
if strcmpi(opt.condstats, 'on') && nc > 1
for g = 1:ng
[F, df, pval] = statcond(data(:,g), 'method', opt.eeglab.method, 'naccu', opt.eeglab.naccu, 'paired', opt.paired{1});
pcond{g} = squeeze(pval);
statscond{g} = squeeze(F);
end
end
if strcmpi(opt.groupstats, 'on') && ng > 1
for c = 1:nc
[F, df, pval] = statcond(data(c,:), 'method', opt.eeglab.method, 'naccu', opt.eeglab.naccu, 'paired', opt.paired{2});
pgroup{c} = squeeze(pval);
statsgroup{c} = squeeze(F);
end
end
if ( strcmpi(opt.groupstats, 'on') || strcmpi(opt.condstats, 'on') ) && ng > 1 && nc > 1
opt.paired = sort(opt.paired); % put 'off' first if present
[F, df, pval] = statcond(data, 'method', opt.eeglab.method, 'naccu', opt.eeglab.naccu, 'paired', opt.paired{1});
for index = 1:length(pval)
pinter{index} = squeeze(pval{index});
statsinter{index} = squeeze(F{index});
end
end
if ~isempty(opt.groupstats) || ~isempty(opt.condstats)
if ~strcmpi(opt.eeglab.mcorrect, 'none'),
disp([ 'Applying ' upper(opt.eeglab.mcorrect) ' correction for multiple comparisons' ]);
for ind = 1:length(pcond), pcond{ind} = mcorrect( pcond{ind} , opt.eeglab.mcorrect ); end
for ind = 1:length(pgroup), pgroup{ind} = mcorrect( pgroup{ind}, opt.eeglab.mcorrect ); end
if ~isempty(pinter),
pinter{1} = mcorrect(pinter{1}, opt.eeglab.mcorrect);
pinter{2} = mcorrect(pinter{2}, opt.eeglab.mcorrect);
pinter{3} = mcorrect(pinter{3}, opt.eeglab.mcorrect);
end
end
if ~isnan(opt.eeglab.alpha)
for ind = 1:length(pcond), pcond{ind} = applythreshold(pcond{ind}, opt.eeglab.alpha); end
for ind = 1:length(pgroup), pgroup{ind} = applythreshold(pgroup{ind}, opt.eeglab.alpha); end
for ind = 1:length(pinter), pinter{ind} = applythreshold(pinter{ind}, opt.eeglab.alpha); end
end
end
else
if ~exist('ft_freqstatistics'), error('Install Fieldtrip-lite to use Fieldtrip statistics'); end
% Fieldtrip statistics
% --------------------
params = {};
if strcmpi(opt.fieldtrip.mcorrect, 'cluster')
params = eval( [ '{' opt.fieldtrip.clusterparam '}' ]);
if isempty(opt.fieldtrip.channelneighbor), opt.fieldtrip.channelneighbor = struct([]); end
params = { params{:} 'neighbours' opt.fieldtrip.channelneighbor }; % channelneighbor is empty if only one channel selected
end
params = { params{:} 'method', opt.fieldtrip.method, 'naccu', opt.fieldtrip.naccu 'mcorrect' opt.fieldtrip.mcorrect 'alpha' opt.fieldtrip.alpha 'numrandomization' opt.fieldtrip.naccu };
params = { params{:} 'structoutput' 'on' }; % before if ~isnan(opt.fieldtrip.alpha), end
pinter = cell(1,3);
statsinter = cell(1,3);
if strcmpi(opt.condstats, 'on') && nc > 1
newdata = data(:,1);
for g = 1:ng
% marginal effect
[F, df, pval] = statcondfieldtrip(data(:,g), 'paired', opt.paired{1}, params{:});
pcond{g} = applymask( F, opt.fieldtrip);
statscond{g} = squeeze(F.stat);
% concatenate data for main statistics
if g > 1 && ng > 1
for c = 1:nc
switch ndims(data{1,1})
case 2, newdata{c,1}(:,end+1:end+size(data{c,g},2)) = data{c,g};
case 3, newdata{c,1}(:,:,end+1:end+size(data{c,g},3)) = data{c,g};
case 4, newdata{c,1}(:,:,:,end+1:end+size(data{c,g},4)) = data{c,g};
end
end
end
end
% main statistics
if ng > 1
[F, df, pval] = statcondfieldtrip(newdata, 'paired', opt.paired{1}, params{:});
pinter{2} = applymask(F, opt.fieldtrip);
statsinter{2} = squeeze(F.stat);
end
else
pcond = {};
end
if strcmpi(opt.groupstats, 'on') && ng > 1
newdata = data(1,:);
for c = 1:nc
% marginal effect
[F, df, pval] = statcondfieldtrip(data(c,:), 'paired', opt.paired{2}, params{:});
pgroup{c} = applymask( F, opt.fieldtrip);
statsgroup{c} = squeeze(F.stat);
% concatenate data for main statistics
if c > 1 && nc > 1
for g = 1:ng
switch ndims(data{1,1})
case 2, newdata{1,g}(:,end+1:end+size(data{c,g},2)) = data{c,g};
case 3, newdata{1,g}(:,:,end+1:end+size(data{c,g},3)) = data{c,g};
case 4, newdata{1,g}(:,:,:,end+1:end+size(data{c,g},4)) = data{c,g};
end
end
end
end
% main statistics
if nc > 1
[F, df, pval] = statcondfieldtrip(newdata, 'paired', opt.paired{1}, params{:});
pinter{1} = applymask(F, opt.fieldtrip);
statsinter{1} = squeeze(F.stat);
end
else
pgroup = {};
end
end
% apply mask for fieldtrip data
% -----------------------------
function p = applymask(F, fieldtrip)
if ~isnan(fieldtrip.alpha), p = squeeze(F.mask);
else
p = squeeze(F.pval);
if ~strcmpi(fieldtrip.mcorrect, 'none')
p(~F.mask) = 1;
end
end
% apply stat threshold to data for EEGLAB stats
% ---------------------------------------------
function newdata = applythreshold(data, threshold)
threshold = sort(threshold);
newdata = zeros(size(data));
for index = 1:length(threshold)
inds = data < threshold(index);
data(inds) = 1;
newdata(inds) = length(threshold)-index+1;
end
% compute correction for multiple comparisons
% -------------------------------------------
function pvals = mcorrect(pvals, method);
switch method
case {'no' 'none'}, return;
case 'bonferoni', pvals = pvals*prod(size(pvals));
case 'holms', [tmp ind] = sort(pvals(:)); [tmp ind2] = sort(ind); pvals(:) = pvals(:).*(prod(size(pvals))-ind2+1);
case 'fdr', pvals = fdr(pvals);
otherwise error(['Unknown method ''' method ''' for correction for multiple comparisons' ]);
end
Datasets
Models
