-
Notifications
You must be signed in to change notification settings - Fork 10
/
snpm_P_FDR.m
148 lines (138 loc) · 5.13 KB
/
snpm_P_FDR.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
function [P] = snpm_P_FDR(Z,df,STAT,n,Ps)
% Returns the corrected FDR P value
% FORMAT [P] = snpm_P_FDR(Z,df,STAT,n,Ps)
%
% Z - height (minimum of n statistics)
% If empty, find all FDR P values for Ps vector (in which case
% STAT is ignored, as Ps must be P-values already).
% df - [df{interest} df{error}]
% STAT - Statisical field
% 'Z' - Gaussian field
% 'T' - T - field
% 'X' - Chi squared field
% 'F' - F - field
% 'P' - P-values
% n - number of component SPMs in conjunction
% Ps - Vector of sorted (ascending) p-values in search volume
%
% P - corrected FDR P value
%
%___________________________________________________________________________
%
% The Benjamini & Hochberch (1995) False Discovery Rate (FDR) procedure
% finds a threshold u such that the expected FDR is at most q. snpm_P_FDR
% returns the smallest q such that Z>u.
%
% If abs(n) > 1, a P-value for a minimum of n values of the statistic
% is returned. If n>0, the P-value assesses the conjunction null
% hypothesis of one or more of the null hypotheses being true. If n<0,
% then the P-value assess the global null of all nulls being true.
%
% FDR Background
%
% For a given threshold on a statistic image, the False Discovery Rate
% is the proportion of suprathreshold voxels which are false positives.
% Recall that the thresholding of each voxel consists of a hypothesis
% test, where the null hypothesis is rejected if the statistic is larger
% than threshold. In this terminology, the FDR is the proportion of
% rejected tests where the null hypothesis is actually true.
%
% A FDR proceedure produces a threshold that controls the expected FDR
% at or below q. The FDR adjusted p-value for a voxel is the smallest q
% such that the voxel would be suprathreshold.
%
% In comparison, a traditional multiple comparisons proceedure
% (e.g. Bonferroni or random field methods) controls Familywise Error
% rate (FWER) at or below alpha. FWER is the *chance* of one or more
% false positives anywhere (whereas FDR is a *proportion* of false
% positives). A FWER adjusted p-value for a voxel is the smallest alpha
% such that the voxel would be suprathreshold.
%
% If there is truely no signal in the image anywhere, then a FDR
% proceedure controls FWER, just as Bonferroni and random field methods
% do. (Precisely, controlling E(FDR) yields weak control of FWE). If
% there is some signal in the image, a FDR method should be more powerful
% than a traditional method.
%
% For careful definition of FDR-adjusted p-values (and distinction between
% corrected and adjusted p-values) see Yekutieli & Benjamini (1999).
%
%
% References
%
% Benjamini & Hochberg (1995), "Controlling the False Discovery Rate: A
% Practical and Powerful Approach to Multiple Testing". J Royal Stat Soc,
% Ser B. 57:289-300.
%
% Benjamini & Yekutieli (2001), "The Control of the false discovery rate
% in multiple testing under dependency". To appear, Annals of Statistics.
% Available at http://www.math.tau.ac.il/~benja
%
% Yekutieli & Benjamini (1999). "Resampling-based false discovery rate
% controlling multiple test procedures for correlated test
% statistics". J of Statistical Planning and Inference, 82:171-196.
%_______________________________________________________________________
% Copyright (C) 2013 The University of Warwick
% Id: snpm_P_FDR.m SnPM13 2013/10/12
% Thomas Nichols
% Based on FIL spm_P_FDR.m 2.6 Thomas Nichols 04/07/02
if n>0
Cnj_n = 1; % Inf on Conjunction Null
else
Cnj_n = abs(n); % Inf on Global Null
end
if isempty(Z)
AllP = 1;
Z = Ps;
STAT = 'P';
else
AllP = 0;
end
% Set Benjamini & Yeuketeli cV for independence/PosRegDep case
%-----------------------------------------------------------------------
cV = 1;
S = length(Ps);
% Calculate p value of Z
%-----------------------------------------------------------------------
if STAT == 'Z'
PZ = (1 - spm_Ncdf(Z)).^Cnj_n;
elseif STAT == 'T'
PZ = (1 - spm_Tcdf(Z,df(2))).^Cnj_n;
elseif STAT == 'X'
PZ = (1 - spm_Xcdf(Z,df(2))).^Cnj_n;
elseif STAT == 'F'
PZ = (1-spm_Fcdf(Z,df)).^Cnj_n;
elseif STAT == 'P'
PZ = Z;
end
%-Calculate FDR p values
%-----------------------------------------------------------------------
% If Z is a value in the statistic image, then the adjusted p-value
% defined in Yekutieli & Benjamini (1999) (eqn 3) is obtained. If Z
% isn't a value in the image, then the adjusted p-value for the next
% smallest statistic value (next largest uncorrected p) is returned.
%-"Corrected" p-values
%-----------------------------------------------------------------------
Qs = Ps*S./(1:S)'*cV;
% -"Adjusted" p-values
%-----------------------------------------------------------------------
P = zeros(size(Z));
if AllP
P(S) = Qs(S);
for i = S-1:-1:1
P(i) = min(Qs(i),P(i+1));
end
else
for i = 1:length(Z)
% Find PZ(i) in Ps, or smallest Ps(j) s.t. Ps(j) >= PZ(i)
%---------------------------------------------------------------
I = min(find(Ps>=PZ(i)));
% -"Adjusted" p-values
%---------------------------------------------------------------
if isempty(I)
P(i) = 1;
else
P(i) = min(Qs(I:S));
end
end
end