-
Notifications
You must be signed in to change notification settings - Fork 9
/
ft_sourcestatistics.m
245 lines (214 loc) · 9.17 KB
/
ft_sourcestatistics.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
function [stat] = ft_sourcestatistics(cfg, varargin)
% FT_SOURCESTATISTICS computes the probability for a given null-hypothesis using
% a parametric statistical test or using a non-parametric randomization test.
%
% Use as
% [stat] = ft_sourcestatistics(cfg, source1, source2, ...)
% where the input data is the result from FT_SOURCEANALYSIS, FT_SOURCEDESCRIPTIVES
% or FT_SOURCEGRANDAVERAGE. The source structures should be spatially alligned
% to each other and should have the same positions for the sourcemodel.
%
% The configuration should contain the following option for data selection
% cfg.parameter = string, describing the functional data to be processed, e.g. 'pow', 'nai' or 'coh'
%
% Furthermore, the configuration should contain:
% cfg.method = different methods for calculating the probability of the null-hypothesis,
% 'montecarlo' uses a non-parametric randomization test to get a Monte-Carlo estimate of the probability,
% 'analytic' uses a parametric test that results in analytic probability,
% 'stats' (soon deprecated) uses a parametric test from the MATLAB statistics toolbox,
%
% The other cfg options depend on the method that you select. You
% should read the help of the respective subfunction FT_STATISTICS_XXX
% for the corresponding configuration options and for a detailed
% explanation of each method.
%
% See also FT_SOURCEANALYSIS, FT_SOURCEDESCRIPTIVES, FT_SOURCEGRANDAVERAGE, FT_MATH,
% FT_STATISTICS_MONTECARLO, FT_STATISTICS_ANALYTIC, FT_STATISTICS_CROSSVALIDATE, FT_STATISTICS_STATS
% FIXME the following needs to be reimplemented
%
% You can restrict the statistical analysis to regions of interest (ROIs)
% or to the average value inside ROIs using the following options:
% cfg.atlas = filename of the atlas
% cfg.roi = string or cell of strings, region(s) of interest from anatomical atlas
% cfg.avgoverroi = 'yes' or 'no' (default = 'no')
% cfg.hemisphere = 'left', 'right', 'both', 'combined', specifying this is required when averaging over regions
% Copyright (C) 2005-2020, 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$
% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin = nargin;
ft_nargout = nargout;
% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar varargin
ft_preamble provenance varargin
ft_preamble trackconfig
ft_preamble randomseed
% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
return
end
% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'required', {'method', 'design'});
cfg = ft_checkconfig(cfg, 'renamed', {'approach', 'method'});
cfg = ft_checkconfig(cfg, 'forbidden', {'transform'});
cfg = ft_checkconfig(cfg, 'forbidden', {'trials'}); % this used to be present until 24 Dec 2014, but was deemed too confusing by Robert
cfg = ft_checkconfig(cfg, 'forbidden', {'channel'});
% check if the input data is valid for this function
for i=1:length(varargin)
varargin{i} = ft_checkdata(varargin{i}, 'datatype', 'source', 'feedback', 'no');
varargin{i} = ft_datatype_source(varargin{i}, 'version', 'upcoming');
end
% set the defaults
cfg.parameter = ft_getopt(cfg, 'parameter'); % default is set below
cfg.correctm = ft_getopt(cfg, 'correctm');
cfg.latency = ft_getopt(cfg, 'latency', 'all');
cfg.avgovertime = ft_getopt(cfg, 'avgovertime', 'no');
cfg.frequency = ft_getopt(cfg, 'frequency', 'all');
cfg.avgoverfreq = ft_getopt(cfg, 'avgoverfreq', 'no');
cfg.parameter = ft_getopt(cfg, 'parameter', 'pow');
if strncmp(cfg.parameter, 'avg.', 4)
cfg.parameter = cfg.parameter(5:end); % remove the 'avg.' part
end
for i=1:length(varargin)
assert(isfield(varargin{i}, cfg.parameter), 'data does not contain parameter "%s"', cfg.parameter);
end
% ensure that the data in all inputs has the same channels, time-axis, etc.
tmpcfg = keepfields(cfg, {'frequency', 'avgoverfreq', 'latency', 'avgovertime', 'avgoverpos', 'parameter', 'showcallinfo', 'select', 'nanmean'});
[varargin{:}] = ft_selectdata(tmpcfg, varargin{:});
% restore the provenance information
[cfg, varargin{:}] = rollback_provenance(cfg, varargin{:});
dimord = getdimord(varargin{1}, cfg.parameter);
dimtok = tokenize(dimord, '_');
dimsiz = getdimsiz(varargin{1}, cfg.parameter, numel(dimtok));
rptdim = find( strcmp(dimtok, 'subj') | strcmp(dimtok, 'rpt') | strcmp(dimtok, 'rpttap'));
datdim = find(~strcmp(dimtok, 'subj') & ~strcmp(dimtok, 'rpt') & ~strcmp(dimtok, 'rpttap'));
datsiz = dimsiz(datdim);
% pass these fields to the low-level functions, they should be removed further down
cfg.dimord = sprintf('%s_', dimtok{datdim});
cfg.dimord = cfg.dimord(1:end-1); % remove trailing _
cfg.dim = dimsiz(datdim);
if isfield(varargin{1}, 'dim')
% the positions can be mapped onto a 3-D volume
if prod(cfg.dim)==prod(varargin{1}.dim)
% one value per grid position, i.e. [nx ny nz]
cfg.origdim = varargin{1}.dim;
else
% multiple value per grid position
cfg.origdim = [varargin{1}.dim cfg.dim(2:end)];
end
end
if isfield(varargin{1}, 'inside')
cfg.inside = varargin{1}.inside;
else
cfg.inside = true(size(varargin{1}.pos,1),1);
end
% also create an inside vector for the reshaped data, which is needed
% when there are multiple values per grid position
cfg.originside = cfg.inside;
cfg.inside = repmat(cfg.inside, prod(cfg.dim)./numel(cfg.inside), 1);
if numel(cfg.dim)==1
cfg.dim(2) = 1; % add a trailing singleton dimensions
end
if isempty(rptdim)
% repetitions are across multiple inputs
dat = nan(prod(dimsiz), length(varargin));
for i=1:length(varargin)
tmp = varargin{i}.(cfg.parameter);
dat(:,i) = tmp(:);
end
else
% repetitions are within inputs
dat = cell(size(varargin));
for i=1:length(varargin)
tmp = varargin{i}.(cfg.parameter);
if rptdim~=1
% move the repetitions to the first dimension
tmp = permute(tmp, [rptdim datdim]);
end
dat{i} = reshape(tmp, size(tmp,1), []);
end
dat = cat(1, dat{:}); % repetitions along 1st dimension
dat = dat'; % repetitions along 2nd dimension
end
dat = dat(cfg.inside,:);
if size(cfg.design,2)~=size(dat,2)
cfg.design = transpose(cfg.design);
end
design = cfg.design;
% fetch function handle to the intermediate-level statistics function
statmethod = ft_getuserfun(cfg.method, 'statistics');
if isempty(statmethod)
ft_error('could not find the corresponding function for cfg.method="%s"\n', cfg.method);
else
ft_info('using "%s" for the statistical testing\n', func2str(statmethod));
end
% check that the design completely describes the data
if size(dat,2) ~= size(cfg.design,2)
ft_error('the length of the design matrix (%d) does not match the number of observations in the data (%d)', size(cfg.design,2), size(dat,2));
end
% determine the number of output arguments
try
% the nargout function in MATLAB 6.5 and older does not work on function handles
num = nargout(statmethod);
catch
num = 1;
end
% perform the statistical test
if num>1
[stat, cfg] = statmethod(cfg, dat, design);
cfg = rollback_provenance(cfg); % ensure that changes to the cfg are passed back to the right level
else
[stat] = statmethod(cfg, dat, design);
end
if ~isstruct(stat)
% only the probability was returned as a single matrix, reformat into a structure
stat = struct('prob', stat);
end
% the statistical output contains multiple elements, e.g. F-value, beta-weights and probability
fn = fieldnames(stat);
for i=1:length(fn)
if numel(stat.(fn{i}))==sum(cfg.inside)
% get the data back onto the inside positions
tmp = nan+zeros(numel(cfg.inside),1);
tmp(cfg.inside) = stat.(fn{i});
stat.(fn{i}) = tmp;
end
if numel(stat.(fn{i}))==prod(datsiz)
% reformat into the same dimensions as the input data
stat.(fn{i}) = reshape(stat.(fn{i}), [datsiz 1]);
end
end
% describe the dimensions of the output data
stat.dimord = cfg.dimord;
% copy the descripive fields into the output
stat = copyfields(varargin{1}, stat, {'freq', 'time', 'pos', 'dim', 'transform', 'tri', 'inside'});
% these were only present to inform the low-level functions
cfg = removefields(cfg, {'dimord', 'tri', 'dim', 'origdim', 'inside', 'originside'});
% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble randomseed
ft_postamble trackconfig
ft_postamble previous varargin
ft_postamble provenance stat
ft_postamble history stat
ft_postamble savevar stat