forked from SnPM-toolbox/SnPM-devel
-
Notifications
You must be signed in to change notification settings - Fork 0
/
snpm_pp.m
1758 lines (1561 loc) · 61.9 KB
/
snpm_pp.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
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function snpm_pp(CWD,varargin)
% SnPM post processing and results display
% FORMAT snpm_pp(CWD)
%
% CWD - Directory containing SnPM results files
%
% If CWD is not specified then user is prompted to locate results file SnPM.mat
%_______________________________________________________________________
%
% snpm_pp is the PostProcessing function for the SnPM nonParametric
% statistical analysis. SnPM statistical analyses are split into three
% stages; Setup, Compute & Assess. This is the third stage.
% Nonparametric randomisation distributions are read in from MatLab
% *.mat files, with which the observed statistic image is assessed
% according to user defined parameters. It is the SnPM equivalent of
% the "Results" section of SPM, albeit with reduced features.
%
% Voxel level corrected p-values are computed from the permutation
% distribution of the maximal statistic. If suprathreshold cluster
% statistics were collected in the computation stage (and the large
% SnPM_ST.mat file hasn't been deleted!), then assessment by
% suprathreshold cluster size is also available, using a user-specified
% primary threshold.
%
% Instructions:
%=======================================================================
%
% You are prompted for the following:
%
% (1) ResultsDir: If the results directory wasn't specified on the command
% line, you are prompted to locate the SnPM results file SnPM.mat.
% The directory in which this file resides is taken to be the
% results directory, which must contain *all* the files listed
% below ("SnPM files required").
%
% Results (spm.ps & any requested image files) are written in the
% present working directory, *not* the directory containing the
% results of the SnPM computations.
%
% ----------------
%
% (2) +/-: Having located and loaded the results files, you are asked to
% chose between "Positive or negative effects?". SnPM, like SPM,
% only implements single tailed tests. Choose "+ve" if you wish to
% assess the statistic image for large values, indicating evidence
% against the null hypothesis in favour of a positive alternative
% (activation, or positive slope in a covariate analysis).
%
% Choose "-ve" to assess the negative contrast, i.e. to look for
% evidence against the null hypothesis in favour of a negative
% alternative (de-activation, or a negative slope in a covariate
% analysis). The "-ve" option negates the statistic image and
% contrast, acting as if the negative of the actual contrast was
% entered.
%
% A two-sided test may be constructed by doing two separate
% analyses, one for each tail, at half the chosen significance
% level, doubling the resulting p-values.
% ( Strictly speaking, this is not equivalent to a rigorous two-sided )
% ( non-parametric test using the permutation distribution of the )
% ( absolute maximum statistic, but it'll do! )
%
% ----------------
%
% (3) WriteFiles: All image files have been written by snpm_cp, and there
% is only one file that may be written here: Users are asked 'write
% filtered statistic image?'.
% ----------------
%
% Next come parameters for the assessment of the statistic image...
%
% (4) alpha: (p-value for filtering)
% First, you are asked 'Use corrected threshold?' You can choose to set
% the threshold as FWE-corrected('FWE'), FDR-corrected('FDR') or
% uncorrected('None'). A FWE threshold controls the chance of one or
% more false positives; a FDR threshold controls the expected
% fraction of false positives among the voxels detected.
%
% Then you specify the threshold, the statistical significance
% level at which you wish to assess the evidence against the null
% hypothesis. In SPM this is called "filtering by corrected
% p-value". SnPM will only show you voxels (& suprathreshold
% regions if you choose) that are significant based on the method you
% choose at level \alpha. I.e. only voxels (& regions) with
% corrected (or uncorrected) p-value less than \alpha are shown to
% you.
%
% Setting \alpha to 1 will show you all voxels with a positive statistic.
%
%
% (5) SpatEx: If you collected supra-threshold cluster statistics during
% the SnPM computation phase, you are offered the option to assess
% the statistic image by supra-threshold cluster size (spatial
% extent).
%
% 5a) ST_Ut: If you chose to asses spatial extent, you are now prompted
% for the primary threshold. This is the threshold applied to the
% statistic image for the identification of supra-threshold
% clusters.
%
% The acceptable range is limited. SnPM has to collect
% suprathreshold information for every relabelling. Rather that
% pre-specify the primary threshold, information is recorded for
% each voxel exceeding a low threshold (set in snpm_cp) for every
% permutation. From this, suprathreshold cluster statistics can be
% generated for any threshold higher than the low recording
% threshold. This presents a lower limit on the possible primary
% threshold.
%
% The upper limit (if specified) corresponds to the statistic value
% at which voxels become individually significant at the chosen
% level (\alpha). There is little point perusing a suprathreshold
% cluster analysis at a threshold at which the voxels are
% individually significant.
%
% If the statistics are t-statistics, then you can also specify the
% threshold via the upper tail probability of the t-distribution.
%
% (NB: For the moment, \alpha=1 precludes suprathreshold analysis, )
% ( since all voxels are significant at \alpha=1. )
%
%
% That's it. SnPM will now compute the appropriate significances,
% reporting its progress in the MatLab command window. Note that
% computing suprathreshold cluster size probabilities can take a long
% time, particularly for low thresholds or large numbers of
% relabellings. Eventually, the Graphics window will come up and the
% results displayed.
%
% - Results
%=======================================================================
%
% The format of the results page is similar to that of SPM:
%
% A Maximum Intensity Projection (MIP) of the statistic image is shown
% top left: Only voxels significant (corrected) at the chosen level
% \alpha are shown. (If suprathreshold cluster size is being assessed,
% then clusters are shown if they have significant size *or* if they
% contain voxels themselves significant at the voxel level.) The MIP is
% labelled SnPM{t} or SnPM{Pseudo-t}, the latter indicating that
% variance smoothing was carried out.
%
% On the top right a graphical representation of the Design matrix is
% shown, with the contrast illustrated above.
%
% The lower half of the output contains the table of p-values and
% statistics, and the footnote of analysis parameters. As with SPM, the
% MIP is tabulated by clusters of voxels, showing the maximum voxel
% within each cluster, along with at most three other local maxima
% within the cluster. The table has the following columns:
%
%
% At Cluster Level
%
% * Pcorrected: If "spatial extent" has been assessed, FWE-corrected
% p-values for region size will be shown. The corrected P-value for
% the suprathreshold cluster size is the probability (conditional on
% the data) of the experiment giving a suprathreshold cluster of size
% as or more extreme anywhere in the statistic image. This is the
% proportion of the permutation distribution of the maximal
% suprathreshold cluster size exceeding (or equalling) the observed
% size of the current cluster.
%
% * k: The size (in voxels) of the cluster.
%
%
% At Voxel Level
%
% * PFWE-corr: FWE-corrected non-parametric P-values. This is the
% probability of the experiment giving a voxel statistic this extreme
% anywhere in the statistic image. This is the proportion of the
% permutation distribution of the maximal statistic exceeding (or
% equalling) the observed statistic value.
%
% * PFDR-corr: FDR-corrected non-parametric P-values. This is the
% smallest FDR level for which the voxel would be significant. A level
% q FDR threshold ensures that, on average, no more than q*100% of the
% suprathreshold voxels will be false positives.
%
% * t / Pseudo-t: The statistic value.
%
% * Puncorrected: uncorrected non-parametric P-values. This is the
% result of a non-parametric permutation test completed at this
% individual voxel. This P-value is is the proportion of the
% permutation distribution (at *this* voxel) exceeding (or
% equalling) the observed statistic value.
%
%
% * {x,y,z} mm: Locations of local maxima.
%
% The SnPM parameters footnote contains the following information:
%
% * Primary threshold: If assessing "spatial extent", the primary
% threshold used for identification of suprathreshold clusters is
% printed. If using t-statistics (as opposed to Pseudo-t's), the
% corresponding upper tail probability is also given.
%
% * Critical STCS: The critical suprathreshold cluster size. This is
% size above which suprathreshold clusters have significant size at
% level \alpha It is computed as the 100(1-alpha)%-ile of the
% permutation distribution of the maximal suprathreshold cluster
% size. Only shown when assessing "spatial extent".
%
% * alpha: The test level specified.
%
% * Critical threshold: The critical statistic level. This is the value
% above which voxels are significant (corrected) at level \alpha. It
% is computed as the 100(1-alpha)%-ile of the permutation
% distribution of the maximal statistic.
%
% * df: The degrees of freedom of the t-statistic. This is printed even
% if
% variance smoothing is used, as a guide.
%
% * Volume & voxel dimensions:
%
% * Design: Description of the design
%
% * Perms: Description of the exchangability and permutations used.
%
%
%
% - SnPM files required:
%=======================================================================
% snpm_pp loads parameters and results from the following files, which
% must all be in the same directory:
% SnPMcfg.mat - SnPM design configuration
% SnPM.mat - SnPM analysis & permutation distribution
% SnPMt.mat - Pointlist of (Pseudo) t-statisic for actual labelling
% XYZ.mat - Co-ordinates of pointlist
% SnPM_ST.mat (*) - Suprathreshold cluster statistics (if required)
%
% Further details of the actual variables required from these files are
% given in the main body of snpm_pp
%
% (*) The SnPM_ST.mat file containing the suprathreshold cluster
% information for each of the relabellings can be very large, and is
% only needed if a suprathreshold cluster size test is required. If
% such an analysis is not required, but suprathreshold cluster stats
% were collected, then this file may be deleted, without compromising
% further voxel-level analyses.
%
%_______________________________________________________________________
% Copyright (C) 2013 The University of Warwick
% Id: snpm_pp.m SnPM13 2013/10/12
% Andrew Holmes, Thomas Nichols, Jun Ding, Darren Gitelman
%-----------------------------functions-called------------------------
% spm_DesMtx
% spm_Tcdf
% spm_clf
% spm_clusters
% spm_figure
% spm_select
% spm_hwrite
% spm_input
% spm_invTcdf
% spm_max
% snpm_mip
% snpm_pp
% spm_str_manip
% spm_t2z
% spm_type
% spm_xyz2e
%-----------------------------functions-called------------------------
%-Variable "decoder" - Following files/variables are required:
%=======================================================================
% NB: Mat files contain additional variables beyond those required here.
% See function that wrote each file for full definitions.
% SnPM design configuration file: SnPMcfg.mat
%-----------------------------------------------------------------------
% - saved by snpm_ui
% H condition partition of DesMtx for correctly labeled data
% C covariate partition of DesMtx for correctly labeled data
% B block partition of DesMtx for correctly labeled data
% G confound partition of DesMtx for correctly labeled data
% HCBGnames string matrix of column names of [H C B G]
% PiCond Permuted conditions matrix, one labelling per row, actual
% labelling on first row
% sPiCond String describing permutations in PiCond
% bhPerms Flag indicating use of "half permutations" trick
% CONT Contrast (only one)
% bVarSm Flag for variance smoothing (Pseudo t-statistics)
% sVarSm Sring describing variance Smoothing (empty if bVarSm=0)
% bST Flag for collection of superthreshold info
% sDesign Description of PlugIn design
%
% SnPM analysis & permutation distribution file: SnPM.mat
%-----------------------------------------------------------------------
% - saved by snpm_cp
% S - Volume analyzed (in voxels)
% V Image file handles (see spm_vol)
% df Residual degrees of freedom of raw t-statistic
% MaxT 2xnPerm matrix of [max;min] t-statistics per perm
% ST_Ut Threshold above which suprathreshold info was collected.
% Voxel locations, t and perm are saved in SnPM_ST.mat for
% t's greater than ST_Ut. ST_Ut=Inf if not saving STCdata
%
% Pointlist of (Pseudo) t-statisic for actual labelling:SnPMt.mat
%-----------------------------------------------------------------------
% - saved by snpm_cp
% SnPMt - 1xS matrix of voxel statistics (t or pseudo-t)
%
% Co-ordinates of pointlist: XYZ.mat
%-----------------------------------------------------------------------
% - saved by snpm_cp
% XYZ - 3xS matrix of co-ordinates [x;y;z] of voxels on SnPMt
%
% Suprathreshold cluster statistics: SnPM_ST.mat
%-----------------------------------------------------------------------
% - saved by snpm_cp
% SnPM_ST - Suprathreshold cluster statistics, see snpm_cp.m
% NB: This file is only required for suprathreshold cluster size analysis
%-Setup
%=======================================================================
global defaults
if isempty(defaults), spm_defaults; end
global SnPMdefs
if isempty(SnPMdefs), snpm_defaults; end
MLver=version;MLver=MLver(1);
fprintf('\nSnPM: snpm_pp\n'),fprintf('%c','='*ones(1,72)),fprintf('\n')
%-Initialise variables & constants
%-----------------------------------------------------------------------
tol = 1e-4; % Tolerance for comparing real numbers
% Two reals with abs(a-b)<tol are considered equal
% ( Reals have to be compared for equality when )
% ( computing FWE-corrected p-values )
Dis = SnPMdefs.Results_distmin; % mm distance between sub regions
Num = SnPMdefs.Results_nbmax; % Maximum number of sub-regions
units = {'mm' 'mm' 'mm'};
%-SetUp figure window
%-----------------------------------------------------------------------
Finter = spm_figure('FindWin','Interactive');
Fgraph = spm_figure('FindWin','Graphics');
if isempty(Fgraph), Fgraph=spm_figure('Create','Graphics'); end
spm_clf(Finter), spm_clf(Fgraph)
set(Finter,'Name','SnPM PostProcess');
%-Get Data
%=======================================================================
% Get analysis directory
if nargin==0
tmp = spm_select(1,'SnPM.mat','Select SnPM.mat for analysis...');
CWD = spm_str_manip(tmp,'hd');
end
if nargin>1
job=varargin{1};
BATCH=true;
else
BATCH=false;
end
%-Skip reports in BATCH; will create them later
if BATCH
Report = {job.Report};
else
% Only can specify multiple reports outside of BATCH mode; otherwise the results will fly by,
% one after another without pausing
Report = {'FWEreport','FDRreport','MIPtable'};
end
%-Load Config file & SnPM permutation data
load(fullfile(CWD,'SnPMcfg'))
load(fullfile(CWD,'SnPM'))
load(fullfile(CWD,'SnPMucp'))
%-Ask whether positive or negative effects be analysed
%-----------------------------------------------------------------------
if BATCH
if STAT == 'T'
bNeg = job.Tsign==-1;
else
bNeg = 0;
end
else
if STAT == 'T'
bNeg = spm_input('Positive or negative effects?',1,'b','+ve|-ve',[0,1],1);
else
bNeg = 0;
str = 'Positive effects';
spm_input('F-statistic, so only:','+1','b',str,1);
end
end
%-Take MaxT for increases or decreases according to bNeg
MaxT = MaxT(:,bNeg+1);
nPerm = size(PiCond,1); %nPerm is consistent with the one in snpm_cp
nPermReal = size(MaxT,1); %different with nPerm when bhPerms==1
[StMaxT, iStMaxT] = sort(MaxT);
%-Load statistic image
%-----------------------------------------------------------------------
load(fullfile(CWD,'SnPMt'))
load(fullfile(CWD,'XYZ'))
XYZ0=XYZ;
%-Negate if looking at negative contrast
%-----------------------------------------------------------------------
if bNeg
SnPMt = -SnPMt;
CONT = -CONT;
end
%-Get ORIGIN, etc
DIM = [V(1).dim(1) V(1).dim(2) V(1).dim(3)];
M=V(1).mat(1:3, 1:3);
VOX=sqrt(diag(M'*M))';
MAT = V(1).mat;
IMAT = inv(MAT);
ORIGIN = IMAT(1:3,4);
% Template vol structure
Vs0 = V(1);
% Vs0 = struct('fname', '',...
% 'dim', [DIM,spm_type('float')],...
% 'mat', MAT,...
% 'pinfo', [1 0 0]',...
% 'descrip', '');
% Process Nonparmaetric P-values
if bNeg
% Here, nPermReal has already been doubled if bhPerms=1
SnPMucp = 1+1/nPermReal-SnPMucp;
end
sSnPMucp = sort(SnPMucp);
%-Write out filtered statistic image? (Get's done later)
%-----------------------------------------------------------------------
if BATCH
if isfield(job.WriteFiltImg,'WF_no')
WrtFlt=0;
else
WrtFlt=1;
WrtFltFn=job.WriteFiltImg.name;
if isempty(spm_str_manip(WrtFltFn, 'e'))
WrtFltFn = [WrtFltFn '.nii'];
end
end
else
WrtFlt = spm_input('Write filtered statistic img?','+1','y/n',[1,0],2);
if WrtFlt
WrtFltFn = 'SnPMt_filtered';
WrtFltFn=spm_input('Filename ?','+1','s',WrtFltFn);
WrtFltFn = [WrtFltFn, '.img'];
end
end
%-Get inference parameters
%=======================================================================
% Map of options from Batch System
% [Root] > .Thr
% Voxel-Level > .Vox
% . Significance > .VoxSig
% . . Uncorrected Nonparametric P | Uncorrected T or F | FDR Corrected | FWE Corrected
% > .Pth .TFth .FDRth .FWEth
% Cluster-Level > .Clus
% . Cluster size statistic > .ClusSize
% . . Cluster-Forming Threshold > .CFth
% . . Significance Level > .ClusSig
% . . . Uncorrected k | FWE Corrected
% > .Cth .FWEthC
%-Get corrected threshold
%-----------------------------------------------------------------------
u = NaN; % Statistic Image threshold
alpha_ucp = NaN; % Uncorrected P-value image threshold
C_MaxT = NaN; % Statistic image threshold set by alph_FWE
C_STCS = NaN; % Cluster size threshold (set directly by uncorrected
% threshold or by alph_FWE)
alph_FWE = NaN; % FWE rate of a specified u threshold
alph_FDR = NaN; % FDR rate of a specified alpha_ucp
if BATCH
bSpatEx = isfield(job.Thr,'Clus');
if ~bSpatEx
% Voxel-wise inference
tmp = fieldnames(job.Thr.Vox.VoxSig);
switch tmp{1}
case 'Pth'
alpha_ucp = BoundCheck(job.Thr.Vox.VoxSig.Pth,[0 1],'Invalid Uncorrected P');
alph_FDR = snpm_P_FDR(alpha_ucp,[],'P',[],sSnPMucp');
case 'TFth'
u = BoundCheck(job.Thr.Vox.VoxSig.TFth,[0 Inf],'Negative Threshold!');
alph_FWE = sum(MaxT > u -tol) / nPermReal;
case 'FDRth'
alph_FDR = BoundCheck(job.Thr.Vox.VoxSig.FDRth,[0 1],'Invalid FDR level');
alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp');
case 'FWEth'
alph_FWE = BoundCheck(job.Thr.Vox.VoxSig.FWEth,[0 1],'Invalid FWE level');
iFWE = ceil((1-alph_FWE)*nPermReal);
if alph_FWE<1
C_MaxT=StMaxT(iFWE);
else
C_MaxT = 0;
end
u = C_MaxT;
end
else
% Cluster-wise inference
if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 && exist(fullfile(CWD,'STCS.mat'))~=2
error(['SnPM:NoClusterInfo', 'ERROR: Cluster-wise inference requested, but no cluster information saved.\n',...
'Re-configure analysis changing "Cluster inference" to "Yes" and re-run.\n'])
end
%%% Sort out the cluster-forming threshold
if pU_ST_Ut==-1 % No threshold was set in snpm_ui.
if isnan(job.Thr.Clus.ClusSize.CFth)
error('SnPM:NoClusterFormingThresh', 'ERROR: Cluster-forming threshold set to NaN in results with "slow" cluster inference method used in compoutation. \nRe-run results specifying a cluster-forming threshold.\n')
end
% Save original ST_Ut
ST_Ut_0 = ST_Ut;
CFth=job.Thr.Clus.ClusSize.CFth;
if (CFth<=0)
error('SnPM:InvalidClusterFormingThresh', 'ERROR: Cluster-forming threshold must be strictly positive.\nRe-run results with a cluster-forming threshold greater than 0.\n')
end
if bVarSm
%-If using pseudo-statistics then can't use (uncorrected)
% upper tail p-values to specify primary threshold
pCFth = NaN;
if (CFth<1)
warning('snpm_pp:pseudoTFormingThresholdP',...
['Pseudo-T cluster-forming threshold defined by '...
'P-value using Gaussian approximation P=' num2str(pU_ST_Ut)...
' -> Z=' num2str(ST_Ut) '; actual Pseudo-T threshold '...
'unknown but may be higher than ' num2str(ST_Ut) '.']);
pCFth = CFth;
CFth = spm_invNcdf(1-CFth);
% error(sprintf('ERROR: Cluster-forming threshold specified as a P-value (%g), but uncorrected P-values are unavailable for the pseudo t (smoothed variance t-test). \nRe-run results with a cluster-forming threshold greater than 1.\n',ST_Ut))
end
if (CFth < ST_Ut)%(CFth>=ST_Ut-tol)
if isnan(pCFth)
error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of %0.2f specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming threshold of %0.2f or higher. (Alternatively, increase SnPMdefs.STprop in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',CFth,ST_Ut,ST_Ut))
else
error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of P=%0.4f (T=%0.2f) specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming P-value threshold of %0.2f or lower. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',pCFth,CFth,ST_Ut,p_ST_Ut))
end
end
else
%-Statistic image is t with df degrees of freedom
p_ST_Ut = STalpha;
if (CFth < 1)
pCFth = CFth;
CFth = spm_invTcdf(1-CFth,df);
else
pCFth = NaN;
if (abs(CFth-ST_Ut)<=tol)
CFth=ST_Ut; % If tmp is very close to ST_Ut, set tmp equal to ST_Ut.
end
end
if (CFth < ST_Ut) %(CFth>=ST_Ut-tol)
if isnan(pCFth) % statistic-value cluster-forming threshold
error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of %0.2f specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming threshold of %0.2f or higher. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',CFth,ST_Ut,ST_Ut))
else
error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of P=%0.4f (T=%0.2f) specified, but statistic image information only saved for %0.2f and greater. \nRe-run results with a cluster-forming P-value threshold of %0.2f or lower. (Alternatively, increase SnPMdefs.STalpha in snpm_defaults.m, re-start SnPM, and re-compute analysis.)\n',pCFth,CFth,ST_Ut,p_ST_Ut))
end
end
end
if (abs(CFth-ST_Ut)<=tol)
CFth = ST_Ut; % If tmp is very close to ST_Ut, set tmp equal to ST_Ut.
end
ST_Ut = CFth;
else % Threshold *was* set in snpm_ui.
if ~isnan(job.Thr.Clus.ClusSize.CFth)
error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of T=%0.2f was already set during analysis configuration; hence, in results, cluster-forming threshold must be left as "NaN".\nRe-run results with cluster-forming threshold set to NaN.\n',ST_Ut))
end
end
u=ST_Ut; % Flag use of a statistic-value threshold
% Inference details...
tmp = fieldnames(job.Thr.Clus.ClusSize.ClusSig);
switch tmp{1}
case 'Cth'
C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth;
case 'PthC'
alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)');
case 'FWEthC'
alph_FWE = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.FWEthC,[0 1],'Invalid FWE level (cluster-level inference)');
iFWE = ceil((1-alph_FWE)*nPermReal);
end
end % END: Cluster-wise inference
else % GUI, interative inference specification
str_img =[STAT,'|P'];
switch spm_input('Results for which img?','+1','b',str_img,[],1)
case 'P' % Use the P-image
bSpatEx = 0; % Cluster-wise inference won't be performed anyway.
str = 'FDR|None';
switch spm_input('Use corrected threshold?','+1','b',str,[],1)
case 'FDR' % False discovery rate
%---------------------------------------------------------------
alph_FDR = spm_input('FDR-Corrected p value threshold','+0','r',0.05,1,[0,1]);
alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp');
otherwise %-Uncorrected: no adjustment
%%%% Now ask: Threshold statistic image or Uncorr P-value Image ?
%%%% If stats image, maybe do ST, no FDR-level of thresh;
%%%% If P-value image, no ST, can find FDR-level of thresh
% p for conjunctions is p of the conjunction SPM
%---------------------------------------------------------------
alpha_ucp = spm_input('Uncorrected p value threshold','+0','r',0.01,1,[0,1]);
alph_FDR = snpm_P_FDR(alpha_ucp,[],'P',[],sSnPMucp');
end
case STAT % Use the T-image
%-Ask whether SupraThreshold cluster size test required
%-----------------------------------------------------------------------
%-To have cluster size inference, need
% 1. Spatial extent information was collected (bST=1),
% 2. SnPM_ST.mat or STCS.mat file exists
bSpatEx = bST & (exist(fullfile(CWD,'SnPM_ST.mat'))==2|exist(fullfile(CWD,'STCS.mat'))==2);
if bSpatEx
str = 'Voxelwise|Clusterwise';
bSpatEx = spm_input('Inference method?','+1','b',str,[1 2],1)==2;
else
str = 'Voxelwise';
spm_input('Inference method:','+1','b',str,1);
end
if ~bSpatEx % Voxel-wise inference
str = 'FWE|None';
switch spm_input('Voxelwise: Use corrected thresh?','+1','b',str,[],1)
case 'FWE' % family-wise false positive rate
%---------------------------------------------------------------
alph_FWE = spm_input('FWE-Corrected p value threshold','+0','r',0.05,1,[0,1]);
iFWE=ceil((1-alph_FWE)*nPermReal);
if alph_FWE<1
C_MaxT=StMaxT(iFWE);
else
C_MaxT = 0;
end
u = C_MaxT;
otherwise %-NB: no adjustment
%%%% Now ask: Threshold statistic image or Uncorr P-value Image ?
%%%% If stats image, maybe do ST, no FDR-level of thresh;
%%%% If P-value image, no ST, can find FDR-level of thresh
% p for conjunctions is p of the conjunction SPM
%---------------------------------------------------------------
if bVarSm, str = 'pseudo t'; else, str = sprintf('t_{%d}',df); end
u = spm_input(['threshold (',str,')'],'+0','r',0.01,1,[0,Inf]);
alph_FWE = sum(MaxT > u -tol) / nPermReal;
end
else % Cluster-wise inference
if pU_ST_Ut==-1 % No threshold was set in snpm_ui.
%-Get primary threshold for STC analysis if requested
%-----------------------------------------------------------------------
% Save original ST_Ut
ST_Ut_0 = ST_Ut;
%-Threshold must be greater or equal to that (ST_Ut) used to collect
% suprathreshold data in snpm_cp
%-If a test level alpha has been set, then it there's no sense in having
% the threshold greater than C_MaxT, above which voxels are individually
% significant
tmp = 0;
if bVarSm
%-If using pseudo-statistics then can't use (uncorrected)
% upper tail p-values to specify primary threshold
while ~(tmp>=ST_Ut-tol)
tmp = spm_input(sprintf(...
'Clus-def thresh(pseudo t>%4.2f)',ST_Ut),'+0');
end
if (abs(tmp-ST_Ut)<=tol)
tmp=ST_Ut; % If tmp is very close to ST_Ut, set tmp equal to ST_Ut.
end
else
%-Statistic image is t with df degrees of freedom
p_ST_Ut = STalpha;
while ~( tmp>=ST_Ut-tol || (tmp>0 && tmp<=p_ST_Ut))
tmp = spm_input(sprintf(...
'Clus-def thresh(p<=%4.2fIt>=%4.2f)',p_ST_Ut,ST_Ut),'+0','r',ST_Ut,1);
end
clear p_ST_Ut
if (tmp < 1)
tmp = spm_invTcdf(1-tmp,df);
else
if (abs(tmp-ST_Ut)<=tol)
tmp=ST_Ut; % If tmp is very close to ST_Ut, set tmp equal to ST_Ut.
end
end
end
ST_Ut = tmp;
end
u=ST_Ut; % Flag use of a statistic-value threshold
str = 'FWE|Uncorr';
switch spm_input('Clusterwise: Use corrected thresh?','+1','b',str,[],1)
case 'FWE' % family-wise false positive rate
%---------------------------------------------------------------
alph_FWE = spm_input('FWE-Corrected p value threshold','+0','r',0.05,1,[0,1]);
iFWE=ceil((1-alph_FWE)*nPermReal);
case 'Uncorr' %Uncorrected cluster size threshold
str = 'ClusterSize|P-value';
switch spm_input('Define uncorrected','+1','b',str,[],1)
case 'ClusterSize'
C_STCS = spm_input('Uncorr cluster size threshold','+0','w',0,1);
case 'P-value'
alpha_ucp = spm_input('Uncorrected p value threshold','+0','r',0.01,1,[0,1]);
end
end
end
end
end
% Workflow for statistical Inference:
%
% Results for T or P
%
% 1) P
%
% case 'FDR': alph_FDR is set by spm_input; alpha_ucp is calculated from
% alph_FDR and is used as the threshold.
% case 'None': alpha_ucp is set by spm_input; alph_FDR is calculated from
% alpha_ucp. alpha_ucp is used as the threshold.
%
% 2) T
%
% a) Voxel-wise
% case 'FWE': alph_FWE is set by spm_input; u=C_MaxT; iFWE is
% calculated. u is used as the threshold.
% case 'None': u is set by spm_input; alph_FWE is calculated from u.
% u is used as the threshold.
%
% b) Cluster-wise
% If cluster defining threshold not set, ask for pU_ST_Ut
% case 'FWE': alph_FWE is set by spm_input; iFWE is calculated from
% alph_FWE. C_STCS is calculated from iFWE and is used
% as the threshold.
% case 'Uncorr':
% i) 'ClusterSize', C_STCS is defined directly and used as
% the threshold;
% ii) 'P-value', alpha_ucp is set by spm_input. C_STCS is
% calculated from alpha_ucp and used as the threshold.
%
%=======================================================================
%- C O M P U T A T I O N
%=======================================================================
set(Finter,'Pointer','Watch')
%-Calculate distribution of Maximum Suprathreshold Cluster size
%-Calculate critical Suprathreshold Cluster Size
%=======================================================================
if bSpatEx
fprintf('Working on spatial extent...\n');
%-Compute suprathreshold voxels - check there are some
%---------------------------------------------------------------
fprintf('\tComputing suprathreshold voxels...');
Q = find(SnPMt > ST_Ut);
SnPMt = SnPMt(Q);
XYZ = XYZ(:,Q);
if isempty(Q)
set(Finter,'Pointer','Arrow')
figure(Fgraph)
axis off
text(0,0.97,CWD,'Fontsize',16,'FontWeight','Bold');
str=sprintf('No voxels above threshold %4.2f\n',ST_Ut);
text(0,0.93,str);
fprintf(['WARNING: ' str])
if length(strmatch('FWEreport',Report))>0
ShowDist(MaxT,C_MaxT,alph_FWE,[],[],[],'max');
if ~BATCH
if spm_input('Review permutation distributions.',1,'bd',...
'Print & Continue|Continue',[1,0],1)
spm_print
end
spm_clf(Fgraph)
end
end
if length(strmatch('FDRreport',Report))>0
axis off
text(0,0.97,'Uncorrected P Permutation Distributions','Fontsize',16,'FontWeight',...
'Bold');
ShowDist(SnPMucp,alpha_ucp,alph_FDR,[],[],[],'uncor');
if ~BATCH
if spm_input('Review permutation distributions.',1,'bd',...
'Print|Done',[1,0],1)
spm_print
end
end
end
return
end
fprintf('done\n')
%-Load & condition statistics
%---------------------------------------------------------------
if pU_ST_Ut==-1 % No threshold was set in snpm_ui.
fprintf('\tLoading & conditioning SupraThreshold statistics...');
try
load(fullfile(CWD,'SnPM_ST'))
catch exception
if strcmp(exception.identifier, 'MATLAB:load:unableToReadMatFile')
warning('SnPM:SnPMSTFileNotLOaded', ...
['SnPM_ST file can not be loaded. Consider using' ...
' ''set cluster-forming threshold now (fast)'' option' ...
' in SnPM ''Specify''.']);
end
% Rethrow exception
throw(exception)
end
%-SnPM_ST stores columns of [x;y;z;abs(t);perm] with perm negative
% where the exceedence was t < -ST_Ut_0
%-Trim statistics according to threshold ST_Ut, if ST_Ut > ST_Ut_0
tmp = find(SnPM_ST(4,:)>ST_Ut);
SnPM_ST = SnPM_ST(:,tmp);
clear tmp;
SnPM_ST_Pos = SnPM_ST(:,SnPM_ST(5,:)>0);
SnPM_ST_Neg = SnPM_ST(:,SnPM_ST(5,:)<0);
SnPM_ST_Neg(5,:) = -SnPM_ST_Neg(5,:);
fprintf('done\n')
%-Calculate distribution of Maximum SupraThreshold Cluster size
%---------------------------------------------------------------
fprintf('\tComputing dist. of max SupraThreshold cluster size: ');
STCS = snpm_STcalc('init',nPerm);
fprintf('\nPerms left: ');
for i = nPerm:-1:1
if (rem(i,10)==0)
fprintf('\b\b\b\b%-4u',i)
drawnow
end
if STAT== 'F'
loop = 1;
else
loop = 1:2;
end
for isPos= loop %1 for positive; 2 for negative
if isPos==1
SnPM_ST = SnPM_ST_Pos;
else
SnPM_ST = SnPM_ST_Neg;
end
tQ = (SnPM_ST(5,:)==i);
if any(tQ)
%-Compute cluster labellings for this perm
Locs_mm = SnPM_ST(1:3,tQ);
Locs_mm (4,:) = 1;
Locs_vox = IMAT * Locs_mm;
% Sometimes Locs_vox are not exactly integers and this raises an
% error later in the code. Here check that the values are
% integers with respect to a level of absolute tolerance (~10^-14)
% and enforce Locs_vox to be integers.
% (As in snpm_cp)
diffWithRounded = max(abs(Locs_vox(:)-round(Locs_vox(:))));
tolerance = 10^-10;
if diffWithRounded > tolerance
Locs_vox_alter = MAT\Locs_mm;
diffWithRounded_alter = max(abs(Locs_vox_alter(:)-round(Locs_vox(:))));
error('SnPM:NonIntegerLocsvox', ['''Locs_vox'' must be integers (difference is ' num2str(diffWithRounded) ...
' or ' num2str(diffWithRounded_alter) ')']);
else
Locs_vox = round(Locs_vox);
end
STCS = snpm_STcalc('update',STCS, SnPM_ST(4,tQ),...
Locs_vox(1:3,:),isPos,i,ST_Ut,df);
end
if i==1
%-Save perm 1 stats for use later - [X;Y;Z;T;perm;STCno]
tmp = spm_clusters(Locs_vox(1:3,:));
if isPos==1
STCstats_Pos = [ SnPM_ST(:,tQ); tmp];
if bNeg==0
STCstats=STCstats_Pos;
end
else
STCstats_Neg = [ SnPM_ST(:,tQ); tmp];
if bNeg==1
STCstats=STCstats_Neg;
end
end
end
end
end
fprintf('\b\b\b\bdone\n');
if bhPerms %Double the STCS variables.
STCS = snpm_STcalc('double',STCS);
end
%-Get the stats from STCS structure that will be used later
STCS_MxK = STCS.MxK(:,bNeg+1);
STCS_K=cat(1,STCS.K{:,bNeg+1});
else % A threshold was set in snpm_ui.
%-Load & condition statistics
%---------------------------------------------------------------
fprintf('\tLoading SupraThreshold statistics...');
load(fullfile(CWD,'STCS'))
STCS_MxK = STCS.MxK(:,bNeg+1);
STCS_K=cat(1,STCS.K{:,bNeg+1});
if (bNeg==0)
load(fullfile(CWD,'SnPM_pp'))
else
load(fullfile(CWD,'SnPM_pp_Neg'))
STCstats = STCstats_Neg;
end
end
%-Compute critical SupraThreshold Cluster size
if isnan(C_STCS)
if ~isnan(alph_FWE)
% STCS_MxK: a vector of maximum cluster sizes of all the permutations.
% STCS_MxK = STCS.MxK(:,bNeg+1);
[StMaxSTCS, iStMaxSTCS] = sort(STCS_MxK);
if alph_FWE < 1
C_STCS = StMaxSTCS(iFWE);
else
C_STCS = 0;
end
elseif ~isnan(alpha_ucp)
% STCS_K: a vector of all the cluster sizes of all the
% permutations.
[StKSTCS, iStKSTCS] = sort(STCS_K);
if alpha_ucp < 1
iucp=ceil((1-alpha_ucp)*length(STCS_K));
C_STCS = StKSTCS(iucp);
else
C_STCS = 0;
end
end
end
%-Check XYZ for points > ST_Ut in perm 1 matches
% XYZ computed above for SnPMt > ST_Ut
%if pU_ST_Ut==-1
% if ~all(all( SnPM_ST(1:3,SnPM_ST(5,:)==1) == XYZ ))
% error('SnPM:InvalidSTXYZ', 'ST XYZ don''t match between STCS & thresh')
% end
%else
if ~all(all( STCstats(1:3,:) == XYZ ))
error('SnPM:InvalidSTXYZ', 'ST XYZ don''t match between STCS & thresh')
end
%end
end
%-Save some time consuming results
%-----------------------------------------------------------------------
if bSpatEx && pU_ST_Ut==-1
save SnPM_pp.mat STCstats_Pos
if STAT == 'T'
save SnPM_pp_Neg.mat STCstats_Neg
end
save STCS.mat STCS
end
%-Filter data at specified corrected p-value alpha
%=======================================================================
if bSpatEx
%-Analysing spatial extent