forked from LIMO-EEG-Toolbox/limo_tools
-
Notifications
You must be signed in to change notification settings - Fork 1
/
limo_mstat_values.m
354 lines (302 loc) · 12.6 KB
/
limo_mstat_values.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
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
function [M, mask, mytitle] = limo_mstat_values(varargin)
% find corrected p values and mask from data under H0
% this file is more multivariate statistics only
%
% FORMAT [M, mask, mytitle] = limo_stat_values2(varargin)
%
% INPUTS
% Type = the type of plots (matters for title)
% FileName = Name of the file selected
% p = p value for thresholding
% MCC = multiple comparisons option
% 1 none
% 2 2D clustering
% 3 1D clusteriung
% 4 uses max values
% 5 uses TFCE
% LIMO = LIMO.mat
% choice = if MCC = 1 but bootstrap was used then we have the choice between empirical or theoretical p values
% choice = 'use empirical p values' or 'use theoretical p values'
% effect_nb = indicates wich effect or regressor to plot
%
% OUTPUTS
% M the stat values
% mask the significant data
% mytitle is the right title to use to make the figure for type 1 and 2
%
% see limo_display_results limo_stat_values
%
% Cyril Pernet v1 13-06-2013
% --------------------------
% Copyright (C) LIMO Team 2010
Type = varargin{1}; % type of plot
FileName = varargin{2}; % Name of the file selected
p = varargin{3}; % p value
MCC = varargin{4}; % multiple comparison option
LIMO = varargin{5}; % LIMO.mat
choice = [];
effect_nb = [];
if nargin >5
choice = varargin{6}; % if MCC = 1 use empirical or theoretical p values
if nargin == 7
effect_nb = varargin{7}; % for regressions only, indicates wich regressor to plot
end
end
load (FileName);
M = []; mask =[]; mytitle=[];
c = clock; disp(' ');
if MCC ~= 1
fprintf('limo_display_results %gh %gmin %gsec: computing statistical correction...\n',c(4),c(5),c(6));
else
fprintf('limo_display_results %gh %gmin %gsec: making figure...\n',c(4),c(5),c(6));
end
%% Deal with each case of FileName
% ---------
%% R2.mat
% ---------
if strcmp(FileName,'R2.mat')
if strcmp(choice,'Roy')
M = squeeze(R2(:,2)); % F values
else
M = squeeze(R2(:,4));
end
MCC_data = 'H0_R2.mat';
% no correction for multiple testing
% -----------------------------------
if MCC == 1
if LIMO.design.bootstrap == 1
try cd('H0');load(MCC_data); cd ..
if strcmp(choice,'Roy')
H0_F_values = squeeze(H0_R2(:,2,:)); clear H0_R2;
else
H0_F_values = squeeze(H0_R2(:,4,:)); clear H0_R2;
end
sorted_values = sort(H0_F_values,2); clear H0_F_values
U = round((1-p)*size(sorted_values,2));
mask = (M >= sorted_values(:,U));
for column = 1:size(M,2)
tmp(column) = sum(M(column)>squeeze(sorted_values(column,:)));
end
M = 1- (tmp ./ size(sorted_values,2)) ; % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on bootstraped F values');
catch ME
if strcmp(choice,'Roy')
mask = R2(:,3) < p;
M = squeeze(R2(:,3)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Roy test');
else
mask = R2(:,5) < p;
M = squeeze(R2(:,5)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Pillai test');
end
end
else
if strcmp(choice,'Roy')
mask = R2(:,3) < p;
M = squeeze(R2(:,3)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Roy test');
else
mask = R2(:,5) < p;
M = squeeze(R2(:,5)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Pillai test');
end
end
mask = mask' ; M = M';
% cluster correction for multiple testing
% ---------------------------------------
elseif MCC == 2 || MCC == 3
% correction using the max
% --------------------------
elseif MCC == 4 % Stat max
% correction using TFCE
% --------------------------
elseif MCC == 5 % Stat max
end
% ---------------------------------
%% Condition_effect (from 1st level)
% ---------------------------------
elseif strncmp(FileName,'Condition_effect',16)
effect_nb = eval(FileName(18:end-4));
if strcmp(choice,'Roy')
M = squeeze(Condition_effect(:,1)); % F values
else
M = squeeze(Condition_effect(:,3));
end
MCC_data = 'H0_Condition_effect.mat';
% no correction for multiple testing
% -----------------------------------
if MCC == 1
if LIMO.design.bootstrap == 1
try cd('H0');load(MCC_data); cd ..
if strcmp(choice,'Roy')
H0_F_values = squeeze(H0_Condition_effect(:,1,:)); clear H0_Condition_effect;
else
H0_F_values = squeeze(H0_Condition_effect(:,3,:)); clear H0_Condition_effect;
end
sorted_values = sort(H0_F_values,2); clear H0_F_values
U = round((1-p)*size(sorted_values,2));
mask = (M >= sorted_values(:,U));
for column = 1:size(M,2)
tmp(column) = sum(M(column)>squeeze(sorted_values(column,:)));
end
M = 1- (tmp ./ size(sorted_values,2)) ; % p values
mytitle = sprintf('Condition_effect : uncorrected threshold \n based on bootstraped F values');
catch ME
if strcmp(choice,'Roy')
mask = Condition_effect(:,2) < p;
M = squeeze(Condition_effect(:,2)); % p values
mytitle = sprintf('Condition_effect : uncorrected threshold \n based on Roy test');
else
mask = Condition_effect(:,4) < p;
M = squeeze(Condition_effect(:,4)); % p values
mytitle = sprintf('Condition_effect : uncorrected threshold \n based on Pillai test');
end
end
else
if strcmp(choice,'Roy')
mask = Condition_effect(:,2) < p;
M = squeeze(Condition_effect(:,2)); % p values
mytitle = sprintf('Condition_effect : uncorrected threshold \n based on Roy test');
else
mask = Condition_effect(:,4) < p;
M = squeeze(Condition_effect(:,4)); % p values
mytitle = sprintf('Condition_effect : uncorrected threshold \n based on Pillai test');
end
end
mask = mask' ; M = M';
% cluster correction for multiple testing
% ---------------------------------------
elseif MCC == 2 || MCC == 3
% correction using the max
% --------------------------
elseif MCC == 4 % Stat max
% correction using TFCE
% --------------------------
elseif MCC == 5 % Stat max
end
% ---------------------------------
%% Covariate_effect (from 1st level)
% ---------------------------------
elseif strncmp(FileName,'Covariate_effect',16)
if strcmp(choice,'Roy')
M = squeeze(R2(:,2)); % F values
else
M = squeeze(R2(:,4));
end
MCC_data = 'H0_R2.mat';
% no correction for multiple testing
% -----------------------------------
if MCC == 1
if LIMO.design.bootstrap == 1
try cd('H0');load(MCC_data); cd ..
if strcmp(choice,'Roy')
H0_F_values = squeeze(H0_R2(:,2,:)); clear H0_R2;
else
H0_F_values = squeeze(H0_R2(:,4,:)); clear H0_R2;
end
sorted_values = sort(H0_F_values,2); clear H0_F_values
U = round((1-p)*size(sorted_values,2));
mask = (M >= sorted_values(:,U));
for column = 1:size(M,2)
tmp(column) = sum(M(column)>squeeze(sorted_values(column,:)));
end
M = 1- (tmp ./ size(sorted_values,2)) ; % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on bootstraped F values');
catch ME
if strcmp(choice,'Roy')
mask = R2(:,3) < p;
M = squeeze(R2(:,3)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Roy test');
else
mask = R2(:,5) < p;
M = squeeze(R2(:,5)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Pillai test');
end
end
else
if strcmp(choice,'Roy')
mask = R2(:,3) < p;
M = squeeze(R2(:,3)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Roy test');
else
mask = R2(:,5) < p;
M = squeeze(R2(:,5)); % p values
mytitle = sprintf('R^2 : uncorrected threshold \n based on Pillai test');
end
end
mask = mask' ; M = M';
% cluster correction for multiple testing
% ---------------------------------------
elseif MCC == 2 || MCC == 3
% correction using the max
% --------------------------
elseif MCC == 4 % Stat max
% correction using TFCE
% --------------------------
elseif MCC == 5 % Stat max
end
else
errordlg2('unidentified FileName - no thresholding done');
end
end
%% subfunctions which do the actual thresholding
function [mask,cluster_p] = local_clustering(M,P,bootM,bootP,LIMO,MCC,p)
% call field trip functions to do 2D or 1D clustering
%
% M = 3D matrix of observed F values (note for a single electrode the format is 1*time frames*trials)
% P = 3D matrix of observed p values (note for a single electrode the format is 1*time frames*trials)
% bootM = 3D matrix of F values for data bootstrapped under H0
% bootP = 3D matrix of F values for data bootstrapped under H0
% LIMO = LIMO structure - information requested is LIMO.data.chanlocs and LIMO.data.neighbouring_matrix
% MCC = 2 (spatial-temporal clustering) or 3 (temporal clustering)
% p = threshold to apply (note this applied to create clusters and to
% threshold the cluster map)
if MCC == 2
nboot = size(bootM,3);
U = round((1-p)*nboot); % bootstrap threshold
if size(bootM,1)>1 % many electrodes
minnbchan = 2;
expected_chanlocs = LIMO.data.chanlocs;
channeighbstructmat = LIMO.data.neighbouring_matrix;
boot_maxclustersum=zeros(nboot,1); % compute bootstrap clusters
for boot=1:nboot
boot_maxclustersum(boot) = limo_getclustersum(bootM(:,:,boot),bootP(:,:,boot),channeighbstructmat,minnbchan,p);
end
sort_boot_maxclustersum = sort(boot_maxclustersum,2);
maxclustersum_th = sort_boot_maxclustersum(U); % the threshold under H0
[mask,cluster_p] = limo_cluster_test2(M,P,maxclustersum_th,channeighbstructmat,minnbchan,p);
elseif size(bootM,1)==1 % one electrode
th = limo_ecluster_make(squeeze(bootM),squeeze(bootP),p);
sigcluster = limo_ecluster_test(squeeze(M),squeeze(P),th,p);
mask = sigcluster.elec; cluster_p = [];
end
elseif MCC == 3
nboot = size(bootM,3);
U = round((1-p)*nboot); % bootstrap threshold
th = limo_ecluster_make(squeeze(bootM),squeeze(bootP),p);
sigcluster = limo_ecluster_test(squeeze(M),squeeze(P),th,p);
mask = sigcluster.elec;
end
end
function [mask,p_val] = max_correction(M,bootM,p)
% correction for multiple testing using the max stat value
% note this works for bootstrapped data under H0 and for TFCE
%
% M = matrix of observed values (note for a single electrode the format is 1*time frames*trials)
% bootM = matrix of F values for data bootstrapped under H0
% p = threshold to apply
nboot = size(bootM,3);
for boot=1:nboot
data = squeeze(bootM(:,:,boot));
maxM(boot) = max(data(:)); % collect highest value in space and time for each boot
end
U=round((1-p).*nboot);
sortmaxM = sort(maxM); % sort bootstraps
maxF_th = sortmaxM(:,U); % get threshold for each parameter
mask = squeeze(M) >= maxF_th;
for row =1:size(M,1)
for column=1:size(M,2)
p_val(row,column) = 1-(sum(squeeze(M(row,column)) >=sortmaxM) / nboot);
end
end
end