forked from fieldtrip/fieldtrip
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathft_statistics_analytic.m
139 lines (131 loc) · 7.46 KB
/
ft_statistics_analytic.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
function [stat, cfg] = ft_statistics_analytic(cfg, dat, design)
% FT_STATISTICS_ANALYTIC performs a parametric statistical test on the data, based on
% a known (i.e. analytic) distribution of the test statistic. 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 = 'analytic'
%
% The configuration options that can be specified are:
% cfg.statistic = string, statistic to compute for each sample or voxel (see below)
% cfg.correctm = string, apply multiple-comparison correction, 'no', 'bonferroni', 'holm', 'hochberg', 'fdr' (default = 'no')
% cfg.alpha = number, critical value for rejecting the null-hypothesis (default = 0.05)
% cfg.tail = number, -1, 1 or 0 (default = 0)
% 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-block variable(s)
%
% The parametric statistic that is computed for each sample (and for
% which the analytic probability of the null-hypothesis is computed) 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.
%
% See also FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS, FT_SOURCESTATISTICS
% FT_STATISTICS_MONTECARLO, FT_STATISTICS_STATS, FT_STATISTICS_MVPA,
% FT_STATISTICS_CROSSVALIDATE
% Copyright (C) 2006-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, 'renamedval', {'correctm', 'bonferoni', 'bonferroni'});
cfg = ft_checkconfig(cfg, 'renamedval', {'correctm', 'holms', 'holm'});
cfg = ft_checkconfig(cfg, 'renamedval', {'correctm', 'none', 'no'});
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
cfg.correctm = ft_getopt(cfg, 'correctm', 'no');
cfg.alpha = ft_getopt(cfg, 'alpha', 0.05);
cfg.tail = ft_getopt(cfg, 'tail', 0);
% 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
% tell the low-level statfun to compute the statistic, critical values, and the probability
cfg.computestat = 'yes';
cfg.computecritval = 'yes';
cfg.computeprob = 'yes';
% perform the statistical test
if nargout(statfun)>1
[stat, cfg] = statfun(cfg, dat, design);
else
[stat] = statfun(cfg, dat, design);
end
% these are only intenced for the low-level statfun, and should not be returned
cfg = rmfield(cfg, 'computestat');
cfg = rmfield(cfg, 'computeprob');
cfg = rmfield(cfg, 'computecritval');
if ~isfield(stat, 'prob')
ft_warning('probability was not computed');
else
switch lower(cfg.correctm)
case 'bonferroni'
ft_notice('performing Bonferoni 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);
case 'no'
ft_notice('not performing a correction for multiple comparisons\n');
stat.mask = stat.prob<=cfg.alpha;
otherwise
ft_error('unsupported option "%s" for cfg.correctm', cfg.correctm);
end
end